text:
using Plots, Statistics, MLDatasets, LinearAlgebra, Clustering, Printf, Dates, Random, Debugger
path1=realpath(dirname(@__FILE__)*"/../..")
include(string(path1,"/src/ANNBN.jl"))

tol1=0.01;ANNBN.isigm1(tol1);ANNBN.isigm1(1.0-tol1)
i_train=1000;i_test=i_train;vars=2
xx_train=rand(i_train,2);yy_train=zeros(i_train)
i1=xx_train[:,1].<0.5;i2=xx_train[:,2].<0.5;i12=(i1+i2).==2;yy_train[i12].=1.0
mi=minimum(yy_train);yy_train.-=mi;ma=maximum(yy_train);yy_train./=ma;yy_train.*=1.0-2tol1;yy_train.+=tol1
scatter3d(xx_train[:,1],xx_train[:,2],yy_train[:,1])
xx_test=rand(i_test,2);yy_test=zeros(i_test)
i1=xx_test[:,1].<0.5;i2=xx_test[:,2].<0.5;i12=(i1+i2).==2;yy_test[i12].=1.0
yy_test.-=mi;yy_test./=ma;yy_test.*=1.0-2tol1;yy_test.+=tol1
scatter3d!(xx_test[:,1],xx_test[:,2],yy_test)

neurons=Int64(floor(i_train/(vars+1)))
# neurons=100
inds_all,n_per_part,items_per_neuron=ANNBN.___clustering(neurons,xx_train,200)
# inds_all,n_per_part,items_per_neuron=ANNBN.___clustering(neurons,yy_train,200)
# inds_all,n_per_part,items_per_neuron=ANNBN.___clustering(neurons,[xx_train yy_train],200)
# plot(yy_train[inds_all])
# minimum(items_per_neuron)
inds_all=(1:i_train)[:,1]
# inds_all=randperm(i_train)
neurons=Int64(i_train/10)
items_per_neuron=10*ones(Int64,neurons)   # (Int64(floor(i_train/(neurons))))
n_per_part=[0;cumsum(items_per_neuron)]

##### Train Sigmoid
a_all,a_layer1,layer1=ANNBN.train_layer_1_sigmoid_fast(neurons,vars,i_train,n_per_part,inds_all,xx_train,yy_train)
predl1=[layer1 ones(i_train)]*a_layer1
predl1_test,layer1_test=ANNBN.predict_new(a_all,a_layer1,xx_test,i_test,neurons)
maetr=mean(abs.(yy_train-predl1))
maete=mean(abs.(yy_test-predl1_test))
scatter3d!(xx_train[:,1],xx_train[:,2],predl1)
scatter3d!(xx_test[:,1],xx_test[:,2],predl1_test)
i1=predl1.<=0.5;predl1[i1].=tol1;i1=predl1.>0.5;predl1[i1].=1.0-tol1;i2=abs.(predl1.-yy_train).>0.01;maetr=sum(i2)/length(i2)
i1=predl1_test.<=0.5;predl1_test[i1].=tol1;i1=predl1_test.>0.5;predl1_test[i1].=1.0-tol1;i2=abs.(predl1_test.-yy_test).<0.01;
maete=1.0-sum(i2)/length(i2)

@time bb=[ones(i_train) xx_train]\ANNBN.isigm1.(yy_train)
predl1=ANNBN.sigm1.([ones(i_train) xx_train]*bb);
maetr=mean(abs.(yy_train-predl1))
# i1=predl1.<=0.5;predl1[i1].=0.0001;i1=predl1.>0.5;predl1[i1].=0.9998;i2=abs.(predl1.-yy_train).<0.01;maetr=sum(i2)/length(i2)
predl1_test=ANNBN.sigm1.([ones(i_test) xx_test]*bb);
maete=mean(abs.(yy_test-predl1_test))
# i1=predl1_test.<=0.5;predl1_test[i1].=0.0001;i1=predl1_test.>0.5;predl1_test[i1].=0.9998;i2=abs.(predl1_test.-yy_test).<0.01;maete=1.0-sum(i2)/length(i2)
scatter3d!(xx_train[:,1],xx_train[:,2],predl1)

# 2nd layer
vars_f=copy(neurons)
xx_train_f=copy(layer1)
neurons_f=copy(neurons)
inds_all_f,n_per_part_f,items_per_neuron_f=ANNBN.___clustering(neurons_f,xx_train_f,200)
# inds_all_f=copy(inds_all)
# n_per_part_f=copy(n_per_part)

@time a_all_f,a_layer1_f,layer1_f=ANNBN.train_layer_1_sigmoid(neurons_f,vars_f,i_train,n_per_part_f,inds_all_f,xx_train_f,yy_train);
predl1_f=[layer1_f ones(i_train)]*a_layer1_f
maetr=mean(abs.(yy_train-predl1_f))
maximum(abs.(yy_train-predl1_f))
# i1=predl1_f.<=0.5;predl1_f[i1].=0.0001;i1=predl1_f.>0.5;predl1_f[i1].=0.9998;i2=abs.(predl1_f.-yy_train).>0.01;maetr=sum(i2)/length(i2)

xx_test_f=copy(layer1_test)
@time predl1_test_f,layer1_test_f=ANNBN.predict_new(a_all_f,a_layer1_f,xx_test_f,i_test,neurons_f);
predl1_test_f=[layer1_test_f ones(i_test)]*a_layer1_f
maete=mean(abs.(yy_test-predl1_test_f))
maximum(abs.(yy_test-predl1_test_f))
# i1=predl1_test_f.<=0.5;predl1_test_f[i1].=0.0001;i1=predl1_test_f.>0.5;predl1_test_f[i1].=0.9998;i2=abs.(predl1_test_f.-yy_test).<0.01;maete=1.0-sum(i2)/length(i2)
# ir=randperm(size(xx_test,1))[1:50]; plot(yy_test[ir]); plot!(predl1_test[ir])
i1=predl1_f.<=0.5;predl1_f[i1].=0.1;i1=predl1_f.>0.5;predl1_f[i1].=0.9;i2=abs.(predl1_f.-yy_train).>0.01;maetr=sum(i2)/length(i2)
i1=predl1_test_f.<=0.5;predl1_test_f[i1].=0.1;i1=predl1_test_f.>0.5;predl1_test_f[i1].=0.9;i2=abs.(predl1_test_f.-yy_test).<0.01;maete=1.0-sum(i2)/length(i2)

# 3nd layer

for i=1:20
    global neurons_f,layer1_f,layer1_test_f
    vars_f=copy(neurons_f)
    xx_train_f=copy(layer1_f)
    neurons_f=copy(neurons_f)
    inds_all_f,n_per_part_f,items_per_neuron_f=ANNBN.___clustering(neurons_f,xx_train_f,200)

    @time a_all_f,a_layer1_f,layer1_f=ANNBN.train_layer_1_sigmoid(neurons_f,vars_f,i_train,n_per_part_f,inds_all_f,xx_train_f,yy_train);
    predl1_f=[layer1_f ones(i_train)]*a_layer1_f
    maetr=mean(abs.(yy_train-predl1_f))
    maximum(abs.(yy_train-predl1_f))

    xx_test_f=copy(layer1_test_f)
    @time predl1_test_f,layer1_test_f=ANNBN.predict_new(a_all_f,a_layer1_f,xx_test_f,i_test,neurons_f);
    # predl1_test_f=[layer1_test_f ones(i_test)]*a_layer1_f
    maete=mean(abs.(yy_test-predl1_test_f))
    maximum(abs.(yy_test-predl1_test_f))

    i1=predl1_f.<=0.5;predl1_f[i1].=0.1;i1=predl1_f.>0.5;predl1_f[i1].=0.9;i2=abs.(predl1_f.-yy_train).>0.01;maetr=sum(i2)/length(i2)
    i1=predl1_test_f.<=0.5;predl1_test_f[i1].=0.1;i1=predl1_test_f.>0.5;predl1_test_f[i1].=0.9;i2=abs.(predl1_test_f.-yy_test).<0.01;maete=1.0-sum(i2)/length(i2)
    println(i,"   >>>>>>>>>>>>>>>>>>>>>>maetr=",maetr,"    >>>>>>>>>>>>>>>>>>>>>>>>>    maete=",maete)
end
# ir=randperm(size(xx_test,1))[1:50]; plot(yy_test[ir]); plot!(predl1_test[ir])
# scatter3d!(xx_train[:,1],xx_train[:,2],predl1_f)
# scatter3d!(xx_test[:,1],xx_test[:,2],predl1_test_f)

using Plots, Statistics, MLDatasets, LinearAlgebra, Clustering, Printf, Dates, Random, MLDatasets, Plots, Printf

train_x, yy_train_all = MLDatasets.MNIST.traindata()
test_x,  yy_test_all  = MLDatasets.MNIST.testdata()
xx_train=MNIST.convert2features(MNIST.traintensor())'
xx_test=MNIST.convert2features(test_x)'
xx_train=convert(Array{Float64,2},xx_train)
xx_test=convert(Array{Float64,2},xx_test)
i_train=size(xx_train,1)
i_test=size(xx_test,1)
vars=size(xx_train,2)
num_to_test=3
yy_train=copy(yy_train_all)
i1=abs.(yy_train.-num_to_test).<0.01;i2=abs.(yy_train.-num_to_test).>0.01;yy_train[i1].=1;;yy_train[i2].=0;
yy_train=0.8*(convert(Array{Float64},yy_train));yy_train.+=0.1

yy_test=copy(yy_test_all)
i1=abs.(yy_test.-num_to_test).<0.01;i2=abs.(yy_test.-num_to_test).>0.01;yy_test[i1].=1;yy_test[i2].=0;
yy_test=0.8*(convert(Array{Float64},yy_test));yy_test.+=0.1

bb=[ones(i_train) xx_train]\ANNBN.isigm1.(yy_train)
bb=[-500;500;500]
predl1=ANNBN.sigm1.([ones(i_train) xx_train]*bb);
maetr=mean(abs.(yy_train-predl1))
# i1=predl1.<=0.5;predl1[i1].=0.0001;i1=predl1.>0.5;predl1[i1].=0.9998;i2=abs.(predl1.-yy_train).<0.01;maetr=sum(i2)/length(i2)
predl1_test=ANNBN.sigm1.([ones(i_test) xx_test]*bb);
maete=mean(abs.(yy_test-predl1_test))
# i1=predl1_test.<=0.5;predl1_test[i1].=0.0001;i1=predl1_test.>0.5;predl1_test[i1].=0.9998;i2=abs.(predl1_test.-yy_test).<0.01;maete=1.0-sum(i2)/length(i2)
scatter3d!(xx_train[:,1],xx_train[:,2],predl1)
scatter3d!(xx_test[:,1],xx_test[:,2],predl1_test)

# xx_train=collect(xx_train)[:,:];xx_test=collect(xx_test)[:,:]
cc1=100.0; @time a_all,a_layer1,layer1=ANNBN.train_layer_1_rbf(neurons,vars,i_train,n_per_part,inds_all,xx_train,yy_train,cc1);
predl1=[layer1 ones(i_train)]*a_layer1
maetr=mean(abs.(yy_train-predl1))
maximum(abs.(yy_train-predl1))
# i1=predl1.<=0.5;predl1[i1].=0.0001;i1=predl1.>0.5;predl1[i1].=0.9998;i2=abs.(predl1.-yy_train).<0.01;maetr=sum(i2)/length(i2)
@time predl1_test,layer1_test=ANNBN.predict_new_rbf(a_all,a_layer1,xx_test,i_test,neurons,n_per_part,inds_all,xx_train,vars,cc1);
maete=mean(abs.(yy_test-predl1_test))
maximum(abs.(yy_test-predl1_test))
# i1=predl1_test.<=0.5;predl1_test[i1].=0.0001;i1=predl1_test.>0.5;predl1_test[i1].=0.9998;i2=abs.(predl1_test.-yy_test).<0.01;maete=1.0-sum(i2)/length(i2)
scatter3d!(xx_train[:,1],xx_train[:,2],predl1)
scatter3d!(xx_test[:,1],xx_test[:,2],predl1_test)

using DecisionTree, Statistics
num_trees=100;
regr_3 = RandomForestClassifier(n_trees=num_trees);
tmp_xx_train=copy(xx_train);tmp_yy_train=copy(yy_train);
tmp_xx_test=copy(xx_test);tmp_yy_test=copy(yy_test);
DecisionTree.fit!(regr_3, tmp_xx_train, tmp_yy_train[:,1]);
predl1 = DecisionTree.predict(regr_3, tmp_xx_train);
predl1_test = DecisionTree.predict(regr_3, tmp_xx_test);
maetr=mean(abs.(yy_train-predl1))
maete=mean(abs.(yy_test-predl1_test))
i1=predl1.<=0.5;predl1[i1].=0.0001;i1=predl1.>0.5;predl1[i1].=0.9998;i2=abs.(predl1.-yy_train).<0.01;maetr=sum(i2)/length(i2)
i1=predl1_test.<=0.5;predl1_test[i1].=0.0001;i1=predl1_test.>0.5;predl1_test[i1].=0.9998;i2=abs.(predl1_test.-yy_test).<0.01;maete=1.0-sum(i2)/length(i2)
scatter3d!(xx_train[:,1],xx_train[:,2],predl1)
scatter3d!(xx_test[:,1],xx_test[:,2],predl1_test)

using PyCall,ScikitLearn
@sk_import ensemble: AdaBoostClassifier
sci_model = fit!(AdaBoostClassifier(n_estimators=100),xx_train,string.(yy_train))
predl1=parse.(Float64,predict(sci_model,xx_train))
predl1_test=parse.(Float64,predict(sci_model,xx_test))
maetr=mean(abs.(yy_train-predl1))
maete=mean(abs.(yy_test-predl1_test))
scatter3d!(xx_train[:,1],xx_train[:,2],predl1)
scatter3d!(xx_test[:,1],xx_test[:,2],predl1_test)

using XGBoost
num_round = 100
bst = xgboost(copy(xx_train), num_round, label = copy(yy_train), eta = 1, max_depth = 7)
predl1 = XGBoost.predict(bst, copy(xx_train))
predl1_test = XGBoost.predict(bst, copy(xx_test))
maetr=mean(abs.(yy_train-predl1))
maete=mean(abs.(yy_test-predl1_test))

i1=predl1.<=0.5;predl1[i1].=0.0001;i1=predl1.>0.5;predl1[i1].=0.9998;i2=abs.(predl1.-yy_train).<0.01;maetr=sum(i2)/length(i2)
i1=predl1_test.<=0.5;predl1_test[i1].=0.0001;i1=predl1_test.>0.5;predl1_test[i1].=0.9998;i2=abs.(predl1_test.-yy_test).<0.01;maete=1.0-sum(i2)/length(i2)

scatter3d!(xx_train[:,1],xx_train[:,2],predl1)
scatter3d!(xx_test[:,1],xx_test[:,2],predl1_test)


text:
@testset "Multidim" begin

using DiffEqNoiseProcess, Random, Statistics

W = WienerProcess(0.0,rand(4,4),rswm=RSWM(adaptivealg=:RSwM3))

Random.seed!(200)
dt = 0.2
calculate_step!(W,dt)

for i in 1:100
  reject_step!(W,0.01)
  accept_step!(W,0.2)
  accept_step!(W,0.1)
end

dt = dt/100

for t in dt:dt:W.t[end]-dt
  W(t)
end

W = WienerProcess!(0.0,rand(4,4),rand(4,4),rswm=RSWM(adaptivealg=:RSwM3))

Random.seed!(200)
dt = 0.2
calculate_step!(W,dt)

for i in 1:100
  reject_step!(W,0.01)
  accept_step!(W,0.2)
  accept_step!(W,0.1)
end

dt = dt/100
for t in dt:dt:W.t[end]-dt
  W(t)
end

end

ext: jl


text:
using SASLib, PooledArrays, BenchmarkTools

versioninfo()

opfn(size) = SASLib.ObjectPool{String}{UInt32}("", size)
pafn(size) = PooledArray{String,UInt32}(fill("", size))
rafn(size) = fill("", size)

size = 100000
dishes = ["food$i" for i in 1:size]

assignitems = (x::AbstractArray, y::AbstractArray) -> x[1:end] = y[1:end]
nsamples = 5000
nseconds = 600

println("** Regular Array **")
bra = @benchmarkable assignitems(x, y) setup=((x,y)=($rafn($size), $dishes)) samples=nsamples seconds=nseconds
display(run(bra))

println("\n\n** Pooled Array (Dict Pool) **")
bpa = @benchmarkable assignitems(x, y) setup=((x,y)=($pafn($size), $dishes)) samples=nsamples seconds=nseconds
display(run(bpa))

println("\n\n** Object Pool **")
bop = @benchmarkable assignitems(x, y) setup=((x,y)=($opfn($size), $dishes)) samples=nsamples seconds=nseconds
display(run(bop))

println()
println()


text:
const DATA_DIR = joinpath(MODULE_DIR, "..", "data")

const COERCE_ADULT = (
    :age => Continuous,
    :workclass => Multiclass,
    :fnlwgt => Continuous,
    :education => Multiclass,
    :education_num => Continuous,
    :marital_status => Multiclass,
    :occupation => Multiclass,
    :relationship => Multiclass,
    :race => Multiclass,
    :sex => Multiclass,
    :capital_gain => Continuous,
    :capital_loss => Continuous,
    :hours_per_week => Continuous,
    :native_country => Multiclass,
    :income_per_year => Multiclass,
)

"""
Checks whether the dataset is already present in data directory. Downloads it if not present.
"""
function ensure_download(url::String, file::String)
    cd(DATA_DIR) # This is to ensue that the dataset is not downloaded to /tmp instead of ./data
    fpath = joinpath(DATA_DIR, file)
    if !isfile(fpath)
        Downloads.download(url, file)
    end
end

"
Macro to read csv file of job data (data/jobs.csv) and convert columns to categorical.
Returns the tuple (X, y, ŷ)
"
macro load_toydata()
    quote
        fpath = joinpath(DATA_DIR, "jobs.csv")
        data = DataFrame(CSV.File(fpath, silencewarnings=true, delim=","); copycols = false)
        transform!(data, names(data)[1:6] .=> categorical, renamecols=false)
        (data[!, names(data)[1:4]], data[!, :Class], data[!, :Pred])
    end
end

"
Macro to create fairness Tensor for data/jobs.csv
The fairness tensor will be created on the basis of the column Job Type.
This column has 3 different values for job types.
"
macro load_toyfairtensor()
    quote
        X, y, ŷ = @load_toydata
        fair_tensor(ŷ, y, X[!, names(X)[4]])
    end
end

"""
Macro to load [COMPAS dataset](https://github.com/propublica/compas-analysis/blob/master/Compas%20Analysis.ipynb)
It is a reduced version of COMPAS Datset with 8 features and 6907 rows. The protected attributes are sex and race.
The available features are used to predict whether a criminal defendant will recidivate(reoffend).

Returns (X, y)
"""
macro load_compas()
    quote
        url = "https://raw.githubusercontent.com/propublica/compas-analysis/master/compas-scores-two-years.csv"
        fname = "compas-scores-two-years.csv"
        ensure_download(url, fname)

        fpath = joinpath(DATA_DIR, fname)
        data = DataFrame(CSV.File(fpath, silencewarnings=true, delim=","); copycols = false)
        data = data[!, ["sex", "age", "age_cat", "race", "c_charge_degree", "priors_count", "days_b_screening_arrest", "decile_score", "is_recid"]]
        dropmissing!(data, disallowmissing=true)
        coerce!(data, Textual => Multiclass)
        data.is_recid = Array(data.is_recid) #Temporary Fix for bug: ambiguous definition while coercing to multiclass for SentinelArrays
        coerce!(data, :is_recid=>Multiclass)
        y, X = unpack(data, ==(:is_recid), col -> true)

        X = data[!, ["sex", "age", "age_cat", "race", "c_charge_degree", "priors_count", "days_b_screening_arrest", "decile_score"]]
        y = data[!, "is_recid"]
        (X, y)
    end
end

"""
Macro to Load the [Adult dataset](https://archive.ics.uci.edu/ml/datasets/adult)
It has 14 features and 32561 rows. The protected attributes are race and sex.
This dataset is used to predict whether income exceeds 50K dollars per year.

Returns (X, y)
"""
macro load_adult()
    quote
        url = "https://archive.ics.uci.edu/ml/machine-learning-databases/adult/adult.data"
        fname = "adult.data"
        cols = ["age", "workclass", "fnlwgt", "education",
            "education_num", "marital_status", "occupation",
            "relationship", "race", "sex", "capital_gain",
            "capital_loss", "hours_per_week", "native_country",
            "income_per_year"
        ]
        ensure_download(url, fname)
        fpath = joinpath(DATA_DIR, fname)
        data = DataFrame(CSV.File(fpath, header=cols, silencewarnings=true, delim=", "); copycols = false)
        # Warning is silenced to supress warnings for lesser number of columns

        data = dropmissing(data, names(data))
        data.income_per_year = map(data.income_per_year) do η
            η == "<=50K" ? 0 : 1
        end

        coerce!(data, COERCE_ADULT...)
        coerce!(data, :income_per_year => Multiclass)
        y, X = unpack(data, ==(:income_per_year), col -> true)
        (X, y)
    end
end

"""
Load the full version of [German credit dataset](https://archive.ics.uci.edu/ml/datasets/Statlog+%28German+Credit+Data%29).
This dataset has 20 features and 1000 rows. The protected attributes are gender_status and age (>25 is priviledged)
Using the 20 features, it classifies the credit decision to a person as good or bad credit risks.

Returns (X, y)
"""
macro load_german()
    quote
        url = "https://archive.ics.uci.edu/ml/machine-learning-databases/statlog/german/german.data"
        fname = "german.data"
        cols = [
            "status", "duration", "credit_history",
            "purpose", "credit_amount", "savings", "employment",
            "installment_rate", "gender_status",
            "other_debtors", "residence_since", "property", "age",
            "installment_plans", "housing", "existing_credits",
            "skill_level", "people_liable", "telephone",
            "foreign_worker", "label"
        ]
        ensure_download(url, fname)
        fpath = joinpath(DATA_DIR, fname)
        df = DataFrame(CSV.File(fpath, header=cols, silencewarnings=true, delim=" "); copycols = false)

        gender_status = Dict(
            "A91" => "male_divorced_separated",
            "A92" => "female_divorced_separated_married",
            "A93" => "male_single",
            "A94" => "male_married_widowed",
            "A95" => "female_single", # There seems to be no row with this value in data
        )

        df.gender_status = map(df.gender_status) do η
            gender_status[η]
        end

        coerce!(df, Textual => Multiclass)
        df.label = Array(df.label)#Temporary Fix for bug: ambiguous definition while coercing to multiclass for SentinelArrays
        coerce!(df, :label => Multiclass)
        y, X = unpack(df, ==(:label), col -> true);
        (X, y)
    end
end

"""
The data is related with direct marketing campaigns of a Portuguese banking institution.
The marketing campaigns were based on phone calls. Often, more than one contact to the same client was required,
in order to access if the product (bank term deposit) would be ('yes') or not ('no') subscribed.
It has 20 features and 41188 rows. The protected attributes is marital.
"""
macro load_bank_marketing()
    quote
        url = "https://archive.ics.uci.edu/ml/machine-learning-databases/00222/bank-additional.zip"
        fname = "bank-marketing.csv"
        fpath = joinpath(DATA_DIR, fname)
        if !isfile(fpath)
            Downloads.download(url, "tempdataset.zip")
            zarchive = ZipFile.Reader("tempdataset.zip")
            zipfile = filter(x->x.name=="bank-additional/bank-additional-full.csv", zarchive.files)[1]
            df = DataFrame(CSV.File(read(zipfile), silencewarnings=true, delim=";"); copycols = false)
            CSV.write(fpath, df)
            close(zarchive)
            Base.Filesystem.rm("tempdataset.zip", recursive=true)
        end
        df = DataFrame(CSV.File(fpath, silencewarnings=true, delim=","); copycols = false)
        df.job = Array(df.job)#Temporary Fix for bug: ambiguous definition while coercing to multiclass for SentinelArrays
        coerce!(df, Textual => Multiclass)
        coerce!(df, :y => Multiclass)
        (y, X) = unpack(df, ==(:y), col->true)
        y = categorical(y)
        levels!(y, ["yes", "no"])
        (X, y)
    end
end

# TODO: Construct protected attribute from the attributes: blackpct, whitepct, asianpct
"""
The per capita violent crimes variable was calculated using population and
the sum of crime variables considered violent crimes in the United States: murder, rape, robbery, and assault.
It has 127 features and 1994 rows. The protected attributes are ....?
"""
macro load_communities_crime()
    quote
        url = "http://archive.ics.uci.edu/ml/machine-learning-databases/communities/communities.data"
        fname = "communities_crime.data"
        ensure_download(url, fname)
        fpath = joinpath(DATA_DIR, fname)
        df = DataFrame(CSV.File(fname, header=false, silencewarnings=true, delim=","); copycols = false)
        df = dropmissing(df, names(df))
        X = df[!, names(df)[1:127]]
        y = df[!, names(df)[128]] .> 0.7
        y = categorical(y)
        (X, y)
    end
end

"""
Student Performance Dataset. It has 395 rows and 30 features.
The target attribute corresponds to grade G1. The target tells whether the student gets grade >= 12.
The protected attribute is sex.
"""
macro load_student_performance()
    quote
        url= "https://archive.ics.uci.edu/ml/machine-learning-databases/00320/student.zip"
        fname = "student-performance.csv"
        fpath = joinpath(DATA_DIR, fname)
        if !isfile(fpath)
            Downloads.download(url, "tempdataset.zip")
            zarchive = ZipFile.Reader("tempdataset.zip")
            zipfile = filter(x->x.name=="student-mat.csv", zarchive.files)[1]
            df = DataFrame(CSV.File(read(zipfile), silencewarnings=true, delim=";"); copycols = false)
            CSV.write(fpath, df)
            close(zarchive)
            Base.Filesystem.rm("tempdataset.zip", recursive=true)
        end
        df = DataFrame(CSV.File(fpath, silencewarnings=true, delim=","); copycols = false)
        coerce!(df, Textual => Multiclass)
        X = df[!, names(df)[1:30]]
        y = df[!, names(df)[31]] .>= 12
        y = categorical(y)
        (X, y)
    end
end


text:
include("1.1 Binomial products/.jl")
include("1.2 Trinomial products/.jl")
include("1.3 Miscellaneous/.jl")

max_stars_repo_path: src/rules/1 Algebraic functions/.jl

max_issues_repo_path: src/rules/1 Algebraic functions/.jl

max_forks_repo_path: src/rules/1 Algebraic functions/.jl


text:
using TypedFun
using Test

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


text:
function similarity{V}(g::AbstractGraph{V})
    similarities = zeros(num_edges(g))
    for e in edges(g)
        e_idx = edge_index(e, g)
        u = source(e, g)
        v = target(e, g)
        u_nei = out_neighbors(u, g)
        v_nei = out_neighbors(v, g)
        uv_common_nei = intersect(u_nei, v_nei)
        uv_all_nei = union(u_nei, v_nei)
        similarities[e_idx] = (length(uv_common_nei) + 2) / length(uv_all_nei)
    end
    similarities
end

function triangle{V}(g::AbstractGraph{V})
    △ = zeros(Int, num_edges(g))
    for e in edges(g)
        e_idx = edge_index(e, g)
        u = source(e, g)
        v = target(e, g)
        u_nei = out_neighbors(u, g)
        v_nei = out_neighbors(v, g)
        uv_nei = intersect(u_nei, v_nei)
        △[e_idx] = length(uv_nei)
    end
    △
end

function quadrangle{V}(g::AbstractGraph{V})
    □ = zeros(Int, num_edges(g))
    for e in edges(g)
        e_idx = edge_index(e, g)
        u = source(e, g)
        v = target(e, g)
        u_nei = out_neighbors(u, g)
        v_nei = out_neighbors(v, g)
        u_unique_nei = setdiff(u_nei, v_nei)
        u_unique_nei = setdiff(u_unique_nei, v)
        v_unique_nei = setdiff(v_nei, u_nei)
        v_unique_nei = setdiff(v_unique_nei, u)
        for w in u_unique_nei
            w_nei = out_neighbors(w, g)
            vw_nei = intersect(v_unique_nei, w_nei)
            □[e_idx] += length(vw_nei)
        end
    end
    □
end

function triquadrsim{V}(g::AbstractGraph{V})
    △ = zeros(Int, num_edges(g))
    □ = zeros(Int, num_edges(g))
    similarities = zeros(num_edges(g))
    for e in edges(g)
        e_idx = edge_index(e, g)
        u = source(e, g)
        v = target(e, g)
        u_nei = out_neighbors(u, g)
        v_nei = out_neighbors(v, g)
        uv_nei = intersect(u_nei, v_nei)
        uv_all_nei = union(u_nei, v_nei)
        similarities[e_idx] = (length(uv_nei) + 2) / length(uv_all_nei)
        △[e_idx] = length(uv_nei)
        u_unique_nei = setdiff(u_nei, uv_nei)
        u_unique_nei = setdiff(u_unique_nei, v)
        v_unique_nei = setdiff(v_nei, uv_nei)
        v_unique_nei = setdiff(v_unique_nei, u)
        for w in u_unique_nei
            w_nei = out_neighbors(w, g)
            vw_nei = intersect(v_unique_nei, w_nei)
            □[e_idx] += length(vw_nei)
        end
    end
    △, □, similarities
end

max_stars_repo_name: afternone/LPA

max_issues_repo_name: afternone/LPA

max_forks_repo_name: afternone/LPA


text:
Dict(
        "hello" => "cava Latt"
       )


text:
function solve()
    x = parse(Int, readline())
    1 - x
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    Segment(p1, p2)

An oriented line segment with end points `p1`, `p2`.
The segment can be called as `s(t)` with `t` between
`0` and `1` to interpolate linearly between its endpoints.

See also [`Line`](@ref).
"""
struct Segment{Dim,T,V<:AbstractVector{Point{Dim,T}}} <: Polytope{1,Dim,T}
  vertices::V
end

isconvex(::Type{<:Segment}) = true

nvertices(::Type{<:Segment}) = 2
nvertices(s::Segment) = nvertices(typeof(s))

Base.minimum(s::Segment) = s.vertices[1]
Base.maximum(s::Segment) = s.vertices[2]
Base.extrema(s::Segment) = s.vertices[1], s.vertices[2]
measure(s::Segment) = norm(s.vertices[2] - s.vertices[1])
Base.length(s::Segment) = measure(s)

function (s::Segment)(t)
  if t < 0 || t > 1
    throw(DomainError(t, "s(t) is not defined for t outside [0, 1]."))
  end
  a, b = s.vertices
  a + t * (b - a)
end

function Base.in(p::Point{Dim,T}, s::Segment{Dim,T}) where {Dim,T}
  a, b = s.vertices
  if Dim == 2
    _0 = zero(T)
  elseif Dim == 3
    _0 = Vec{3,T}(0, 0, 0)
  else
    throw(ErrorException("not implemented"))
  end
  # given collinear points (a, b, p), the point p intersects
  # segment ab if and only if vectors satisfy 0 ≤ ap ⋅ ab ≤ ||ab||²
  iscollinear = isapprox((b - a) × (p - a), _0, atol=atol(T)^2)
  iscollinear && zero(T) ≤ (b - a) ⋅ (p - a) ≤ (b - a) ⋅ (b - a)
end


text:
const alphabet = ['a','b','c','d','e','f','g',
		  'h','i','j','k','l','m','n',
		  'o','p','q','r','s','t','u',
		  'v','w','x','y','z']
"""
Project a walk into the alphabet.
"""
function alphabet_projection(walk::Tuple)
    projected = Vector{Char}()
    map = Dict()
    map_idx = 1
    for w in walk
        if !haskey(map, w)
            map[w] = alphabet[map_idx]
            map_idx += 1
        end
        push!(projected, map[w])
    end
    return projected
end

"""
Accept a set of walks, assumed to be of the same length.
Project each walk into a motif. If no weights are provided,
each has weight 1.
"""
function count_motifs(walks::Vector{Tuple}, weights=ones(length(walks)))
    motifs = Dict()
    widx = 1
    for walk in walks
        p = alphabet_projection(walk)
        if !haskey(motifs, p)
            motifs[p] = 0
        end
        motifs[p] += weights[widx]
        widx += 1
    end
    return motifs
end


text:
using Test
using LinearAlgebra, SparseArrays
using LuxurySparse

@testset "IMatrix" begin
    include("IMatrix.jl")
end

@testset "PermMatrix" begin
    include("PermMatrix.jl")
end

@testset "kronecker" begin
    include("kronecker.jl")
end

@testset "linalg" begin
    include("linalg.jl")
end

@testset "staticize" begin
    include("staticize.jl")
end

@testset "broadcast" begin
    include("broadcast.jl")
end


text:
"""
Main module for `FourierFlows.jl` -- an ecosystem for solving partial differential equations
on periodic domains using Fourier-based pseudospectral methods."

# Exports
$(EXPORTS)
"""
module FourierFlows

export
  Device,
  CPU,
  GPU,
  ArrayType,

  cxtype,
  fltype,
  innereltype,

  AbstractVars,
  AbstractParams,

  AbstractGrid,
  ZeroDGrid,
  OneDGrid,
  TwoDGrid,
  ThreeDGrid,
  dealias!,
  gridpoints,

  Diagnostic,
  resize!,
  update!,
  increment!,

  Output,
  saveoutput,
  saveproblem,
  groupsize,
  savediagnostic,

  @zeros,
  @devzeros,
  @createarrays,
  @superzeros,
  devzeros,
  superzeros,
  supersize,

  AbstractTimeStepper,
  TimeStepper,
  ForwardEulerTimeStepper,
  FilteredForwardEulerTimeStepper,
  RK4TimeStepper,
  FilteredRK4TimeStepper,
  ETDRK4TimeStepper,
  FilteredETDRK4TimeStepper,
  AB3TimeStepper,
  FilteredAB3TimeStepper,
  stepforward!,
  step_until!

using
  FFTW,
  JLD2,
  Statistics,
  Interpolations,
  Reexport,
  DocStringExtensions

@reexport using FFTW: fft, ifft, rfft, irfft
@reexport using CUDA

import Base: resize!, getindex, setindex!, push!, append!, show, summary

using Base: fieldnames
using FFTW: fftfreq, rfftfreq

"Abstract supertype for grids."
abstract type AbstractGrid{T, A} end

"Abstract supertype for timesteppers."
abstract type AbstractTimeStepper{T} end

"Abstract supertype for parameters."
abstract type AbstractParams end

"Abstract supertype for variables."
abstract type AbstractVars end

"Abstract supertype for diagnostics."
abstract type AbstractDiagnostic end

"Abstract supertype for device."
abstract type Device end

"CPU device."
struct CPU <: Device end

"GPU device."
struct GPU <: Device end

# The main show
include("problem.jl")
include("domains.jl")
include("utils.jl")
include("diagnostics.jl")
include("output.jl")
include("timesteppers.jl")

# Physics
include("diffusion.jl")

# CUDA functionality
if CUDA.has_cuda()
  include("CuFourierFlows.jl")
end

function __init__()
  threads = Threads.nthreads()
  if threads > 1
    @info "FourierFlows will use $threads threads"
    FFTW.set_num_threads(threads)
  end

  if CUDA.has_cuda()
    @debug "CUDA-enabled GPU(s) detected:"
    for (gpu, dev) in enumerate(CUDA.devices())
      @debug "$dev: $(CUDA.name(dev))"
    end

    CUDA.allowscalar(false)
  end
end

function show(io::IO, vars::AbstractVars)
  names = propertynames(vars)
  showstring = ""
  for name in names[1:end-1]
    field = getproperty(vars, name)
    showstring = string(showstring, "  ├───── variable: " * string(name) * " -> ", summary(field), "\n")
  end
  name = names[end]
  field = getproperty(vars, name)
  showstring = string(showstring, "  └───── variable: " * string(name) * " -> ", summary(field), "\n")

  return print(io, "Variables\n", showstring)
end

summary(::Function) = "Function"

function show(io::IO, params::AbstractParams)
  names = propertynames(params)
  showstring = ""
  for name in names[1:end-1]
    field = getproperty(params, name)
    showstring = string(showstring, "  ├───── parameter: " * string(name) * " -> ", summary(field), "\n")
  end
  name = names[end]
  field = getproperty(params, name)
  showstring = string(showstring, "  └───── parameter: " * string(name) * " -> ", summary(field), "\n")

  return print(io, "Parameters\n", showstring)
end

end # module


text:
const content =
    vbox(
        link("content.jl",title(1, "Content")),
        "Rendering everyday Julia objects",
        "Markdown",
        "Font",
        "TeX/mathematics",
        "Images",
        "Plots",
        "Vector graphics",
    )

const layout =
    vbox(
        link("layout.jl", title(1, "Layout")),
        "Width and Height",
        md"`hbox` and `vbox`",
        "Padding",
        "Absolute positioning",
        "Packing space",
        "Fill color",
        "Borders and Colors",
    )

const widgets_toc =
    vbox(
        title(1, "Widgets"),
        "Basic input widgets",
        "Code mirror",
        "Camera widget",
    )

const interaction =
    vbox(
        title(1, "Interaction"),
        "Signals",
        "A counter",
        "A series of counters",
        "Animations",
        "Combining signals",
    )

function main(window)
    vbox(
        title(2, "Examples"),
        hbox(content, hskip(2em), layout),
        hbox(interaction, hskip(2em), widgets_toc),
    ) |> Escher.pad(2em)

end


text:
@testset "Moments" begin
    f = zeros(5,5,9)
    f1D = zeros(30,3)
    height = zeros(5,5)
    hei = zeros(30)
    velx = zeros(5,5)
    vely = zeros(5,5)
    vel = zeros(30)
    @testset "No velocity" begin
        f[:,:,1] .= 1.0
        Swalbe.moments!(height, velx, vely, f)
        @test all(height .== 1.0)
        f1D[:,1] .= 1.0
        Swalbe.moments!(hei, vel, f1D)
        @test all(hei .== 1.0)
    end
    @testset "Velocity" begin
        # 2D
        f[:,:,1] .= 1.0
        f[:,:,2] .= 0.1
        Swalbe.moments!(height, velx, vely, f)
        @test all(height .== 1.1)
        @test all(velx .== 0.1/1.1)
        # 1D
        f1D[:,1] .= 1.0
        f1D[:,2] .= 0.1
        Swalbe.moments!(hei, vel, f1D)
        @test all(hei .== 1.1)
        @test all(vel .== 0.1/1.1)
        # 2D
        f[:,:,3] .= 0.2
        Swalbe.moments!(height, velx, vely, f)
        @test all(height .== 1.3)
        @test all(velx .== 0.1/1.3)
        @test all(vely .== 0.2/1.3)
        # 1D
        f1D[:,3] .= 0.2
        Swalbe.moments!(hei, vel, f1D)
        @test all(hei .== 1.3)
        @test all(vel .== -0.1/1.3)
        # 2D
        f[:,:,4] .= -0.2
        Swalbe.moments!(height, velx, vely, f)
        @test all(height .== 1.1)
        @test all(velx .≈ 0.3/1.1)
        @test all(vely .== 0.2/1.1)
        f[:,:,5] .= -0.1
        Swalbe.moments!(height, velx, vely, f)
        @test all(height .== 1.0)
        @test all(velx .≈ 0.3/1.0)
        @test all(vely .≈ 0.3/1.0)
        f[:,:,6] .= 0.1
        Swalbe.moments!(height, velx, vely, f)
        @test all(height .== 1.1)
        @test all(velx .≈ 0.4/1.1)
        @test all(vely .≈ 0.4/1.1)
    end
end


text:
using SignatureGB
using Singular
P, (x1, x2, x3, x4, x5, x6, x7, x8, x9, x10)  = Singular.PolynomialRing(Fp(65521), ["x1", "x2", "x3", "x4", "x5", "x6", "x7", "x8", "x9", "x10"]);
gb = SignatureGB.f5([x1, x2, x3, x4, x5, x6, x7, x8, x9, x10]);F = [4761*x1^2*x3^2+49*x1^2*x5^2-140*x1^2*x5*x7+4*x1^2*x6^2+100*x1^2*x7^2-200*x1*x10^2*x6+10350*x1*x10*x3^2+770*x1*x10*x4*x5-1100*x1*x10*x4*x7+490*x1*x2*x4*x5-700*x1*x2*x4*x7+9246*x1*x3^2*x7+2208*x1*x3*x4^2-826*x1*x3*x5^2+1180*x1*x3*x5*x7-184*x1*x3*x6*x7-12834*x1*x3*x6*x9-1932*x1*x3*x9^2-98*x1*x5*x6*x9+140*x1*x6*x7*x9-16*x1*x6*x9^2+11909*x10^4+18236*x10^3*x2+1552*x10^3*x4+16878*x10^3*x6+8836*x10^2*x2^2+1504*x10^2*x2*x4+16356*x10^2*x2*x6+23329*x10^2*x3^2-10520*x10^2*x3*x7+784*x10^2*x3*x8+3960*x10^2*x3*x9+3089*x10^2*x4^2+1392*x10^2*x4*x6+7569*x10^2*x6^2+2206*x10^2*x7^2-3696*x10^2*x7*x9+16*x10^2*x8^2+884*x10^2*x9^2+13500*x10*x2*x3*x5-4900*x10*x2*x3*x6+3850*x10*x2*x4^2-12600*x10*x2*x5*x7+3300*x10*x2*x5*x9-200*x10*x2*x6*x8-4700*x10*x2*x7^2+1034*x10*x3^2*x7+2400*x10*x3*x4^2-6490*x10*x3*x4*x5+2160*x10*x3*x4*x9+2156*x10*x3*x5*x9-13950*x10*x3*x6*x9-368*x10*x3*x7*x8+17280*x10*x3*x8^2-2100*x10*x3*x9^2-770*x10*x4*x6*x9-400*x10*x4*x7^2-2016*x10*x4*x7*x9+528*x10*x4*x9^2+88*x10*x5*x8*x9-4350*x10*x6*x7^2-16128*x10*x7*x8^2+4224*x10*x8^2*x9+1225*x2^2*x4^2+5625*x2^2*x5^2+625*x2^2*x6^2-4130*x2*x3*x4*x5+2300*x2*x3*x6*x7+1800*x2*x4*x5*x9-490*x2*x4*x6*x9-550*x2*x5*x6*x9+14400*x2*x5*x8^2+3706*x3^2*x5^2+8721*x3^2*x7^2+2144*x3*x4^2*x7+360*x3*x4*x5*x9-690*x3*x5^3+826*x3*x5*x6*x9+1980*x3*x5*x7^2-1012*x3*x5*x7*x9-12462*x3*x6*x7*x9-1508*x3*x7*x9^2+256*x4^4-2976*x4^2*x6*x9-160*x4^2*x9^2-552*x4*x5^2*x9+1584*x4*x7^2*x9+2304*x4*x8^2*x9+529*x5^4-3036*x5^2*x7^2+121*x5^2*x9^2+8698*x6^2*x9^2+2604*x6*x9^3+4981*x7^4+9216*x8^4+212*x9^4+1164*x10^2*x6+1128*x10*x2*x6-16072*x10*x3*x4+96*x10*x4*x6-656*x10*x4*x8+1044*x10*x6^2+4100*x2*x4*x6+7544*x3*x4*x7-390*x3*x5*x6-1804*x4*x5*x9-312*x4*x6*x9+598*x5^2*x6-2016*x6*x7^2+176*x1*x6-4400*x10^2+1470*x3*x5-4048*x3*x7+6724*x4^2+1176*x4*x9-2254*x5^2+205*x6^2+6468*x7^2-352*x9^2-1274*x6+4337, 490*x1*x4*x5-700*x1*x4*x7+18236*x10^3+17672*x10^2*x2+1504*x10^2*x4+16356*x10^2*x6+13500*x10*x3*x5-4900*x10*x3*x6+3850*x10*x4^2-12600*x10*x5*x7+3300*x10*x5*x9-200*x10*x6*x8-4700*x10*x7^2+2450*x2*x4^2+11250*x2*x5^2+1250*x2*x6^2-4130*x3*x4*x5+2300*x3*x6*x7+1800*x4*x5*x9-490*x4*x6*x9-550*x5*x6*x9+14400*x5*x8^2+1128*x10*x6+4100*x4*x6, 9522*x1^2*x3+20700*x1*x10*x3+18492*x1*x3*x7+2208*x1*x4^2-826*x1*x5^2+1180*x1*x5*x7-184*x1*x6*x7-12834*x1*x6*x9-1932*x1*x9^2+46658*x10^2*x3-10520*x10^2*x7+784*x10^2*x8+3960*x10^2*x9+13500*x10*x2*x5-4900*x10*x2*x6+2068*x10*x3*x7+2400*x10*x4^2-6490*x10*x4*x5+2160*x10*x4*x9+2156*x10*x5*x9-13950*x10*x6*x9-368*x10*x7*x8+17280*x10*x8^2-2100*x10*x9^2-4130*x2*x4*x5+2300*x2*x6*x7+7412*x3*x5^2+17442*x3*x7^2+2144*x4^2*x7+360*x4*x5*x9-690*x5^3+826*x5*x6*x9+1980*x5*x7^2-1012*x5*x7*x9-12462*x6*x7*x9-1508*x7*x9^2-16072*x10*x4+7544*x4*x7-390*x5*x6+1470*x5-4048*x7, 770*x1*x10*x5-1100*x1*x10*x7+490*x1*x2*x5-700*x1*x2*x7+4416*x1*x3*x4+1552*x10^3+1504*x10^2*x2+6178*x10^2*x4+1392*x10^2*x6+7700*x10*x2*x4+4800*x10*x3*x4-6490*x10*x3*x5+2160*x10*x3*x9-770*x10*x6*x9-400*x10*x7^2-2016*x10*x7*x9+528*x10*x9^2+2450*x2^2*x4-4130*x2*x3*x5+1800*x2*x5*x9-490*x2*x6*x9+4288*x3*x4*x7+360*x3*x5*x9+1024*x4^3-5952*x4*x6*x9-320*x4*x9^2-552*x5^2*x9+1584*x7^2*x9+2304*x8^2*x9-16072*x10*x3+96*x10*x6-656*x10*x8+4100*x2*x6+7544*x3*x7-1804*x5*x9-312*x6*x9+13448*x4+1176*x9, 98*x1^2*x5-140*x1^2*x7+770*x1*x10*x4+490*x1*x2*x4-1652*x1*x3*x5+1180*x1*x3*x7-98*x1*x6*x9+13500*x10*x2*x3-12600*x10*x2*x7+3300*x10*x2*x9-6490*x10*x3*x4+2156*x10*x3*x9+88*x10*x8*x9+11250*x2^2*x5-4130*x2*x3*x4+1800*x2*x4*x9-550*x2*x6*x9+14400*x2*x8^2+7412*x3^2*x5+360*x3*x4*x9-2070*x3*x5^2+826*x3*x6*x9+1980*x3*x7^2-1012*x3*x7*x9-1104*x4*x5*x9+2116*x5^3-6072*x5*x7^2+242*x5*x9^2-390*x3*x6-1804*x4*x9+1196*x5*x6+1470*x3-4508*x5, 8*x1^2*x6-200*x1*x10^2-184*x1*x3*x7-12834*x1*x3*x9-98*x1*x5*x9+140*x1*x7*x9-16*x1*x9^2+16878*x10^3+16356*x10^2*x2+1392*x10^2*x4+15138*x10^2*x6-4900*x10*x2*x3-200*x10*x2*x8-13950*x10*x3*x9-770*x10*x4*x9-4350*x10*x7^2+1250*x2^2*x6+2300*x2*x3*x7-490*x2*x4*x9-550*x2*x5*x9+826*x3*x5*x9-12462*x3*x7*x9-2976*x4^2*x9+17396*x6*x9^2+2604*x9^3+1164*x10^2+1128*x10*x2+96*x10*x4+2088*x10*x6+4100*x2*x4-390*x3*x5-312*x4*x9+598*x5^2-2016*x7^2+176*x1+410*x6-1274, -140*x1^2*x5+200*x1^2*x7-1100*x1*x10*x4-700*x1*x2*x4+9246*x1*x3^2+1180*x1*x3*x5-184*x1*x3*x6+140*x1*x6*x9-10520*x10^2*x3+4412*x10^2*x7-3696*x10^2*x9-12600*x10*x2*x5-9400*x10*x2*x7+1034*x10*x3^2-368*x10*x3*x8-800*x10*x4*x7-2016*x10*x4*x9-8700*x10*x6*x7-16128*x10*x8^2+2300*x2*x3*x6+17442*x3^2*x7+2144*x3*x4^2+3960*x3*x5*x7-1012*x3*x5*x9-12462*x3*x6*x9-1508*x3*x9^2+3168*x4*x7*x9-6072*x5^2*x7+19924*x7^3+7544*x3*x4-4032*x6*x7-4048*x3+12936*x7, 784*x10^2*x3+32*x10^2*x8-200*x10*x2*x6-368*x10*x3*x7+34560*x10*x3*x8+88*x10*x5*x9-32256*x10*x7*x8+8448*x10*x8*x9+28800*x2*x5*x8+4608*x4*x8*x9+36864*x8^3-656*x10*x4, -12834*x1*x3*x6-3864*x1*x3*x9-98*x1*x5*x6+140*x1*x6*x7-32*x1*x6*x9+3960*x10^2*x3-3696*x10^2*x7+1768*x10^2*x9+3300*x10*x2*x5+2160*x10*x3*x4+2156*x10*x3*x5-13950*x10*x3*x6-4200*x10*x3*x9-770*x10*x4*x6-2016*x10*x4*x7+1056*x10*x4*x9+88*x10*x5*x8+4224*x10*x8^2+1800*x2*x4*x5-490*x2*x4*x6-550*x2*x5*x6+360*x3*x4*x5+826*x3*x5*x6-1012*x3*x5*x7-12462*x3*x6*x7-3016*x3*x7*x9-2976*x4^2*x6-320*x4^2*x9-552*x4*x5^2+1584*x4*x7^2+2304*x4*x8^2+242*x5^2*x9+17396*x6^2*x9+7812*x6*x9^2+848*x9^3-1804*x4*x5-312*x4*x6+1176*x4-704*x9, -400*x1*x10*x6+10350*x1*x3^2+770*x1*x4*x5-1100*x1*x4*x7+47636*x10^3+54708*x10^2*x2+4656*x10^2*x4+50634*x10^2*x6+17672*x10*x2^2+3008*x10*x2*x4+32712*x10*x2*x6+46658*x10*x3^2-21040*x10*x3*x7+1568*x10*x3*x8+7920*x10*x3*x9+6178*x10*x4^2+2784*x10*x4*x6+15138*x10*x6^2+4412*x10*x7^2-7392*x10*x7*x9+32*x10*x8^2+1768*x10*x9^2+13500*x2*x3*x5-4900*x2*x3*x6+3850*x2*x4^2-12600*x2*x5*x7+3300*x2*x5*x9-200*x2*x6*x8-4700*x2*x7^2+1034*x3^2*x7+2400*x3*x4^2-6490*x3*x4*x5+2160*x3*x4*x9+2156*x3*x5*x9-13950*x3*x6*x9-368*x3*x7*x8+17280*x3*x8^2-2100*x3*x9^2-770*x4*x6*x9-400*x4*x7^2-2016*x4*x7*x9+528*x4*x9^2+88*x5*x8*x9-4350*x6*x7^2-16128*x7*x8^2+4224*x8^2*x9+2328*x10*x6+1128*x2*x6-16072*x3*x4+96*x4*x6-656*x4*x8+1044*x6^2-8800*x10];
sort!(F, by = p -> leading_monomial(p))
println(@elapsed SignatureGB.f5(F, verbose=true));exit()


text:
# RUN: julia -e "import Brutus; Brutus.lit(:emit_lowered)" --startup-file=no %s 2>&1 | FileCheck %s

function gauss(N)
    acc = 0
    for i in 1:N
        acc += i
    end
    return acc
end
emit(gauss, Int64)

# CHECK: Core.MethodMatch(Tuple{typeof(Main.Main.gauss), Int64}, svec(), gauss(N) in Main.Main at /home/mccoy/Dev/brutus/test/Codegen/lower/gauss.jl:3, true)after translating to MLIR in JLIR dialect:module  {
# CHECK-NEXT:   func nested @"Tuple{typeof(Main.gauss), Int64}"(%arg0: !jlir<"typeof(Main.gauss)">, %arg1: !jlir.Int64) -> !jlir.Int64 attributes {llvm.emit_c_interface} {
# CHECK-NEXT:     "jlir.goto"()[^bb1] : () -> ()
# CHECK-NEXT:   ^bb1:  // pred: ^bb0
# CHECK-NEXT:     %0 = "jlir.constant"() {value = #jlir<"#<intrinsic #29 sle_int>">} : () -> !jlir<"typeof(Core.IntrinsicFunction)">
# CHECK-NEXT:     %1 = "jlir.constant"() {value = #jlir<"1">} : () -> !jlir.Int64
# CHECK-NEXT:     %2 = "jlir.call"(%0, %1, %arg1) : (!jlir<"typeof(Core.IntrinsicFunction)">, !jlir.Int64, !jlir.Int64) -> !jlir.Bool
# CHECK-NEXT:     %3 = "jlir.constant"() {value = #jlir.ifelse} : () -> !jlir<"typeof(ifelse)">
# CHECK-NEXT:     %4 = "jlir.constant"() {value = #jlir<"0">} : () -> !jlir.Int64
# CHECK-NEXT:     %5 = "jlir.call"(%3, %2, %arg1, %4) : (!jlir<"typeof(ifelse)">, !jlir.Bool, !jlir.Int64, !jlir.Int64) -> !jlir.Int64
# CHECK-NEXT:     %6 = "jlir.constant"() {value = #jlir<"#<intrinsic #27 slt_int>">} : () -> !jlir<"typeof(Core.IntrinsicFunction)">
# CHECK-NEXT:     %7 = "jlir.constant"() {value = #jlir<"1">} : () -> !jlir.Int64
# CHECK-NEXT:     %8 = "jlir.call"(%6, %5, %7) : (!jlir<"typeof(Core.IntrinsicFunction)">, !jlir.Int64, !jlir.Int64) -> !jlir.Bool
# CHECK-NEXT:     "jlir.gotoifnot"(%8)[^bb3, ^bb2] {operand_segment_sizes = dense<[1, 0, 0]> : vector<3xi32>} : (!jlir.Bool) -> ()
# CHECK-NEXT:   ^bb2:  // pred: ^bb1
# CHECK-NEXT:     %9 = "jlir.constant"() {value = #jlir.true} : () -> !jlir.Bool
# CHECK-NEXT:     %10 = "jlir.undef"() : () -> !jlir.Int64
# CHECK-NEXT:     %11 = "jlir.undef"() : () -> !jlir.Int64
# CHECK-NEXT:     "jlir.goto"(%9, %10, %11)[^bb4] : (!jlir.Bool, !jlir.Int64, !jlir.Int64) -> ()
# CHECK-NEXT:   ^bb3:  // pred: ^bb1
# CHECK-NEXT:     %12 = "jlir.constant"() {value = #jlir.false} : () -> !jlir.Bool
# CHECK-NEXT:     %13 = "jlir.constant"() {value = #jlir<"1">} : () -> !jlir.Int64
# CHECK-NEXT:     %14 = "jlir.constant"() {value = #jlir<"1">} : () -> !jlir.Int64
# CHECK-NEXT:     "jlir.goto"(%12, %13, %14)[^bb4] : (!jlir.Bool, !jlir.Int64, !jlir.Int64) -> ()
# CHECK-NEXT:   ^bb4(%15: !jlir.Bool, %16: !jlir.Int64, %17: !jlir.Int64):  // 2 preds: ^bb2, ^bb3
# CHECK-NEXT:     %18 = "jlir.constant"() {value = #jlir<"#<intrinsic #44 not_int>">} : () -> !jlir<"typeof(Core.IntrinsicFunction)">
# CHECK-NEXT:     %19 = "jlir.call"(%18, %15) : (!jlir<"typeof(Core.IntrinsicFunction)">, !jlir.Bool) -> !jlir.Bool
# CHECK-NEXT:     %20 = "jlir.constant"() {value = #jlir<"0">} : () -> !jlir.Int64
# CHECK-NEXT:     %21 = "jlir.constant"() {value = #jlir<"0">} : () -> !jlir.Int64
# CHECK-NEXT:     "jlir.gotoifnot"(%19, %21, %16, %17, %20)[^bb10, ^bb5] {operand_segment_sizes = dense<[1, 1, 3]> : vector<3xi32>} : (!jlir.Bool, !jlir.Int64, !jlir.Int64, !jlir.Int64, !jlir.Int64) -> ()
# CHECK-NEXT:   ^bb5(%22: !jlir.Int64, %23: !jlir.Int64, %24: !jlir.Int64):  // 2 preds: ^bb4, ^bb9
# CHECK-NEXT:     %25 = "jlir.constant"() {value = #jlir<"#<intrinsic #2 add_int>">} : () -> !jlir<"typeof(Core.IntrinsicFunction)">
# CHECK-NEXT:     %26 = "jlir.call"(%25, %24, %22) : (!jlir<"typeof(Core.IntrinsicFunction)">, !jlir.Int64, !jlir.Int64) -> !jlir.Int64
# CHECK-NEXT:     %27 = "jlir.constant"() {value = #jlir<"===">} : () -> !jlir<"typeof(===)">
# CHECK-NEXT:     %28 = "jlir.call"(%27, %23, %5) : (!jlir<"typeof(===)">, !jlir.Int64, !jlir.Int64) -> !jlir.Bool
# CHECK-NEXT:     "jlir.gotoifnot"(%28)[^bb7, ^bb6] {operand_segment_sizes = dense<[1, 0, 0]> : vector<3xi32>} : (!jlir.Bool) -> ()
# CHECK-NEXT:   ^bb6:  // pred: ^bb5
# CHECK-NEXT:     %29 = "jlir.undef"() : () -> !jlir.Int64
# CHECK-NEXT:     %30 = "jlir.undef"() : () -> !jlir.Int64
# CHECK-NEXT:     %31 = "jlir.constant"() {value = #jlir.true} : () -> !jlir.Bool
# CHECK-NEXT:     "jlir.goto"(%29, %30, %31)[^bb8] : (!jlir.Int64, !jlir.Int64, !jlir.Bool) -> ()
# CHECK-NEXT:   ^bb7:  // pred: ^bb5
# CHECK-NEXT:     %32 = "jlir.constant"() {value = #jlir<"#<intrinsic #2 add_int>">} : () -> !jlir<"typeof(Core.IntrinsicFunction)">
# CHECK-NEXT:     %33 = "jlir.constant"() {value = #jlir<"1">} : () -> !jlir.Int64
# CHECK-NEXT:     %34 = "jlir.call"(%32, %23, %33) : (!jlir<"typeof(Core.IntrinsicFunction)">, !jlir.Int64, !jlir.Int64) -> !jlir.Int64
# CHECK-NEXT:     %35 = "jlir.constant"() {value = #jlir.false} : () -> !jlir.Bool
# CHECK-NEXT:     "jlir.goto"(%34, %34, %35)[^bb8] : (!jlir.Int64, !jlir.Int64, !jlir.Bool) -> ()
# CHECK-NEXT:   ^bb8(%36: !jlir.Int64, %37: !jlir.Int64, %38: !jlir.Bool):  // 2 preds: ^bb6, ^bb7
# CHECK-NEXT:     %39 = "jlir.constant"() {value = #jlir<"#<intrinsic #44 not_int>">} : () -> !jlir<"typeof(Core.IntrinsicFunction)">
# CHECK-NEXT:     %40 = "jlir.call"(%39, %38) : (!jlir<"typeof(Core.IntrinsicFunction)">, !jlir.Bool) -> !jlir.Bool
# CHECK-NEXT:     "jlir.gotoifnot"(%40, %26)[^bb10, ^bb9] {operand_segment_sizes = dense<[1, 1, 0]> : vector<3xi32>} : (!jlir.Bool, !jlir.Int64) -> ()
# CHECK-NEXT:   ^bb9:  // pred: ^bb8
# CHECK-NEXT:     "jlir.goto"(%36, %37, %26)[^bb5] : (!jlir.Int64, !jlir.Int64, !jlir.Int64) -> ()
# CHECK-NEXT:   ^bb10(%41: !jlir.Int64):  // 2 preds: ^bb4, ^bb8
# CHECK-NEXT:     "jlir.return"(%41) : (!jlir.Int64) -> ()
# CHECK-NEXT:   }
# CHECK-NEXT: }

# CHECK: module  {
# CHECK-NEXT:   func nested @"Tuple{typeof(Main.gauss), Int64}"(%arg0: !jlir<"typeof(Main.gauss)">, %arg1: i64) -> i64 attributes {llvm.emit_c_interface} {
# CHECK-NEXT:     %c1_i64 = constant 1 : i64
# CHECK-NEXT:     %c0_i64 = constant 0 : i64
# CHECK-NEXT:     %false = constant false
# CHECK-NEXT:     %true = constant true
# CHECK-NEXT:     %0 = cmpi sle, %c1_i64, %arg1 : i64
# CHECK-NEXT:     %1 = select %0, %arg1, %c0_i64 : i64
# CHECK-NEXT:     %2 = "jlir.convertstd"(%1) : (i64) -> !jlir.Int64
# CHECK-NEXT:     %3 = cmpi slt, %1, %c1_i64 : i64
# CHECK-NEXT:     cond_br %3, ^bb1, ^bb2(%false, %c1_i64, %c1_i64 : i1, i64, i64)
# CHECK-NEXT:   ^bb1:  // pred: ^bb0
# CHECK-NEXT:     %4 = "jlir.undef"() : () -> !jlir.Int64
# CHECK-NEXT:     %5 = "jlir.undef"() : () -> !jlir.Int64
# CHECK-NEXT:     %6 = "jlir.convertstd"(%4) : (!jlir.Int64) -> i64
# CHECK-NEXT:     %7 = "jlir.convertstd"(%5) : (!jlir.Int64) -> i64
# CHECK-NEXT:     br ^bb2(%true, %6, %7 : i1, i64, i64)
# CHECK-NEXT:   ^bb2(%8: i1, %9: i64, %10: i64):  // 2 preds: ^bb0, ^bb1
# CHECK-NEXT:     %11 = xor %8, %true : i1
# CHECK-NEXT:     cond_br %11, ^bb3(%9, %10, %c0_i64 : i64, i64, i64), ^bb7(%c0_i64 : i64)
# CHECK-NEXT:   ^bb3(%12: i64, %13: i64, %14: i64):  // 2 preds: ^bb2, ^bb6
# CHECK-NEXT:     %15 = "jlir.convertstd"(%13) : (i64) -> !jlir.Int64
# CHECK-NEXT:     %16 = addi %14, %12 : i64
# CHECK-NEXT:     %17 = "jlir.==="(%15, %2) : (!jlir.Int64, !jlir.Int64) -> !jlir.Bool
# CHECK-NEXT:     %18 = "jlir.convertstd"(%17) : (!jlir.Bool) -> i1
# CHECK-NEXT:     cond_br %18, ^bb4, ^bb5
# CHECK-NEXT:   ^bb4:  // pred: ^bb3
# CHECK-NEXT:     %19 = "jlir.undef"() : () -> !jlir.Int64
# CHECK-NEXT:     %20 = "jlir.undef"() : () -> !jlir.Int64
# CHECK-NEXT:     %21 = "jlir.convertstd"(%19) : (!jlir.Int64) -> i64
# CHECK-NEXT:     %22 = "jlir.convertstd"(%20) : (!jlir.Int64) -> i64
# CHECK-NEXT:     br ^bb6(%21, %22, %true : i64, i64, i1)
# CHECK-NEXT:   ^bb5:  // pred: ^bb3
# CHECK-NEXT:     %23 = addi %13, %c1_i64 : i64
# CHECK-NEXT:     br ^bb6(%23, %23, %false : i64, i64, i1)
# CHECK-NEXT:   ^bb6(%24: i64, %25: i64, %26: i1):  // 2 preds: ^bb4, ^bb5
# CHECK-NEXT:     %27 = xor %26, %true : i1
# CHECK-NEXT:     cond_br %27, ^bb3(%24, %25, %16 : i64, i64, i64), ^bb7(%16 : i64)
# CHECK-NEXT:   ^bb7(%28: i64):  // 2 preds: ^bb2, ^bb6
# CHECK-NEXT:     return %28 : i64
# CHECK-NEXT:   }
# CHECK-NEXT: }

# CHECK:   llvm.func @"Tuple{typeof(Main.gauss), Int64}"(%arg0: !llvm.ptr<struct<"struct_jl_value_type", opaque>>, %arg1: i64) -> i64 attributes {llvm.emit_c_interface, sym_visibility = "nested"} {
# CHECK-NEXT:     %0 = llvm.mlir.constant({{[0-9]+}} : i64) : i64
# CHECK-NEXT:     %1 = llvm.mlir.constant({{[0-9]+}} : i64) : i64
# CHECK-NEXT:     %2 = llvm.mlir.constant(false) : i1
# CHECK-NEXT:     %3 = llvm.mlir.constant(true) : i1
# CHECK-NEXT:     %4 = llvm.icmp "sle" %0, %arg1 : i64
# CHECK-NEXT:     %5 = llvm.select %4, %arg1, %1 : i1, i64
# CHECK-NEXT:     %6 = llvm.icmp "slt" %5, %0 : i64
# CHECK-NEXT:     llvm.cond_br %6, ^bb1, ^bb2(%2, %0, %0 : i1, i64, i64)
# CHECK-NEXT:   ^bb1:  // pred: ^bb0
# CHECK-NEXT:     %7 = llvm.mlir.undef : i64
# CHECK-NEXT:     llvm.br ^bb2(%3, %7, %7 : i1, i64, i64)
# CHECK-NEXT:   ^bb2(%8: i1, %9: i64, %10: i64):  // 2 preds: ^bb0, ^bb1
# CHECK-NEXT:     %11 = llvm.xor %8, %3  : i1
# CHECK-NEXT:     llvm.cond_br %11, ^bb3(%9, %10, %1 : i64, i64, i64), ^bb7(%1 : i64)
# CHECK-NEXT:   ^bb3(%12: i64, %13: i64, %14: i64):  // 2 preds: ^bb2, ^bb6
# CHECK-NEXT:     %15 = llvm.add %14, %12  : i64
# CHECK-NEXT:     %16 = llvm.icmp "eq" %13, %5 : i64
# CHECK-NEXT:     llvm.cond_br %16, ^bb4, ^bb5
# CHECK-NEXT:   ^bb4:  // pred: ^bb3
# CHECK-NEXT:     %17 = llvm.mlir.undef : i64
# CHECK-NEXT:     llvm.br ^bb6(%17, %17, %3 : i64, i64, i1)
# CHECK-NEXT:   ^bb5:  // pred: ^bb3
# CHECK-NEXT:     %18 = llvm.add %13, %0  : i64
# CHECK-NEXT:     llvm.br ^bb6(%18, %18, %2 : i64, i64, i1)
# CHECK-NEXT:   ^bb6(%19: i64, %20: i64, %21: i1):  // 2 preds: ^bb4, ^bb5
# CHECK-NEXT:     %22 = llvm.xor %21, %3  : i1
# CHECK-NEXT:     llvm.cond_br %22, ^bb3(%19, %20, %15 : i64, i64, i64), ^bb7(%15 : i64)
# CHECK-NEXT:   ^bb7(%23: i64):  // 2 preds: ^bb2, ^bb6
# CHECK-NEXT:     llvm.return %23 : i64
# CHECK-NEXT:   }
# CHECK-NEXT:   llvm.func @"_mlir_ciface_Tuple{typeof(Main.gauss), Int64}"(%arg0: !llvm.ptr<struct<"struct_jl_value_type", opaque>>, %arg1: i64) -> i64 attributes {llvm.emit_c_interface, sym_visibility = "nested"} {
# CHECK-NEXT:     %0 = llvm.call @"Tuple{typeof(Main.gauss), Int64}"(%arg0, %arg1) : (!llvm.ptr<struct<"struct_jl_value_type", opaque>>, i64) -> i64
# CHECK-NEXT:     llvm.return %0 : i64
# CHECK-NEXT:   }
# CHECK-NEXT: }

max_stars_repo_name: JuliaLabs/brutus

max_issues_repo_name: JuliaLabs/brutus

max_forks_repo_name: JuliaLabs/brutus


text:
using RvSpectML
using CSV

function calc_rvs_from_ccf_total(ccfs::AbstractArray{T1,2}, pipeline::PipelinePlan; v_grid::AbstractVector{T2}, times::AbstractVector{T3},
                                 alg_fit_rv::AbstractMeasureRvFromCCF = RVFromCCF.MeasureRvFromCCFGaussian(),
                                 bin_nightly::Bool = true, bin_consecutive::Integer = 0 , recalc::Bool = false, verbose::Bool = true) where {T1<:Real, T2<:Real, T3<:Real }
   @assert length(v_grid) == size(ccfs,1)
   need_to!(pipeline_plan, :rvs_ccf_total)
   if need_to(pipeline_plan, :rvs_ccf_total)
      if verbose println("# Measuring RVs from CCF.")  end
      @assert !need_to(pipeline_plan,:ccf_total)
      #fit_gaussian_to_ccf = RVFromCCF.MeasureRvFromCCFGaussian()
      #fit_quadratic_to_ccf = RVFromCCF.MeasureRvFromCCFQuadratic()
      #rvs_ccf = RVFromCCF.measure_rvs_from_ccf(v_grid,ccfs,alg=fit_gaussian_to_ccf)
      #rvs_ccf = RVFromCCF.measure_rvs_from_ccf(v_grid,ccfs,alg=fit_quadratic_to_ccf)
      rvs_ccf = measure_rvs_from_ccf(v_grid,ccfs,alg=alg_fit_rv)

      if bin_nightly
         rms_rv_nightly = bin_rvs_nightly(times=times,rvs=rvs_ccf.-mean(rvs_ccf))
         rms_rv_within_night = rms_rvs_within_night(times=times,rvs=rvs_ccf.-mean(rvs_ccf))
         if verbose   println("# RMS of RVs: ", std(rvs_ccf), "  nightly RVs: ", std(rms_rv_nightly), "   within night ", (rms_rv_within_night) )   end
      elseif 1 <= bin_consecutive <= floor(size(ccfs,2)//2)
         rms_rv_binned= bin_rvs_consecutive(rvs_ccf.-mean(rvs_ccf),bin_consecutive)
         if verbose   println("# RMS of RVs: ", std(rvs_ccf), "  binned RVs: ", std(rms_rv_binned) )   end
      else
         if verbose   println("# RMS of RVs: ", std(rvs_ccf))   end
      end

      if save_data(pipeline_plan, :rvs_ccf_total)
         CSV.write(joinpath(output_dir,target_subdir * "_rvs_ccf.csv"),DataFrame("Time [MJD]"=>times,"CCF RV [m/s]"=>rvs_ccf))
      end
      set_cache!(pipeline, :rvs_ccf_total, rvs_ccf )
      dont_need_to!(pipeline_plan, :rvs_ccf_total)
   end
   if has_cache(pipeline,:rvs_ccf_total) return read_cache(pipeline,:rvs_ccf_total)
   else   @error("Invalid pipeline state.")          end
end


text:
function show_histogram(
    bins, indices;
    truevalue = nothing, xlabel = L"N [-]", ylabel = L"p(N)"
)
    nbins = length(bins)
    counts = zeros(nbins)
    for i in 1:length(indices)
        counts[indices[i]] += 1
    end
    p = plot(
            bins,
            counts./sum(counts),
            xlabel=xlabel,
            ylabel=ylabel,
            legend=false
           )
    if !isnothing(truevalue)
        plot!([truevalue], seriestype="vline",legend=false)
    end
    return p
end

function show_EPSP_trace(times, epsps)
    return plot(
        times,
        epsps,
        xlabel="Time [s]",
        ylabel="EPSC [A]",
        legend=false
    )
end

function flatten(A::AbstractMatrix)
    n,m = size(A)
    return reshape(A, n*m)
end

function posterior_plot(
    fstate, times, epsps;
    truemodel = nothing,
    truestate = nothing,
    showstates = false
)
    Nrng = Array(fstate.model.Nrng)
    prng = Array(fstate.model.prng)
    qrng = Array(fstate.model.qrng)
    σrng = Array(fstate.model.σrng)
    τrng = Array(fstate.model.τrng)

    Nind = Array(fstate.model.Nind)
    pind = Array(fstate.model.pind)
    qind = Array(fstate.model.qind)
    σind = Array(fstate.model.σind)
    τind = Array(fstate.model.τind)

    pE = show_EPSP_trace(times, epsps)
    pN = show_histogram(Nrng, Nind,
            xlabel = L"N [-]", ylabel = L"p(N)")
    !isnothing(truemodel) && plot!([truemodel.N[1]], seriestype="vline",legend=false)
    pp = show_histogram(prng, pind,
            xlabel = L"p [-]", ylabel = L"p(p)")
    !isnothing(truemodel) && plot!([truemodel.p[1]], seriestype="vline",legend=false)
    pq = show_histogram(qrng, qind,
            xlabel = L"q [A]", ylabel = L"p(q)")
    !isnothing(truemodel)g && plot!([truemodel.q[1]], seriestype="vline",legend=false)
    pσ = show_histogram(σrng, σind,
            xlabel = L"\sigma [A]", ylabel = L"p(\sigma)")
    !isnothing(truemodel) && plot!([truemodel.σ[1]], seriestype="vline",legend=false)
    pτ = show_histogram(τrng, τind,
            xlabel = L"\tau [s]", ylabel = L"p(\tau)")
    !isnothing(truemodel) && plot!([truemodel.τ[1]], seriestype="vline",legend=false)
    if showstates
        pn = histogram(flatten(Array(fstate.state.n)), bins=1:20,
                normalize = :probability,
                xlabel = L"n", ylabel = L"p(n)")
        !isnothing(truestate) && plot!([truestate.n[1,1]], seriestype="vline",legend=false)
        pk = histogram(flatten(Array(fstate.state.k)), bins=1:20,
                normalize = :probability,
                xlabel = L"k", ylabel = L"p(k)")
        !isnothing(truestate) && plot!([truestate.k[1,1]], seriestype="vline",legend=false)
        display(plot(pE, pN, pp, pq, pσ, pτ, pn, pk, layout = (4, 2)))
        return
    end
    display(plot(pE, pN, pp, pq, pσ, pτ, layout = (3, 2)))
end


text:
module CSP

import Base: Condition, open, close, send, recv, map,
             eltype, join_eltype, reduce, filter, merge, select

export Channel, Transport, ChannelClosedError, open, close, send, recv

type ChannelClosedError <: Exception
end

type Channel{T}
    isopen::Bool
    recv_blocker::Condition
    send_blocker::Condition
    bufsize::Int
    buffer::Vector{T}
end

Channel(T::Type=Any, bufsize=0) =
    Channel{T}(true, Condition(), Condition(), bufsize, T[])

open(c::Channel) = c.isopen = true
close(c::Channel) = c.isopen = false
isopen(c::Channel) = c.isopen
bufsize(c::Channel) = c.bufsize
eltype{T}(c::Channel{T}) = T

function send(c::Channel, val)
    if !isopen(c)
        throw(ChannelClosedError())
    end

    if length(c.buffer) < bufsize(c)
        push!(c.buffer, val)
        notify(c.recv_blocker, all=false)
    else
        if !isempty(c.recv_blocker.waitq)
            # if procs are already waiting to recv
            push!(c.buffer, val)
	    notify(c.recv_blocker, all=false)
        else
            # otherwise, block until someone wants to recv
            push!(c.buffer, val)
            wait(c.send_blocker)
	    notify(c.recv_blocker, all=false)
        end
    end
end

function recv(c::Channel)
    # Read out if there is something in the buffer
    if length(c.buffer) > 0
        try
            v = shift!(c.buffer)
            # tell others that you have freed a buffer slot
            notify(c.send_blocker, all=false)
            return v
        catch
            # Continue and wait if cannot shift!
        end
    end
    # Or block till there is something
    if !isopen(c)
        warn("waiting to receive from a closed channel!")
    end

    wait(c.recv_blocker)
    v = shift!(c.buffer)
    notify(c.send_blocker, all=false)
    return v
end

typealias Selected{T} (Channel{T}, T)

function select(chans; bufsize=0)
    out = Channel(Selected, bufsize)
    for c in chans
	@async send(out, recv(c))
    end
    out
end

function map(f::Base.Callable, inp::Channel; typ=eltype(inp), bufsize=bufsize(inp))
    send(out, v)
    @async while true
        out = Channel(typ, bufsize=bufsize)
        send(out, f(recv(inp)))
    end
    out
end

function reduce{T}(f::Base.Callable, v0::T, inp::Channel; bufsize=bufsize(inp))
    out = Channel(T, bufsize)
    store = v0
    @async while true
        send(out, f(store, recv(inp)))
    end
    out
end

function filter{T}(pred::Function, inp::Channel{T}; bufsize=bufsize(inp))
    out = Channel(T, bufsize)
    @async while true
        v = recv(inp)
        if pred(v)
            send(out, v)
	end
    end
    out
end

merge(
        chans...;
        typ=reduce(join_eltype, map(eltype, chans)),
        bufsize=max(map(bufsize, chans)...)
    ) =
    map(x -> x[2], select(chans), typ=typ, bufsize=bufsize)

abstract Transport

function pipe(t::Transport, c::Channel=Channel())
    @async while true
        send(c, recv(t)) # Assumes recv blocks
    end
end

function pipe(c1::Channel, c2::Channel)
    @async while true
        send(c1, recv(c2)) # Assumes recv blocks
    end
end

function pipe(c::Channel, t::Transport)
    @async while true
        send(t, recv(c))
    end
end

(>>>)(a::Union(Channel, Transport), b::Union(Channel, Transport)) = pipe(a, b)

end # module


text:
export Coad1D, set!, step!

const r₁ = (0.1)*1e-6 # minimum droplet size bin, μm->m
# NOTE: We use subscript "1" because Julia is a 1-indexed language. By convention
# .     we reserve the nought subscript for general parameters.
const rm = (100_000.0)*1e-6 # maximum droplet size bin, μm->m
const gmin = 1e-60 # lower bound on permissible bin mass density

"""
    struct Coad1D{
      FT <: AbstractFloat,
      I <: Int,
      FTA1D <: AbstractArray{FT, 1},
      FTA2D <: AbstractArray{FT, 2},
      IA2D <: AbstractArray{I, 2}
    } <: AbstractCoadModel

A collison/coalescence model type.

Typing inspired by: https://github.com/CliMA/RRTMGP.jl/blob/master/src/optics/RTE.jl
"""
struct Coad1D{
  FT <: AbstractFloat,
  I <: Int,
  FTA1D <: AbstractArray{FT, 1},
  FTA2D <: AbstractArray{FT, 2},
  IA2D <: AbstractArray{I, 2}
} <: AbstractCoadModel
    # Grid spacing
    α::FT
    # Collison kernel
    kernel::Symbol

    # Grid definition
    xᵢ::FTA1D
    rᵢ::FTA1D
    Δlnr::FT

    # Initial conditions - don't set by default
    gᵢ::FTA1D

    # Cached computatonal values
    c::FTA2D
    ima::IA2D # courant limits and collison idx LUT
    ck::FTA2D # smoothed collision kernel LUT

    # Do we need to include anything else here?
    function Coad1D(;n::Int64, kernel::Symbol)
        α = 2^(1/n)

        m = ceil(Integer, 1 + 3*log(rm / r₁)/log(α)) # number of mass bins to populate
        Δlnr = log(α) / 3  # constant grid distance of logarithmic grid

        println("Instatiating model with $m mass bins...")
        rᵢ = r₁*(α.^((collect(1:m) .- 1)./3)) # meter
        xᵢ = mass_from_r.(rᵢ) # kg
        gᵢ = similar(xᵢ) .* 0

        println("  Pre-caching grid Courant limits...")
        c, ima = courant(xᵢ)
        println("  Pre-caching $kernel kernel lookup table...")
        ck = kernels(xᵢ, kernel)

        FT = eltype(xᵢ)
        I = eltype(ima)
        FTA1D = typeof(xᵢ)
        FTA2D = typeof(c)
        IA2D = typeof(ima)

        return new{
          FT, I, FTA1D, FTA2D, IA2D
        }(
          α, kernel, xᵢ, rᵢ, Δlnr, gᵢ, c, ima, ck
        )
    end
end

function Base.show(io::IO, model::Coad1D)
  m = length(model.xᵢ)
  print(io, "Coad1D Model ($m radii ∈ ($(r₁*1e6), $(rm*1e6)) μm, kernel=$(model.kernel))")
end

"""
    courant(x::AbstractArray{FT, 1}) where {FT <: AbstractFloat}

Given a binned mass grid x, pre-compute the Courant limits of the B98 advection
flux algorithm associated with each combination of potential collisions sampled
from bins on this grid.
"""
@inline function courant(x::AbstractArray{FT, 1}) where {FT <: AbstractFloat}

  # Set up return matrices; first inspect length of our mass grid, 'x'
  m = length(x)
  # FT = eltype(x)
  c = zeros(FT, m, m) # Courant numbers for limiting advection flux
  ima = zeros(Int64, m, m) # This is basically keeping track of the left bound of the

  # Compute Courant limits
  for i ∈ 1:m
    for j ∈ i:m

      x0 = x[i] + x[j]  # Summed / total mass from collision
      for k ∈ j:m
        if (x[k] ≥ x0) && (x[k-1] < x0) # There is probably an easier way to exploit the size of the collision here than linear searching for bounding masses
          if (c[i, j] < (1 - 1e-8))
            kk = k - 1
            c[i, j] = log(x0 / x[k-1])
          else
            c[i, j] = 0.0
            kk = k
          end
          ima[i, j] = min(m - 1, kk)
          break
        end
      end # k loop

      # Copy over diagonal for symmetry
      c[j, i] = c[i, j]
      ima[j, i] = ima[i, j]

    end # j loop
  end # i loop

  return c, ima
end

"""
    set!(model::Coad1D, dist::AbstractSizeDist{FT}) where {FT <: AbstractFloat}

Set the state vector (mass concentraton, model.gᵢ) according to some specified
size distribution (number concentration).
"""
function set!(model::Coad1D, dist::AbstractSizeDist{FT}) where {FT <: AbstractFloat}
  # Initial droplet distribution; g(y, t) = 3x² * n(x, t), eq. 2 from B98
  model.gᵢ .= (3 .* model.xᵢ.^2) .* nc.(Ref(dist), model.xᵢ) # kg / m3
end

"""
    step!(model::Coad1D, Δt::FT) where {FT <: AbstractFloat}

Integrate one step of the SCE forward in time by increment Δt. The model's state
will be updated to reflect the change in mass distribution gᵢ. = g(xᵢ, t+Δt).
"""
@inline function step!(model::Coad1D, Δt::FT) where {FT <: AbstractFloat}

  # Lower and Upper integration limit i0, i1
  i0, i1 = find_bounds(model.gᵢ, gmin=gmin)

  # Main collision/coalescence loop
  @inbounds for i ∈ i0:i1
    @inbounds for j ∈ i:i1
      k = model.ima[i, j]  # Get pre-computed index of target coalescence bin
      kp = k + 1

      # NOTE - Here we handle a weird edge condition in the initialization of
      # ima, where the computed target bin index could actually be invalid (0)
      # in the original code. This only seems to happen for the very largest
      # bins, which are going to overshoot the grid anyway (leak/diffuse mass)
      # but we need this sanity check to avoid indexing outside of our array.
      if k < 1
        continue
      end

      # We did not pre-scale the collison kernels LUT by Δt and Δy in this
      # version so we have to account for them here.
      x0 = model.ck[i, j] * model.gᵢ[i] * model.gᵢ[j] * Δt * model.Δlnr
      x0 = min(x0, model.gᵢ[i] * model.xᵢ[j])
      if j != k # Not sure what's going on here.
        x0 = min(x0, model.gᵢ[j] * model.xᵢ[i])
      end

      gsi = x0 / model.xᵢ[j]
      gsj = x0 / model.xᵢ[i]
      gsk = gsi + gsj
      model.gᵢ[i] = model.gᵢ[i] - gsi
      model.gᵢ[j] = model.gᵢ[j] - gsj
      gk = model.gᵢ[k] + gsk

      if gk > gmin

        # ORIGINAL
        # x1 = log(gᵢ[kp] / gk + 1e-60)

        # MODIFIED - Apply a limiter to avoid negative args to log.
        # We note that this may not strictly obey the formulation of the flux
        # algorithm, but this tends to work okay in practice.
        log_arg = model.gᵢ[kp] / gk + 1e-60
        log_arg = max(1e-60, log_arg)
        x1 = log(log_arg)

        # We did not pre-scale the cached Courant limits for this version so we
        # have to apply them when estimating the flux directly.
        flux = gsk / x1 * (exp(0.5 * x1) - exp(x1 * (0.5 - (model.c[i, j]/3/model.Δlnr))))
        flux = min(flux, gk)
        model.gᵢ[k] = gk - flux
        model.gᵢ[kp] = model.gᵢ[kp] + flux
      end

    end # j
  end # i

end

"""
    find_bounds(g::AbstractArray{FT, 1}; gmin::FT=1e-60) where {FT <: AbstractFloat}

Given a binned discretization of a droplet mass distribution, find the range of
bin indices which have non-negligble mass. This function is intended to help
accelerate other computations which require analyzing a process on the full mass
grid by minimizing the number of bins which need to be inspected.
"""
@inline function find_bounds(g::AbstractArray{FT, 1}; gmin::FT=1e-60) where {FT <: AbstractFloat}

  m = length(g)

  # TODO: refactor since this can be wrapped in a single array function
  # This basically sets a "focus" in the array where we have mass that needs
  # to get collided / advected around, so we don't waste cycles on empty bins.
  # In practice seems to be a limiter on numerical issues.
  i0 = 1
  @inbounds for i ∈ 1:m-1
    i0 = i
    if g[i] > gmin
      break
    end
  end
  i1 = m-1
  @inbounds for i ∈ m-1:-1:1
    i1 = i
    if g[i] > gmin
      break
    end
  end

  return i0, i1
end


text:
function ChainRulesCore.rrule(d::Flux.Dense, x)
    θ = d.W*x + d.b

    Y, ∂Ȳ = frule(broadcast, d.σ, θ)

    function dense_pullback(ȳ)
        ΔȲ = ∂Ȳ(nothing, nothing, ȳ)

        ∂W = InplaceableThunk(
            @thunk(ΔȲ*transpose(x)),
            x̄->ΔȲ*transpose(x)+x̄
        )
        ∂b = InplaceableThunk(
            @thunk(ΔȲ),
            x̄->ΔȲ+x̄
        )
        ∂σ = Zero()

        ∂dense = (W = ∂W, b = ∂b, σ = ∂σ)

        return (∂dense, @thunk(transpose(d.W)*ΔȲ))
    end
    return Y, dense_pullback
end


text:
export Lens, set, get, modify
export @lens
export set, get, modify
using ConstructionBase
export setproperties
export constructorof

import Base: get
using Base: getproperty

"""
    Lens

A `Lens` allows to access or replace deeply nested parts of complicated objects.

# Example
```jldoctest
julia> using Setfield

julia> struct T;a;b; end

julia> obj = T("AA", "BB")
T("AA", "BB")

julia> lens = @lens _.a
(@lens _.a)

julia> get(obj, lens)
"AA"

julia> set(obj, lens, 2)
T(2, "BB")

julia> obj
T("AA", "BB")

julia> modify(lowercase, obj, lens)
T("aa", "BB")
```

# Interface
Concrete subtypes of `Lens` have to implement
* `set(obj, lens, val)`
* `get(obj, lens)`

These must be pure functions, that satisfy the three lens laws:

```jldoctest; output = false, setup = :(using Setfield; (≅ = (==)); obj = (a="A", b="B"); lens = @lens _.a; val = 2; val1 = 10; val2 = 20)
@assert get(set(obj, lens, val), lens) ≅ val
        # You get what you set.
@assert set(obj, lens, get(obj, lens)) ≅ obj
        # Setting what was already there changes nothing.
@assert set(set(obj, lens, val1), lens, val2) ≅ set(obj, lens, val2)
        # The last set wins.

# output

```
Here `≅` is an appropriate notion of equality or an approximation of it. In most contexts
this is simply `==`. But in some contexts it might be `===`, `≈`, `isequal` or something
else instead. For instance `==` does not work in `Float64` context, because
`get(set(obj, lens, NaN), lens) == NaN` can never hold. Instead `isequal` or
`≅(x::Float64, y::Float64) = isequal(x,y) | x ≈ y` are possible alternatives.

See also [`@lens`](@ref), [`set`](@ref), [`get`](@ref), [`modify`](@ref).
"""
abstract type Lens end

"""
    modify(f, obj, l::Lens)

Replace a deeply nested part `x` of `obj` by `f(x)`. See also [`Lens`](@ref).
"""
function modify end

"""
    get(obj, l::Lens)

Access a deeply nested part of `obj`. See also [`Lens`](@ref).
"""
function get end

"""
    set(obj, l::Lens, val)

Replace a deeply nested part of `obj` by `val`. See also [`Lens`](@ref).
"""
function set end

@inline function modify(f, obj, l::Lens)
    old_val = get(obj, l)
    new_val = f(old_val)
    set(obj, l, new_val)
end

struct IdentityLens <: Lens end
get(obj, ::IdentityLens) = obj
set(obj, ::IdentityLens, val) = val

struct PropertyLens{fieldname} <: Lens end

function get(obj, l::PropertyLens{field}) where {field}
    getproperty(obj, field)
end

@inline function set(obj, l::PropertyLens{field}, val) where {field}
    patch = (;field => val)
    setproperties(obj, patch)
end

struct ComposedLens{LO, LI} <: Lens
    outer::LO
    inner::LI
end

"""
    compose([lens₁, [lens₂, [lens₃, ...]]])

Compose `lens₁`, `lens₂` etc. There is one subtle point here:
While the two composition orders `(lens₁ ∘ lens₂) ∘ lens₃` and `lens₁ ∘ (lens₂ ∘ lens₃)` have equivalent semantics,
their performance may not be the same. The compiler tends to optimize right associative composition
(second case) better then left associative composition.

The compose function tries to use a composition order, that the compiler likes. The composition order is therefore not part of the stable API.
"""
function compose end
compose() = IdentityLens()
compose(l::Lens) = l
compose(::IdentityLens, ::IdentityLens) = IdentityLens()
compose(::IdentityLens, l::Lens) = l
compose(l::Lens, ::IdentityLens) = l
compose(outer::Lens, inner::Lens) = ComposedLens(outer, inner)
function compose(l1::Lens, ls::Lens...)
    compose(l1, compose(ls...))
end

"""
    lens₁ ∘ lens₂

Compose lenses `lens₁`, `lens₂`, ..., `lensₙ` to access nested objects.

# Example
```jldoctest
julia> using Setfield

julia> obj = (a = (b = (c = 1,),),);

julia> la = @lens _.a
       lb = @lens _.b
       lc = @lens _.c
       lens = la ∘ lb ∘ lc
(@lens _.a.b.c)

julia> get(obj, lens)
1
```
"""
Base.:∘(l1::Lens, l2::Lens) = compose(l1, l2)

function get(obj, l::ComposedLens)
    inner_obj = get(obj, l.outer)
    get(inner_obj, l.inner)
end

function set(obj,l::ComposedLens, val)
    inner_obj = get(obj, l.outer)
    inner_val = set(inner_obj, l.inner, val)
    set(obj, l.outer, inner_val)
end

struct IndexLens{I <: Tuple} <: Lens
    indices::I
end

Base.@propagate_inbounds function get(obj, l::IndexLens)
    getindex(obj, l.indices...)
end
Base.@propagate_inbounds function set(obj, l::IndexLens, val)
    setindex(obj, val, l.indices...)
end

struct DynamicIndexLens{F} <: Lens
    f::F
end

Base.@propagate_inbounds get(obj, I::DynamicIndexLens) = obj[I.f(obj)...]

Base.@propagate_inbounds set(obj, I::DynamicIndexLens, val) =
    setindex(obj, val, I.f(obj)...)

"""
    FunctionLens(f)
    @lens f(_)

Lens with [`get`](@ref) method definition that simply calls `f`.
[`set`](@ref) method for each function `f` must be implemented manually.
Use `methods(set, (Any, Setfield.FunctionLens, Any))` to get a list of
supported functions.

Note that `FunctionLens` flips the order of composition; i.e.,
`(@lens f(_)) ∘ (@lens g(_)) == @lens g(f(_))`.

# Example
```jldoctest
julia> using Setfield

julia> obj = ((1, 2), (3, 4));

julia> lens = (@lens first(_)) ∘ (@lens last(_))
(@lens last(first(_)))

julia> get(obj, lens)
2

julia> set(obj, lens, '2')
((1, '2'), (3, 4))
```

# Implementation

To use `myfunction` as a lens, define a `set` method with the following
signature:

```julia
Setfield.set(obj, ::typeof(@lens myfunction(_)), val) = ...
```

`typeof` is used above instead of `FunctionLens` because how actual
type of `@lens myfunction(_)` is implemented is not the part of stable
API.
"""
struct FunctionLens{f} <: Lens end
FunctionLens(f) = FunctionLens{f}()

get(obj, ::FunctionLens{f}) where f = f(obj)


text:
module QuHamiltonian

include("pauli_matrix.jl")

include("sites.jl")
include("region.jl")
include("local_operators.jl")
include("placeholder.jl")

include("sub_hamiltonian.jl")
include("symbolic.jl")
include("eval.jl")

include("conversions.jl")

# include("precompile.jl")
# _precompile_()

end # module

max_stars_repo_name: Roger-luo/QuHamiltonian

max_issues_repo_name: Roger-luo/QuHamiltonian

max_forks_repo_name: Roger-luo/QuHamiltonian


text:
module DeckEncoder # LoRDeckCodes

using ..LoRDeckCodes: Base32, VarintTranslator, Deck, CardCodeAndCount, ArgumentException

const DECK_CODE_FORMAT = 1
const MAX_KNOWN_VERSION = 1
const CARD_CODE_LENGTH = 7

const FactionCodeToIntIdentifier = Dict{String,Int}(
    "DE" => 0,
    "FR" => 1,
    "IO" => 2,
    "NX" => 3,
    "PZ" => 4,
    "SI" => 5,
    "BW" => 6,
    "SH" => 7,
    "MT" => 9
)

const IntIdentifierToFactionCode = Dict{Int,String}(
    0 => "DE",
    1 => "FR",
    2 => "IO",
    3 => "NX",
    4 => "PZ",
    5 => "SI",
    6 => "BW",
    7 => "SH",
    9 => "MT"
)

"""
    encode_deck(cards::Vector{CardCodeAndCount})::String
"""
function encode_deck(cards::Vector{CardCodeAndCount})::String
    encode_deck(cards, MAX_KNOWN_VERSION)
end

"""
    encode_deck(deck::Deck)::String
"""
function encode_deck(deck::Deck)::String
    encode_deck(deck.cards, deck.version)
end

function encode_deck(cards::Vector{CardCodeAndCount}, version::UInt8)::String
    stream = IOBuffer()
    !isvalid(cards) && throw(ArgumentException("The provided deck contains invalid card codes."))
    format_and_version = UInt8((DECK_CODE_FORMAT << 4) | version)
    write(stream, format_and_version)
    of3 = Vector{CardCodeAndCount}()
    of2 = Vector{CardCodeAndCount}()
    of1 = Vector{CardCodeAndCount}()
    ofN = Vector{CardCodeAndCount}()
    for card in cards
        if card.count == 3
            push!(of3, card)
        elseif card.count == 2
            push!(of2, card)
        elseif card.count == 1
            push!(of1, card)
        elseif card.count < 1
            throw(ArgumentException(string("Invalid count of ", card.count, " for card ", card.code)))
        else
            push!(ofN, card)
        end
    end
    groupedOf3s = getGroupedOfs(of3)
    groupedOf2s = getGroupedOfs(of2)
    groupedOf1s = getGroupedOfs(of1)
    encodeGroupOf(stream, sortGroupOf(groupedOf3s))
    encodeGroupOf(stream, sortGroupOf(groupedOf2s))
    encodeGroupOf(stream, sortGroupOf(groupedOf1s))
    encodeNOfs(stream, sort(ofN, by = c -> c.code))
    seekstart(stream)
    Base32.encode(read(stream))
end

function getGroupedOfs(list::Vector{CardCodeAndCount})::Vector{Vector{CardCodeAndCount}}
    result = Vector{Vector{CardCodeAndCount}}()
    while !isempty(list)
        currentSet = Vector{CardCodeAndCount}()
        firstCardCodeAndCount = first(list)
        push!(currentSet, firstCardCodeAndCount)
        deleteat!(list, 1)
        for i in length(list):-1:1
            currentCardCodeAndCount = list[i]
            if currentCardCodeAndCount.set == firstCardCodeAndCount.set && currentCardCodeAndCount.faction == firstCardCodeAndCount.faction
                push!(currentSet, currentCardCodeAndCount)
                deleteat!(list, i)
            end
        end
        push!(result, currentSet)
    end
    result
end

function sortGroupOf(group::Vector{Vector{CardCodeAndCount}})::Vector{Vector{CardCodeAndCount}}
    groupOf = sort(group, by = g -> length(g))
    for i = 1:length(groupOf)
        groupOf[i] = sort(groupOf[i], by = c -> c.code)
    end
    groupOf
end

function encodeGroupOf(stream::IOBuffer, groupOf::Vector{Vector{CardCodeAndCount}})
    write(stream, VarintTranslator.get_varint(length(groupOf)))
    for currentList in groupOf
        write(stream, VarintTranslator.get_varint(length(currentList)))
        currentCardCodeAndCount = first(currentList)
        write(stream, VarintTranslator.get_varint(currentCardCodeAndCount.set))
        write(stream, VarintTranslator.get_varint(FactionCodeToIntIdentifier[currentCardCodeAndCount.faction]))
        for cd in currentList
            write(stream, VarintTranslator.get_varint(cd.number))
        end
    end
end

function encodeNOfs(stream::IOBuffer, nOfs::Vector{CardCodeAndCount})
    for card in nOfs
        write(stream, VarintTranslator.get_varint(card.count))
        write(stream, VarintTranslator.get_varint(card.set))
        write(stream, VarintTranslator.get_varint(FactionCodeToIntIdentifier[card.faction]))
        write(stream, VarintTranslator.get_varint(card.number))
    end
end

function isvalid(card::CardCodeAndCount)::Bool
    length(card.code) != CARD_CODE_LENGTH && return false
    try
        parse(Int, card.code[1:2]) # set
        parse(Int, card.code[5:7]) # number
    catch err
        return false
    end
    !haskey(FactionCodeToIntIdentifier, card.faction) && return false
    card.count < 1 && return false
    true
end

function isvalid(cards::Vector{CardCodeAndCount})::Bool
    for card in cards
        !isvalid(card) && return false
    end
    true
end

function isvalid(; deckcode::Union{Nothing,AbstractString} = nothing, cardcode::Union{Nothing,AbstractString} = nothing)::Bool
    if deckcode isa AbstractString
        try
            decode_deck(deckcode)
        catch
            return false
        end
        return true
    elseif cardcode isa AbstractString
        return isvalid(CardCodeAndCount(cardcode, 1))
    else
        return false
    end
end

function decode_deck(deckcode::AbstractString)::Deck
    InvalidDeckCode = ArgumentException("Invalid deck code")
    isempty(deckcode) && throw(InvalidDeckCode)
    bytes = UInt8[]
    try
        bytes = Base32.decode(String(deckcode))
    catch
        throw(InvalidDeckCode)
    end
    stream = IOBuffer(bytes)
    (firstbyte,) = read(stream, 1)
    format = firstbyte >> 4
    version = firstbyte & 0x0f
    version > MAX_KNOWN_VERSION && throw(ArgumentException("'The provided code requires a higher version of this library; please update."))
    cards = Vector{CardCodeAndCount}()
    for i in 3:-1:1
        numGroupOfs = VarintTranslator.pop_varint(stream)
        for j in 1:numGroupOfs
            numOfsInThisGroup = VarintTranslator.pop_varint(stream)
            set = VarintTranslator.pop_varint(stream)
            faction_id = VarintTranslator.pop_varint(stream)
            for k in 1:numOfsInThisGroup
                 number = VarintTranslator.pop_varint(stream)
                 push!(cards, CardCodeAndCount(set, IntIdentifierToFactionCode[faction_id], number, i))
            end
        end
    end
    while !eof(stream)
        count = VarintTranslator.pop_varint(stream)
        set = VarintTranslator.pop_varint(stream)
        faction_id = VarintTranslator.pop_varint(stream)
        number = VarintTranslator.pop_varint(stream)
        push!(cards, CardCodeAndCount(set, IntIdentifierToFactionCode[faction_id], number, count))
    end
    Deck(cards, version)
end

"""
    Deck(deckcode::AbstractString)::Deck
"""
function Deck(deckcode::AbstractString)::Deck
    decode_deck(deckcode)
end

function Deck(cards::Vector{CardCodeAndCount})::Deck
    Deck(cards, MAX_KNOWN_VERSION)
end

end # module LoRDeckCodes.DeckEncoder


text:
# Test for GaussRadau()

n = 1
x, w = gaussradau(n)
@test_approx_eq x[1] -1.
@test_approx_eq w[1] 2.

n = 2
x, w = gaussradau(n)
@test_approx_eq x[1] -1.
@test_approx_eq x[2] 1./3.
@test_approx_eq w[1] .5
@test_approx_eq w[2] 1.5

tol = 1e-14
n = 42
ntests = 4
x, w = gaussradau(n)
@test length(x) == n && length(w) == n
@test abs(x[37] - 0.908847278001044) < tol
@test abs(w[37] - 0.031190846817016) < tol
@test x[1] == -1

@test_approx_eq dot( w,exp(x)) exp(1)-exp(-1)


text:
# Inference algorithm for generalized bilinear model (GBM) with g(E(Y)) = XA' + BZ' + XCZ' + UDV'.
#
# Copyright (c) 2020: Jeffrey W. Miller.
# This file is released under the MIT "Expat" License.

module Infer

export infer

using LinearAlgebra: diag
import LinearAlgebra
Identity = LinearAlgebra.I
Diagonal = LinearAlgebra.Diagonal

include("common.jl"); using .Common
include("outcome_nb.jl"); using .Outcome_NB

nchoose2(n) = div(n*(n-1),2)
block(i,K) = ((i-1)*K+1):(i*K)

# _________________________________________________________________________________________________

# Compute standard errors for estimates of the parameters of a Negative-Binomial generalized bilinear model (NB-GBM).
#
# INPUTS:
#    Y, X, Z, A, B, C, D, U, V, S, T, and omega (see gbm_estimation for parameter descriptions).
#
# OUTPUTS:
#    se_A = J-by-K matrix of (approximate) standard errors for A estimates
#    se_B = I-by-L matrix of (approximate) standard errors for B estimates
#    se_C = K-by-L matrix of (approximate) standard errors for C estimates
#    se_U = I-by-M matrix of (approximate) standard errors for U estimates
#    se_V = J-by-M matrix of (approximate) standard errors for V estimates
#    se_S = I-by-1 matrix of (approximate) standard errors for S estimates
#    se_T = J-by-1 matrix of (approximate) standard errors for T estimates
#
function infer(Y,X,Z,A,B,C,D,U,V,S,T,omega; prior=Prior(), Offset=0.0*Y)
    I,K = size(X)
    J,L = size(Z)
    M = length(D)

    p = prior

    Mu = zeros(I,J)
    W = zeros(I,J)
    E = zeros(I,J)
    logMu = X*A' + B*Z' + X*C*Z' + U*(D.*V') + Offset
    r = exp.(-(S.+T'.+omega))
    compute_MuWE!(Mu,W,E,Y,logMu,r)

    # Compute Fisher information matrices
    F_A = [X'*(W[:,j].*X) for j=1:J]
    F_A_inv = [inv(F_A[j] + p.lambda_A*Identity) for j=1:J]

    F_B = [Z'*(W[i,:].*Z) for i=1:I]
    F_B_inv = [inv(F_B[i] + p.lambda_B*Identity) for i=1:I]

    F_C = hvcat(L, permutedims([X'*((W*(Z[:,l1].*Z[:,l2])).*X) for l1=1:L, l2=1:L])...)
    F_C_inv = inv(F_C + p.lambda_C*Identity)

    UD = U.*D'
    VD = V.*D'
    F_U = [VD'*(W[i,:].*VD) for i=1:I]
    F_V = [UD'*(W[:,j].*UD) for j=1:J]
    # F_D = [sum((U[:,m1].*U[:,m2]).*W.*(V[:,m1].*V[:,m2])') for m1=1:M, m2=1:M]
    F_D = compute_F_D(U,V,W)
    F_U_inv = [inv(F_U[i] + p.lambda_U*Identity) for i=1:I]
    F_V_inv = [inv(F_V[j] + p.lambda_V*Identity) for j=1:J]
    F_D_inv = inv(F_D + p.lambda_D*Identity)

    if true
        # Compute full matrix for (U,V,Q) only, and propagate uncertainty to A,B,C.

        # Uncertainty in U and V
        F_U_by_V = zeros(I*M,J*M)
        for i = 1:I, j = 1:J
            for m1 = 1:M, m2 = 1:M
                F_U_by_V[(i-1)*M+m1, (j-1)*M+m2] = W[i,j]*V[j,m1]*D[m1]*U[i,m2]*D[m2]
            end
        end

        Q_U,Q_V = compute_constraint_Jacobian_UV(X,Z,U,V)

        # F_UiQ = vcat([F_U_inv[i]*view(Q_U,block(i,M),:) for i=1:I]...)
        F_UiQ = zeros(size(Q_U))
        for i = 1:I
            F_UiQ[block(i,M),:] = F_U_inv[i]*view(Q_U,block(i,M),:)
        end
        F_VUUiQ = F_U_by_V'*F_UiQ
        inv_QFUQ = inv(Q_U'*F_UiQ)
        # F_UiUV = vcat([F_U_inv[i]*view(F_U_by_V,block(i,M),:) for i=1:I]...)
        F_UiUV = zeros(I*M,J*M)
        for i = 1:I
            F_UiUV[block(i,M),:] = F_U_inv[i]*view(F_U_by_V,block(i,M),:)
        end
        F_VUUiUV = F_U_by_V'*F_UiUV  # <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<  This is the expensive term?  O(I * J^2 * M^3)
        Temp = -F_VUUiUV + F_VUUiQ * inv_QFUQ * F_VUUiQ'
        for j=1:J; Temp[block(j,M),block(j,M)] += F_V[j] + p.lambda_V*Identity; end

        C_V = inv([Temp  Q_V ; Q_V' zeros(size(Q_V,2),size(Q_V,2))])[1:J*M,1:J*M]

        diag_T = vcat(map(diag,F_U_inv)...) - vec(sum(F_UiQ'.*(inv_QFUQ * F_UiQ'), dims=1))
        F_UVT = F_UiUV' - F_VUUiQ * inv_QFUQ * F_UiQ'
        var_U_approx = diag_T + vec(sum(F_UVT.*(C_V*F_UVT),dims=1))
        var_V_approx = diag(C_V)

        # dW_dMu_Mu = Mu.*(r./(r.+Mu)).^2
        # dE_dMu_Mu = -Mu.*r.*(r.+Y)./(r.+Mu).^2
        # grad_U = E*(V.*D')
        # grad_V = E'*(U.*D')
        # var_A_baseline = vcat(map(diag,F_A_inv)...)
        # var_B_baseline = vcat(map(diag,F_B_inv)...)
        # var_V_from_A = compute_var_V_from_A(U.*D',X,var_A_baseline,F_V_inv,grad_V,dW_dMu_Mu,dE_dMu_Mu)
        # var_V_from_B = compute_var_V_from_B(U.*D',Z,var_B_baseline,F_V_inv,grad_V,dW_dMu_Mu,dE_dMu_Mu)
        # var_U_from_B = compute_var_U_from_B(V.*D',Z,var_B_baseline,F_U_inv,grad_U,dW_dMu_Mu,dE_dMu_Mu)
        # var_V_approx += var_V_from_A + var_V_from_B
        # var_U_approx += var_U_from_B

        # Propagate uncertainty in U and V through to D
        var_D_approx = diag(F_D_inv)  # TODO: NOT YET IMPLEMENTED

        # Propagate uncertainty in U and V through to A and B
        grad_A = E'*X
        grad_B = E*Z
        dW_dMu_Mu = Mu.*(r./(r.+Mu)).^2
        dE_dMu_Mu = -Mu.*r.*(r.+Y)./(r.+Mu).^2
        var_B_from_U = compute_var_B_from_U(Z,V,D,var_U_approx,F_B_inv,grad_B,dW_dMu_Mu,dE_dMu_Mu)
        var_B_from_V = compute_var_B_from_V(Z,U,D,var_V_approx,F_B_inv,grad_B,dW_dMu_Mu,dE_dMu_Mu)
        var_A_from_V = compute_var_A_from_V(X,U,D,var_V_approx,F_A_inv,grad_A,dW_dMu_Mu,dE_dMu_Mu)
        var_A_from_U = compute_var_A_from_U(X,V,D,var_U_approx,F_A_inv,grad_A,dW_dMu_Mu,dE_dMu_Mu)
        var_A_approx = vcat(map(diag,F_A_inv)...) + var_A_from_U + var_A_from_V
        var_B_approx = vcat(map(diag,F_B_inv)...) + var_B_from_U + var_B_from_V

        # Propagate uncertainty in A and B through to C
        grad_C = vec(X'*E*Z)
        var_C_from_A = compute_var_C_from_A(X,Z,F_A_inv,F_C_inv,grad_C,dW_dMu_Mu,dE_dMu_Mu)
        var_C_from_B = compute_var_C_from_B(X,Z,F_B_inv,F_C_inv,grad_C,dW_dMu_Mu,dE_dMu_Mu)
        var_C_approx = diag(F_C_inv) + var_C_from_A + var_C_from_B

    else

        # For testing purposes, compute the full Fisher information matrix
        F_A_full = zeros(J*K,J*K); for j = 1:J; F_A_full[block(j,K),block(j,K)] = F_A[j]; end
        F_B_full = zeros(I*L,I*L); for i = 1:I; F_B_full[block(i,L),block(i,L)] = F_B[i]; end
        F_C_full = F_C

        F_U_full = zeros(I*M,I*M); for i = 1:I; F_U_full[block(i,M),block(i,M)] = F_U[i]; end
        F_V_full = zeros(J*M,J*M); for j = 1:J; F_V_full[block(j,M),block(j,M)] = F_V[j]; end
        F_D_full = F_D

        F_A_by_B = hvcat(I, permutedims([W[i,j]*X[i,:]*Z[j,:]' for j=1:J, i=1:I])...)
        F_A_by_C = hvcat(L, permutedims([Z[j,l]*X'*(W[:,j].*X) for j=1:J, l=1:L])...)
        F_A_by_U = hvcat(I, permutedims([W[i,j]*X[i,:]*(V[j,:].*D)' for j=1:J, i=1:I])...)
        F_A_by_V = hvcat(J, permutedims([X'*(W[:,j1].*U.*D')*(j1==j2) for j1=1:J, j2=1:J])...)
        F_A_by_D = vcat([X'*(W[:,j].*U.*V[j,:]') for j=1:J]...)

        F_B_by_C = hvcat(L, permutedims([(Z'*(W[i,:].*Z[:,l]))*X[i,:]' for i=1:I, l=1:L])...)
        F_B_by_U = hvcat(I, permutedims([Z'*(W[i1,:].*V.*D')*(i1==i2)  for i1=1:I, i2=1:I])...)
        F_B_by_V = hvcat(J, permutedims([W[i,j]*Z[j,:]*(U[i,:].*D)'   for i=1:I, j=1:J])...)
        F_B_by_D = vcat([Z'*(W[i,:].*V.*U[i,:]') for i=1:I]...)

        F_C_by_U = hvcat(I, permutedims([X[i,:]*((W[i,:].*Z[:,l])'*(V.*D'))  for l=1:L, i=1:I])...)
        F_C_by_V = hvcat(J, permutedims([X'*(W[:,j].*Z[j,l].*U.*D') for l=1:L, j=1:J])...)
        F_C_by_D = hvcat(M, permutedims([(X.*U[:,m])'*W*(Z[:,l].*V[:,m]) for l=1:L, m=1:M])...)

        F_U_by_V = hvcat(J, permutedims([W[i,j]*(V[j,:].*D)*(U[i,:].*D)' for i=1:I, j=1:J])...)
        F_U_by_D = vcat([(V.*D')'*(W[i,:].*V.*U[i,:]') for i=1:I]...)
        F_V_by_D = vcat([(U.*D')'*(W[:,j].*U.*V[j,:]') for j=1:J]...)

        Q,Q_A,Q_B,Q_U,Q_V = compute_constraint_Jacobian(X,Z,A,B,U,V)

        F_ABC = [F_A_full  F_A_by_B  F_A_by_C
                 F_A_by_B' F_B_full  F_B_by_C
                 F_A_by_C' F_B_by_C' F_C_full]
        F_UVD = [F_U_full  F_U_by_V  F_U_by_D
                 F_U_by_V' F_V_full  F_V_by_D
                 F_U_by_D' F_V_by_D' F_D_full]
        F_ABC_by_UVD = [F_A_by_U F_A_by_V F_A_by_D
                        F_B_by_U F_B_by_V F_B_by_D
                        F_C_by_U F_C_by_V F_C_by_D]
        F_ABCUDV = [F_ABC         F_ABC_by_UVD
                    F_ABC_by_UVD' F_UVD       ]
        N = size(Q,2)
        F_full = [F_ABCUDV  Q
                  Q'        zeros(N,N)]
        F_full_inv = inv(F_full)  # todo: Should we be using the priors here?

        S_A = 1:J*K
        S_B = (1:I*L).+maximum(S_A)
        S_C = (1:K*L).+maximum(S_B)
        S_U = (1:I*M).+maximum(S_C)
        S_V = (1:J*M).+maximum(S_U)
        S_D = (1:M).+maximum(S_V)
        S_Q = (1:N).+maximum(S_D)

        var_A_approx = diag(F_full_inv)[S_A]
        var_B_approx = diag(F_full_inv)[S_B]
        var_C_approx = diag(F_full_inv)[S_C]
        var_U_approx = diag(F_full_inv)[S_U]
        var_V_approx = diag(F_full_inv)[S_V]
        var_D_approx = diag(F_full_inv)[S_D]

    end

    se_A = sqrt.(max.(eps(0.0),reshape(var_A_approx, K,J))')
    se_B = sqrt.(max.(eps(0.0),reshape(var_B_approx, L,I))')
    se_C = sqrt.(max.(eps(0.0),reshape(var_C_approx, K,L)))
    se_U = sqrt.(max.(eps(0.0),reshape(var_U_approx, M,I))')
    se_V = sqrt.(max.(eps(0.0),reshape(var_V_approx, M,J))')
    se_D = sqrt.(max.(eps(0.0),var_D_approx))

    # Delta propagation for S and T
    D1,D2 = compute_logdispersion_derivatives(Y,Mu,r)
    g_S = vec(sum(D1; dims=2)) .- p.lambda_S*(S .- p.mu_S)
    g_T = vec(sum(D1; dims=1)) .- p.lambda_T*(T .- p.mu_T)
    h_S = vec(sum(D2; dims=2)) .- p.lambda_S
    h_T = vec(sum(D2; dims=1)) .- p.lambda_T

    # Compute observed information for S and T (use observed information since expected info can't be analytically computed)
    F_S = [(h_S[i]<=0 ? -h_S[i] : NaN) for i = 1:I]
    F_T = [(h_T[j]<=0 ? -h_T[j] : NaN) for j = 1:J]

    var_S_from_U,var_S_from_V = compute_var_S_from_UV(Y,Mu,W,E,r,U,V,D,g_S,h_S,se_U,se_V)
    var_S_from_B,var_S_from_A = compute_var_S_from_UV(Y,Mu,W,E,r,X,Z,1,g_S,h_S,se_B,se_A)
    var_T_from_U,var_T_from_V = compute_var_T_from_UV(Y,Mu,W,E,r,U,V,D,g_T,h_T,se_U,se_V)
    var_T_from_B,var_T_from_A = compute_var_T_from_UV(Y,Mu,W,E,r,X,Z,1,g_T,h_T,se_B,se_A)

    # Calculate standard errors for S and T
    se_S = sqrt.(1.0./F_S + var_S_from_U + var_S_from_V + var_S_from_A + var_S_from_B)
    se_T = sqrt.(1.0./F_T + var_T_from_U + var_T_from_V + var_T_from_A + var_T_from_B)

    # Original approach: Just use F_S and F_T
    # h_S = vec(sum(D2; dims=2)) .- p.lambda_S
    # D1,D2 = compute_logdispersion_derivatives(Y,Mu,r)
    # h_T = vec(sum(D2; dims=1)) .- p.lambda_T
    # F_S = [(h_S[i]<=0 ? -h_S[i] : NaN) for i = 1:I]
    # F_T = [(h_T[j]<=0 ? -h_T[j] : NaN) for j = 1:J]
    # se_S = sqrt.(1.0./F_S)
    # se_T = sqrt.(1.0./F_T)

    if any(isnan.(se_A)); @warn("One or more invalid standard errors in se_A."); end
    if any(isnan.(se_B)); @warn("One or more invalid standard errors in se_B."); end
    if any(isnan.(se_C)); @warn("One or more invalid standard errors in se_C."); end
    if any(isnan.(se_D)); @warn("One or more invalid standard errors in se_D."); end
    if any(isnan.(se_U)); @warn("One or more invalid standard errors in se_U."); end
    if any(isnan.(se_V)); @warn("One or more invalid standard errors in se_V."); end
    if any(isnan.(se_S)); @warn("One or more invalid standard errors in se_S."); end
    if any(isnan.(se_T)); @warn("One or more invalid standard errors in se_T."); end

    return se_A,se_B,se_C,se_U,se_V,se_S,se_T
end

# __________________________________________________________________________________________________
# Functions for inference in NB-GBM

# Compute X'*diag(w)*X and store it in XtWX
function compute_XtWX!(XtWX,X,w)
    I,K = size(X)
    for k1 = 1:K, k2 = 1:K
        value = 0.0
        for i = 1:I
            value += X[i,k1]*w[i]*X[i,k2]
        end
        XtWX[k1,k2] = value
    end
end

function compute_F_D(U,V,W)
    I,M = size(U)
    J,M = size(V)
    F_D = zeros(M,M)
    for m1 = 1:M, m2 = 1:M
        value = 0.0
        for j = 1:J
            uw = 0.0
            for i = 1:I
                uw += U[i,m1]*U[i,m2]*W[i,j]
            end
            value += uw*V[j,m1]*V[j,m2]
        end
        F_D[m1,m2] = value
    end
    return F_D
end

# Jacobian (transposed) of identifiability constraints on U and V
# (i.e., X'*U=0, Z'*V=0, U'*U=I, and V'*V=I)
function compute_constraint_Jacobian_UV(X,Z,U,V)
    I,M = size(U)
    J,M = size(V)

    Q_U_x = kron(X, Diagonal(ones(M)))  # Constraint: vec(U'*X)=0
    Q_V_z = kron(Z, Diagonal(ones(M)))  # Constraint: vec(V'*Z)=0

    # Orthogonality contraints on U and V
    Nc = nchoose2(M+1)
    Q_U_u = zeros(I*M,Nc)
    Q_V_v = zeros(J*M,Nc)
    for m = 1:M
        cols = [nchoose2(m).+(1:m) ; nchoose2.(m.+(1:(M-m))).+m]
        Q_U_u[(1:M:I*M).+(m-1), cols] = U
        Q_V_v[(1:M:J*M).+(m-1), cols] = V
        # Note: the factor of 2 for the m'th column can be ignored since the contraints are equivalent.
    end

    Q_U = [Q_U_u Q_U_x]
    Q_V = [Q_V_v Q_V_z]
    return Q_U,Q_V
end

function compute_var_C_from_A(X,Z,F_A_inv,F_C_inv,grad_C,dW_dMu_Mu,dE_dMu_Mu)
    I,K = size(X)
    J,L = size(Z)
    dC_dA = zeros(K*L,J*K)  # Jacobian
    dW_da = zeros(I)
    dE_da = zeros(I)
    ZjZjt = zeros(L,L)
    XtWX = zeros(K,K)
    for j = 1:J
        for l1 = 1:L, l2 = 1:L
            ZjZjt[l1,l2] = Z[j,l1]*Z[j,l2]
        end
        for k = 1:K
            for i = 1:I
                dW_da[i] = dW_dMu_Mu[i,j]*X[i,k]
                dE_da[i] = dE_dMu_Mu[i,j]*X[i,k]
            end
            compute_XtWX!(XtWX,X,dW_da)  # Compute XtWX = X'*(dW_da.*X)
            dF_da = kron(ZjZjt, XtWX)
            dC_dA[:,(j-1)*K+k] = F_C_inv*(-dF_da*(F_C_inv*grad_C) + vec((X'*dE_da)*Z[j,:]'))
        end
    end
    F_A_inv_dC_dAt = zeros(J*K,K*L)
    for j = 1:J
        F_A_inv_dC_dAt[block(j,K),:] = F_A_inv[j]*dC_dA[:,((j-1)*K+1):(j*K)]'
    end
    return vec(sum(dC_dA'.*F_A_inv_dC_dAt,dims=1))
end

function compute_var_C_from_B(X,Z,F_B_inv,F_C_inv,grad_C,dW_dMu_Mu,dE_dMu_Mu)
    I,K = size(X)
    J,L = size(Z)
    dC_dB = zeros(K*L,I*L)  # Jacobian
    dW_db = zeros(J)
    dE_db = zeros(J)
    XiXit = zeros(K,K)
    ZtWZ = zeros(L,L)
    for i = 1:I
        for k1 = 1:K, k2 = 1:K
            XiXit[k1,k2] = X[i,k1]*X[i,k2]
        end
        for l = 1:L
            for j = 1:J
                dW_db[j] = dW_dMu_Mu[i,j]*Z[j,l]
                dE_db[j] = dE_dMu_Mu[i,j]*Z[j,l]
            end
            compute_XtWX!(ZtWZ,Z,dW_db)  # Compute ZtWZ = Z'*(dW_db.*Z)
            dF_db = kron(ZtWZ, XiXit)
            dC_dB[:,(i-1)*L+l] = F_C_inv*(-dF_db*(F_C_inv*grad_C) + vec(X[i,:]*(dE_db'*Z)))
        end
    end
    F_B_inv_dC_dBt = zeros(I*L,K*L)
    for i = 1:I
        F_B_inv_dC_dBt[block(i,L),:] = F_B_inv[i]*dC_dB[:,((i-1)*L+1):(i*L)]'
    end
    return vec(sum(dC_dB'.*F_B_inv_dC_dBt,dims=1))
end

function compute_var_A_from_V(X,U,D,var_V,F_A_inv,grad_A,dW_dMu_Mu,dE_dMu_Mu)
    J,K = size(grad_A)
    I,M = size(U)
    var_A_from_V = zeros(J*K)
    dA_dV = zeros(K,M)
    dW_dv = zeros(I)
    dE_dv = zeros(I)
    XtWX = zeros(K,K)
    for j = 1:J
        for m = 1:M
            for i = 1:I
                dW_dv[i] = dW_dMu_Mu[i,j]*U[i,m]*D[m]
                dE_dv[i] = dE_dMu_Mu[i,j]*U[i,m]*D[m]
            end
            compute_XtWX!(XtWX,X,dW_dv)  # Compute XtWX = X'*(dW_dv.*X)
            dA_dV[:,m] = (-F_A_inv[j]*XtWX)*(F_A_inv[j]*grad_A[j,:]) + F_A_inv[j]*(X'*dE_dv)
        end
        var_A_from_V[block(j,K)] = sum(dA_dV'.*(var_V[block(j,M)].*dA_dV'), dims=1)
    end
    return var_A_from_V
end

function compute_var_A_from_U(X,V,D,var_U,F_A_inv,grad_A,dW_dMu_Mu,dE_dMu_Mu)
    I,K = size(X)
    J,M = size(V)
    dA_dU = zeros(J*K,I*M)  # Jacobian
    for j = 1:J
        Temp = -F_A_inv[j]*(X.*dW_dMu_Mu[:,j])' .* (X*F_A_inv[j]*grad_A[j,:])' + F_A_inv[j]*(X.*dE_dMu_Mu[:,j])'
        for m = 1:M
            dA_dU[block(j,K), (1:M:I*M).+(m-1)] = Temp*V[j,m]*D[m]
        end
    end
    return vec(sum(dA_dU'.*(var_U.*dA_dU'),dims=1))
end

function compute_var_B_from_U(Z,V,D,var_U,F_B_inv,grad_B,dW_dMu_Mu,dE_dMu_Mu)
    return compute_var_A_from_V(Z,V,D,var_U,F_B_inv,grad_B,dW_dMu_Mu',dE_dMu_Mu')
end

function compute_var_B_from_V(Z,U,D,var_V,F_B_inv,grad_B,dW_dMu_Mu,dE_dMu_Mu)
    return compute_var_A_from_U(Z,U,D,var_V,F_B_inv,grad_B,dW_dMu_Mu',dE_dMu_Mu')  # <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< double-check that this is correct
end

function compute_var_V_from_A(UD,X,var_A,F_V_inv,grad_V,dW_dMu_Mu,dE_dMu_Mu)
    return compute_var_A_from_V(UD,X,ones(size(X,2)),var_A,F_V_inv,grad_V,dW_dMu_Mu,dE_dMu_Mu)
end

function compute_var_V_from_B(UD,Z,var_B,F_V_inv,grad_V,dW_dMu_Mu,dE_dMu_Mu)
    return compute_var_A_from_U(UD,Z,ones(size(Z,2)),var_B,F_V_inv,grad_V,dW_dMu_Mu,dE_dMu_Mu)
end

function compute_var_U_from_B(VD,Z,var_B,F_U_inv,grad_U,dW_dMu_Mu,dE_dMu_Mu)
    return compute_var_B_from_U(VD,Z,ones(size(Z,2)),var_B,F_U_inv,grad_U,dW_dMu_Mu,dE_dMu_Mu)
end

function compute_var_S_from_UV(Y,Mu,W,E,r,U,V,D,g_S,h_S,se_U,se_V)
    Q = -W.*E./r
    R = -2*W.*E./(r.+Mu)

    # Delta propagation to account for uncertainty in S due to U
    dg_du = Q*(V.*D')
    dh_du = -dg_du + R*(V.*D')
    ds_du = -dg_du./h_S + dh_du.*(g_S./h_S.^2)
    var_S_from_U = vec(sum(ds_du.*(ds_du.*(se_U.^2)); dims=2))

    # Delta propagation to account for uncertainty in S due to V
    I = size(U,1)
    var_S_from_V = zeros(I)
    for i = 1:I
        dg_dv = Q[i,:].*(U[i,:].*D)'
        dh_dv = -dg_dv + R[i,:].*(U[i,:].*D)'
        ds_dv = -dg_dv/h_S[i] + dh_dv*(g_S[i]/h_S[i]^2)
        var_S_from_V[i] = sum(ds_dv.*(ds_dv.*(se_V.^2)))
    end

    return var_S_from_U,var_S_from_V
end

function compute_var_T_from_UV(Y,Mu,W,E,r,U,V,D,g_T,h_T,se_U,se_V)
     var_T_from_V,var_T_from_U = compute_var_S_from_UV(Y',Mu',W',E',r',V,U,D,g_T,h_T,se_V,se_U)
     return var_T_from_U,var_T_from_V
end

# _________________________________________________________________________________________________
# The following function (compute_constraint_Jacobian) is used only in the full Fisher information matrix version.

# Jacobian (transposed) of identifiability constraints on A, B, U, and V
# (i.e., Z'*A=0, X'*B=0, X'*U=0, Z'*V=0, U'*U=I, and V'*V=I)
function compute_constraint_Jacobian(X,Z,A,B,U,V)
    I,K = size(X)
    J,L = size(Z)
    M = size(U,2)

    # Constraints on A and B
    Q_A = kron(Z, Diagonal(ones(K)))  # if we combine with F_A, these might fall into the same block diagonal structure as F_A
    Q_B = kron(X, Diagonal(ones(L)))  # ... F_B ...

    Q_U,Q_V = compute_constraint_Jacobian_UV(X,Z,U,V)

    S_A = 1:J*K
    S_B = (1:I*L).+maximum(S_A)
    S_C = (1:K*L).+maximum(S_B)
    S_U = (1:I*M).+maximum(S_C)
    S_V = (1:J*M).+maximum(S_U)
    S_D = (1:M).+maximum(S_V)

    Q = zeros(maximum(S_D), K*L+K*L+K*M+Nc+L*M+Nc)
    Q[S_A,1:K*L] = Q_A
    Q[S_B,(1:K*L).+K*L] = Q_B
    Q[S_U,(1:(K*M+Nc)).+K*L.+K*L] = Q_U
    Q[S_V,(1:(L*M+Nc)).+K*L.+K*L.+K*M.+Nc] = Q_V

    return Q,Q_A,Q_B,Q_U,Q_V
end

end # module Infer