text:
### A Pluto.jl notebook ###
# v0.12.4

using Markdown
using InteractiveUtils

# ╔═╡ d006b66c-0ea6-11eb-0e41-a37bc0f2689e
using Pkg

# ╔═╡ f33f0d64-0eb0-11eb-13ce-5dcce414b595
using TextAnalysis

# ╔═╡ 761f51f0-0eb7-11eb-256e-c92b976c2337
using Plots

# ╔═╡ b1cea02c-0ebf-11eb-0bab-4747d1edf62d
using Genie

# ╔═╡ acc93070-0ec0-11eb-1cff-e337990c9b6b
using HTTP

# ╔═╡ d849e946-0ec1-11eb-1606-25ff1efe4233
using LightGraphs

# ╔═╡ 14f119cc-0ec3-11eb-3bca-994fb5f6b123
using GraphPlot

# ╔═╡ febda3f6-0ec1-11eb-02d1-05a4caaddc82
using SimpleWeightedGraphs

# ╔═╡ ab229316-0ec4-11eb-0b00-a756c44cd0f2
using Franklin

# ╔═╡ 53a6a35a-0ec9-11eb-3aec-e152cac5c4e4
using Images, FileIO

# ╔═╡ 378c94c0-0e32-11eb-2781-898c22ee1f9a
html"<button onclick=present()>Present</button>"

# ╔═╡ 13eb6dec-0e33-11eb-3895-0925c4c233b0
md"# Introduction to Julia"

# ╔═╡ ada1bad6-0e9c-11eb-060b-13019859ffef
md"""# The problem

- specialists work on algorithms (python, r, matlab)

- programmers make it production ready (c++, c#, java)

- deploy
"""

# ╔═╡ 48c5b260-0e9d-11eb-16b9-6fe39c17024a
md"""# The solution

- specialists and programmers work collaborate (julia)
- deploy

"""

# ╔═╡ 4f6b8102-0e9e-11eb-3af1-d36898b19ae5
md"""# Come for the syntax, stay for the speed"""

# ╔═╡ f6e47e10-0ea4-11eb-13c6-8dc21748213a
md"""# Flux.jl

- julia replacement for pytorch, tensorflow
- examples available at https://github.com/FluxML/model-zoo
"""

# ╔═╡ b5cbc9a8-0ea9-11eb-30de-9b7b32ac35a4
Pkg.activate(".")

# ╔═╡ d76b1ba2-0ea9-11eb-19b4-af7445adb3e2
Pkg.instantiate()

# ╔═╡ e014f00c-0ea9-11eb-3cfe-15ec3e904b16
md"""
Using Flux is as simple as doing:

```jl
] add Flux

using Flux
```
"""

# ╔═╡ ff6a8e58-0eae-11eb-14cd-85194adc6cb3
md"""# TextAnalysis.jl
"""

# ╔═╡ 27f8d1a2-0eb1-11eb-14db-955bc45260cb
str = "this is a new string"

# ╔═╡ 7156320e-0eb1-11eb-34a7-3bcd88f43ba7
sd = StringDocument(str)

# ╔═╡ c4a6a524-0eb1-11eb-079e-d5d65072a92b
tokens(sd)

# ╔═╡ cc093e80-0eb1-11eb-1202-0510253dadc2
ngrams(sd)

# ╔═╡ d4563f52-0eb1-11eb-1baf-f591153c1edb
ngrams(sd, 2)

# ╔═╡ fb74f386-0eb5-11eb-3e63-099045ccbd8c
md"""# """

# ╔═╡ fb81ac10-0eb1-11eb-26a4-e9eb81a04204
str_2 = "this is a bad string !( some thing here*^"

# ╔═╡ 175ace94-0eb2-11eb-134f-6d7b29122634
sd_2 = StringDocument(str_2)

# ╔═╡ af7469a6-0eb2-11eb-017c-51456043a04f
prepare!(sd_2, strip_punctuation)

# ╔═╡ b9a6ccac-0eb2-11eb-255f-fdf6a9d3a6fa
text(sd_2)

# ╔═╡ b725994a-0eb2-11eb-1a93-2126fe4ca0b1
md"""# Plots.jl"""

# ╔═╡ 7abc8296-0eb1-11eb-2143-299bb5273810
x = 1:10; y = rand(10)

# ╔═╡ f5600526-0eba-11eb-243f-ffe32c8a4c1e
plot(x, y)

# ╔═╡ f558d724-0ef1-11eb-3777-3f417a72b58d
md"# "

# ╔═╡ fb6f9b84-0eba-11eb-07c0-0d85d8377eb5
plot(Plots.fakedata(50, 5), w = 2)

# ╔═╡ a9151da2-0ebb-11eb-238b-4de8054a23c9
theme(:solarized)

# ╔═╡ 335968c6-0ebc-11eb-310f-31a13429441d
md"""# DifferentialEquations.jl

- solve ODEs through julia!
- one of the most popular julia packages
"""

# ╔═╡ 749ac066-0ebf-11eb-3f96-5d670cf38d3a
md"""# Genie.jl

- full-stack mvc web framework in julia
"""

# ╔═╡ c5d3aa24-0ebf-11eb-3947-054ca8659285
route("/hello") do
	"Hello, world!"
end

# ╔═╡ 92a9a896-0ec0-11eb-3bfb-6178d10f49c0
up()

# ╔═╡ b2478a42-0ec0-11eb-2f93-774dd80dc210
HTTP.get("http://localhost:8000/hello")

# ╔═╡ c8286982-0ec0-11eb-32ff-9183ffab6633
down()

# ╔═╡ e8f7a2b8-0ec0-11eb-1136-f520768cc6e8
md"""# LightGraphs.jl

- provides a framework for graph framrworks like SimpleWeightedGraphs.jl
"""

# ╔═╡ e54f8498-0ec1-11eb-1756-6f29b7d0897e
g_1 = path_graph(6)

# ╔═╡ ed0f5cfc-0ec2-11eb-2e68-71834cf8fe1d
gplot(g_1)

# ╔═╡ 065fc5a8-0ef2-11eb-326d-a706ed46bff2
md"# "

# ╔═╡ 0758d3f0-0ec2-11eb-1cce-dd57ef093259
g_2 = SimpleWeightedGraph(10)

# ╔═╡ 27cd3d5a-0ec3-11eb-0830-8f08536cd744
gplot(g_2)

# ╔═╡ 304922fa-0ec3-11eb-04c5-23b1e8f04957
md"""# Documenter.jl

- generate documentation for your julia packages quickly

```
using Documenter, Example

makedocs(sitename="My documentation")
```
"""

# ╔═╡ a2f8da16-0ec3-11eb-2d75-2d996650a7a8
md"""# Franklin.jl

- static site generator in julia
"""

# ╔═╡ be3d1dae-0ec4-11eb-3712-1f370b47457b
newsite("My new site", template="pure-sm")

# ╔═╡ 486c8028-0ec5-11eb-0ee3-d16eaf57ae77
readdir()

# ╔═╡ bcb6705c-0ef1-11eb-1ab2-93c4095bdd37
md"""# Pluto.jl

### This presentation is made with Pluto.jl!
"""

# ╔═╡ e66f3dc4-0ec5-11eb-18bd-05905909377a
md"""# Hush, there are a lot of these

- checkout https://juliaobserver.com/packages
- or just do a quick google or github search
- join slack and ask on #helpdesk
- everyone is very active and ready to help
"""

# ╔═╡ 20bbfee0-0ec6-11eb-2503-7d18df89f6d0
md"""# Julia, is it better than python?

- judge it yourselves
- here are some stats for you guys
"""

# ╔═╡ a3992086-0ec6-11eb-3e22-3745ab03e956
load("./images/picture_1.png")

# ╔═╡ 4926e5d8-0ef2-11eb-2550-bd608cb98950
md"# "

# ╔═╡ aab7782c-0ecb-11eb-2d06-8b1ee9d60b5a
load("./images/picture_2.png")

# ╔═╡ 3554232c-0ecc-11eb-1929-47354e855959
md"""# """

# ╔═╡ 220ab8bc-0ecc-11eb-3c41-33b1e564486a
load("./images/picture_3.png")

# ╔═╡ 53cc2dd8-0ef2-11eb-2b79-8572c61d90e9
md"# "

# ╔═╡ 24f89fee-0ecc-11eb-1525-cddda435e969
load("./images/picture_4.png")

# ╔═╡ 4225747a-0ecc-11eb-19fc-61c792a7b733
md"""# """

# ╔═╡ 28926e96-0ecc-11eb-215c-61b582c64036
load("./images/picture_5.png")

# ╔═╡ 2baf6408-0ecc-11eb-22e1-d1aed1102d96
md"""# Final words

- once should def look at julia as an alternative to their tasks
- it has come very reliable and fast and is expected to grow a lot more in future
- is a perfect solution if you work in data science field
- has a lot more than what I summarized in these very few slides
- enjoy and hope you liked this mini presentation
"""

# ╔═╡ 12b49df0-0ecd-11eb-0ad7-592765b2578b
md"""# Any questions?

ps: please don't ask me tough ones :P
"""

# ╔═╡ 26aa052c-0ecd-11eb-1df6-db2143e34471
println("thank you!")

# ╔═╡ Cell order:
# ╟─378c94c0-0e32-11eb-2781-898c22ee1f9a
# ╟─13eb6dec-0e33-11eb-3895-0925c4c233b0
# ╟─ada1bad6-0e9c-11eb-060b-13019859ffef
# ╟─48c5b260-0e9d-11eb-16b9-6fe39c17024a
# ╟─4f6b8102-0e9e-11eb-3af1-d36898b19ae5
# ╟─f6e47e10-0ea4-11eb-13c6-8dc21748213a
# ╠═d006b66c-0ea6-11eb-0e41-a37bc0f2689e
# ╠═b5cbc9a8-0ea9-11eb-30de-9b7b32ac35a4
# ╠═d76b1ba2-0ea9-11eb-19b4-af7445adb3e2
# ╟─e014f00c-0ea9-11eb-3cfe-15ec3e904b16
# ╟─ff6a8e58-0eae-11eb-14cd-85194adc6cb3
# ╠═f33f0d64-0eb0-11eb-13ce-5dcce414b595
# ╠═27f8d1a2-0eb1-11eb-14db-955bc45260cb
# ╠═7156320e-0eb1-11eb-34a7-3bcd88f43ba7
# ╠═c4a6a524-0eb1-11eb-079e-d5d65072a92b
# ╠═cc093e80-0eb1-11eb-1202-0510253dadc2
# ╠═d4563f52-0eb1-11eb-1baf-f591153c1edb
# ╟─fb74f386-0eb5-11eb-3e63-099045ccbd8c
# ╠═fb81ac10-0eb1-11eb-26a4-e9eb81a04204
# ╠═175ace94-0eb2-11eb-134f-6d7b29122634
# ╠═af7469a6-0eb2-11eb-017c-51456043a04f
# ╠═b9a6ccac-0eb2-11eb-255f-fdf6a9d3a6fa
# ╟─b725994a-0eb2-11eb-1a93-2126fe4ca0b1
# ╠═761f51f0-0eb7-11eb-256e-c92b976c2337
# ╠═7abc8296-0eb1-11eb-2143-299bb5273810
# ╠═f5600526-0eba-11eb-243f-ffe32c8a4c1e
# ╟─f558d724-0ef1-11eb-3777-3f417a72b58d
# ╠═fb6f9b84-0eba-11eb-07c0-0d85d8377eb5
# ╠═a9151da2-0ebb-11eb-238b-4de8054a23c9
# ╟─335968c6-0ebc-11eb-310f-31a13429441d
# ╟─749ac066-0ebf-11eb-3f96-5d670cf38d3a
# ╠═b1cea02c-0ebf-11eb-0bab-4747d1edf62d
# ╠═c5d3aa24-0ebf-11eb-3947-054ca8659285
# ╠═92a9a896-0ec0-11eb-3bfb-6178d10f49c0
# ╠═acc93070-0ec0-11eb-1cff-e337990c9b6b
# ╠═b2478a42-0ec0-11eb-2f93-774dd80dc210
# ╠═c8286982-0ec0-11eb-32ff-9183ffab6633
# ╟─e8f7a2b8-0ec0-11eb-1136-f520768cc6e8
# ╠═d849e946-0ec1-11eb-1606-25ff1efe4233
# ╠═e54f8498-0ec1-11eb-1756-6f29b7d0897e
# ╠═14f119cc-0ec3-11eb-3bca-994fb5f6b123
# ╠═ed0f5cfc-0ec2-11eb-2e68-71834cf8fe1d
# ╟─065fc5a8-0ef2-11eb-326d-a706ed46bff2
# ╠═febda3f6-0ec1-11eb-02d1-05a4caaddc82
# ╠═0758d3f0-0ec2-11eb-1cce-dd57ef093259
# ╠═27cd3d5a-0ec3-11eb-0830-8f08536cd744
# ╟─304922fa-0ec3-11eb-04c5-23b1e8f04957
# ╟─a2f8da16-0ec3-11eb-2d75-2d996650a7a8
# ╠═ab229316-0ec4-11eb-0b00-a756c44cd0f2
# ╠═be3d1dae-0ec4-11eb-3712-1f370b47457b
# ╠═486c8028-0ec5-11eb-0ee3-d16eaf57ae77
# ╟─bcb6705c-0ef1-11eb-1ab2-93c4095bdd37
# ╟─e66f3dc4-0ec5-11eb-18bd-05905909377a
# ╟─20bbfee0-0ec6-11eb-2503-7d18df89f6d0
# ╟─53a6a35a-0ec9-11eb-3aec-e152cac5c4e4
# ╟─a3992086-0ec6-11eb-3e22-3745ab03e956
# ╟─4926e5d8-0ef2-11eb-2550-bd608cb98950
# ╟─aab7782c-0ecb-11eb-2d06-8b1ee9d60b5a
# ╟─3554232c-0ecc-11eb-1929-47354e855959
# ╟─220ab8bc-0ecc-11eb-3c41-33b1e564486a
# ╟─53cc2dd8-0ef2-11eb-2b79-8572c61d90e9
# ╟─24f89fee-0ecc-11eb-1525-cddda435e969
# ╟─4225747a-0ecc-11eb-19fc-61c792a7b733
# ╟─28926e96-0ecc-11eb-215c-61b582c64036
# ╟─2baf6408-0ecc-11eb-22e1-d1aed1102d96
# ╟─12b49df0-0ecd-11eb-0ad7-592765b2578b
# ╠═26aa052c-0ecd-11eb-1df6-db2143e34471


text:
# Main module for BayesianMixtures package
module BayesianMixtures

include("MFM.jl")
include("RandomNumbers.jl")

include("Normal.jl")
include("MVN.jl")
include("MVNaaC.jl")
include("MVNaaN.jl")
include("MVNaaRJ.jl")
include("NormalNonoptimized.jl")

using Random
using Statistics
using SpecialFunctions
lgamma_(x) = logabsgamma(x)[1]

# ===================================================================
# ===================================================================
# ================== Functions to generate results ==================
# ===================================================================
# ===================================================================

# Create an options object to specify model, data, and MCMC parameters.
function options(
        mode, # "Normal", "MVN", "MVNaaC", "MVNaaN", or "MVNaaRJ"
        model_type, # "MFM" or "DPM"
        x, # data
        n_total; # total number of MCMC sweeps to run the sampler
        n_keep=n_total, # number of MCMC sweeps to keep after thinning
        n_burn=round(Int,n_total/10), # number of MCMC sweeps (out of n_total) to discard as burn-in
        verbose=true, # display information or not
        use_hyperprior=true, # put hyperprior on the base distn parameters or not
        t_max=40, # a guess at an upper bound on # of clusters that will be encountered during MCMC

        # MFM options:
        gamma=1.0, # Dirichlet_k(gamma,...,gamma)
        log_pk="k -> log(0.1)+(k-1)*log(0.9)", # string representation of log(p(k))
            # (where p(k) is the log of the prior on # of components, K)

        # DPM options:
        alpha_random=true, # put prior on alpha (DPM concentration parameter) or not
        alpha=1.0, # value of alpha (initial value if alpha_random=true)

        # Jain-Neal split-merge options:
        use_splitmerge=true, # use split-merge or not
        n_split=5, # number of intermediate sweeps for split launch state
        n_merge=5,  #                 "         "       merge    "     "

        # RJMCMC options:
        k_max=t_max # a guess at an upper bound on # of components that will be encountered during MCMC
    )

    # Compute partition distribution values
    n = length(x)
    if model_type=="MFM"
        lpk = eval(Meta.parse(log_pk))
        log_pk_fn(k) = Base.invokelatest(lpk,k)
        log_v = MFM.coefficients(log_pk_fn,gamma,n,t_max+1)
        a = b = gamma
    elseif model_type=="DPM"
        log_v = float(1:t_max+1)*log(alpha) .- lgamma_(alpha+n) .+ lgamma_(alpha)
        a,b = 1.,0.
    else
        error("Invalid model_type: $model_type.")
    end
    if mode=="MVNaaRJ"; @assert(model_type=="MFM", "RJMCMC is not implemented for DPMs."); end

    n_keep = min(n_keep,n_total)
    module_ = getfield(BayesianMixtures,Symbol(mode))
    return module_.Options(mode, model_type, x, n_total, n_keep, n_burn, verbose,
                           use_hyperprior, t_max, gamma, log_pk, alpha_random, alpha,
                           use_splitmerge, n_split, n_merge, k_max, a, b, log_v, n)
end

# Run the MCMC sampler with the specified options.
function run_sampler(options)
    o = options
    n,n_total,n_keep = o.n,o.n_total,o.n_keep
    module_ = getfield(BayesianMixtures,Symbol(o.mode))

    # Short run to precompile (calling precompile doesn't seem to work...)
    module_.sampler(o,1,1)

    if o.verbose
        println(o.mode, " ", o.model_type)
        println("n = $n, n_total = $n_total, n_keep = $n_keep")
        print("Running... ")
    end

    # Main run
    elapsed_time = (@elapsed t_r,N_r,z_r,theta_r,keepers = module_.sampler(o,n_total,n_keep))
    time_per_step = elapsed_time/(n_total*n)

    if o.verbose
        println("complete.")
        println("Elapsed time = $elapsed_time seconds")
        println("Time per step ~ $time_per_step seconds")
    end

    # # Run profiler
    # @profile sampler(x,n_total,n_keep,o)
    # Profile.print(format = :flat)
    # Profile.clear()

    return module_.Result(o,t_r,N_r,z_r,theta_r,keepers,elapsed_time,time_per_step)
end

# ===================================================================
# ===================================================================
# ================== Functions to analyze results ===================
# ===================================================================
# ===================================================================

using Pkg
packages_installed = [pkg.name for pkg in collect(values(Pkg.dependencies()))]
can_save = ("JLD" in packages_installed)
if can_save; using JLD; end

# Save result to file.
function save_result(filename,result)
    @assert(can_save, "Results cannot be saved to file since JLD is not installed.")
    JLD.save(filename,"result",result)
end

# Load result from file.
function load_result(filename)
    @assert(can_save, "Results cannot be loaded from file since JLD is not installed.")
    return JLD.load(filename,"result")
end

# Compute running average of posterior CDF on t (# of clusters)
function t_running(result)
    o = result.options
    n_total = o.n_total
    T = zeros(n_total,o.t_max)
    for i = 1:n_total; T[i,result.t[i]] = 1; end
    cdfs = cumsum(cumsum(T,dims=1),dims=2)./(1:n_total)
    return cdfs
end

# Compute histogram with the specified bin edges,
# where x[i] is in bin j if edges[j] < x[i] <= edges[j+1].
function histogram(x, edges=[]; n_bins=50, weights=ones(length(x)))
    if isempty(edges)
        mn,mx = minimum(x),maximum(x)
        r = mx-mn
        edges = range(mn-r/n_bins, stop=mx+r/n_bins, length=n_bins)
    else
        n_bins = length(edges)-1
    end
    counts = zeros(Float64,n_bins)
    for i=1:length(x)
        for j=1:n_bins
            if (edges[j] < x[i] <= edges[j+1])
                counts[j] += weights[i]
                break
            end
        end
    end
    return counts,edges
end

# Compute posterior on t (# of clusters)
function t_posterior(result)
    o = result.options
    use = o.n_burn+1:o.n_total
    counts,edges = histogram(result.t[use], 0:o.t_max) # bins are: (0,1],(1,2],...,(t_max-1,t_max]
    return counts/length(use)
end

# Compute posterior on k (# of components)
function k_posterior(result; upto=result.options.t_max)
    o = result.options
    @assert(o.model_type=="MFM", "The posterior on k is not defined for the DPM.")
    lpk = eval(Meta.parse(o.log_pk))
    log_pk_fn(k) = Base.invokelatest(lpk,k)
    p_kt = MFM.p_kt(log_pk_fn,o.gamma,o.n,upto,o.t_max)
    return p_kt*t_posterior(result)
end

# Compute density estimate at a point x (using Green's pseudo-predictive density).
function density_estimate(x,result)
    o = result.options
    n,b = o.n,o.b
    t,N,theta,keepers = result.t, result.N, result.theta, result.keepers
    n_keep = length(keepers)
    loglik = getfield(BayesianMixtures,Symbol(o.mode)).log_likelihood
    r = 0.0
    n_used = 0
    for ik = 1:n_keep
        i = keepers[ik]
        if i>o.n_burn
            for c in findall(N[:,i].>0)
                r += exp(loglik(x,theta[c,ik])) * (N[c,i]+b)/(n + b*t[i])
            end
            n_used += 1
        end
    end
    return r/n_used
end

rand_categorical(p) = (s=rand(); i=0; while s>0; i+=1; s-=p[i]; end; i)
rand_dirichlet(a) = (w=RandomNumbers.gamma.(a,1); w/sum(w))
counts(a,k) = (N=zeros(Int,k); for ai in a; N[ai]+=1; end; N)
logsumexp(x) = (m = maximum(x); m == -Inf ? -Inf : log.(sum(exp.(x.-m))) + m)

# Convert MFM result into posterior samples of k,w,a,theta,
# where k = number of components
#       w = mixture weights
#       a = assignments of data points to mixture components
#       theta = mixture component parameters
function sample_mixture_parameters(result,kmax)
    r,o = result,result.options
    @assert(o.model_type=="MFM", "This function is not defined for the DPM.")
    lpk = eval(Meta.parse(o.log_pk))
    log_pk_fn(k) = Base.invokelatest(lpk,k)
    p_kt = MFM.p_kt(log_pk_fn,o.gamma,o.n,kmax,o.t_max)

    module_ = getfield(BayesianMixtures,Symbol(o.mode))
    Theta = module_.Theta
    H = module_.construct_hyperparameters(o)

    k = zeros(Int16,o.n_keep)
    a = zeros(Int16,o.n,o.n_keep)
    w = Array{Float64,1}[]
    theta = Array{Theta,1}[]
    for i = 1:o.n_keep
        # sample k from p(k|t)
        ikeep = r.keepers[i]
        k[i] = rand_categorical(p_kt[:,r.t[ikeep]])

        # sample a from p(a|k,C) where C = clustering
        uz = sort(unique(r.z[:,i]))
        assign = zeros(Int,o.t_max+3)
        assign[uz] = randperm(k[i])[1:length(uz)]
        a[:,i] = assign[r.z[:,i]]

        # sample w from p(w|k,a)
        push!(w, rand_dirichlet(o.gamma .+ counts(a[:,i],k[i])))

        # sample from p(theta|k,a,phi)
        theta_i = [module_.new_theta(H)::Theta for j=1:k[i]]
        for j = 1:k[i]
            c_j = findfirst(isequal(j),assign)
            if c_j==nothing
                module_.prior_sample!(theta_i[j],H)
            else
                theta_i[j] = deepcopy(r.theta[c_j,i])
            end
        end
        push!(theta,theta_i)
    end
    return k,a,w,theta
end

# Compute log-likelihood of the data for a sequence of mixture weights w and component parameters theta.
function log_likelihoods(x,w,theta,mode)
    module_ = getfield(BayesianMixtures,Symbol(mode))
    log_f = module_.log_likelihood
    m = length(w)
    n = length(x)
    loglik = zeros(m)
    for iter = 1:m
        k = length(w[iter])
        ll = 0.0
        for i = 1:n
            ll += logsumexp(Float64[log(w[iter][j]) + log_f(x[i],theta[iter][j]) for j = 1:k])
        end
        loglik[iter] = ll
    end
    return loglik
end

# Compute the posterior similarity matrix (probability that i and j are in same cluster).
function similarity_matrix(result)
    o = result.options
    n = o.n
    z = result.z
    C = zeros(Int,n,n)
    n_used = 0
    for index = 1:o.n_keep
        if result.keepers[index]>o.n_burn
            for i = 1:n, j = 1:n
                if z[i,index]==z[j,index]; C[i,j] += 1; end
            end
            n_used += 1
        end
    end
    return C/n_used
end

# Compute cross-correlation between two vectors
function xcorr_(x,y)
    m,n = length(x),length(y)
    z = zeros(m+n-1)
    for k = 1:(m+n-1)
        s = 0.0
        for i = max(1,k-n+1):min(k,m)
            s += x[i]*y[n-k+i]
        end
        z[k] = s
    end
    return z
end

# Compute autocorrelation curve for the sequence y=[y_1,...,y_T]
function autocorrelation(y)
    T = length(y)
    u = ones(typeof(y[1]),T)
    d = [1:T; T-1:-1:1]
    autocov = xcorr_(y,y)./d - (xcorr_(y,u)./d).*(xcorr_(u,y)./d)
    a = autocov[T:-1:1]
    return a/a[1]
end
#autocorrelation_time(auto,c) = (tau=1+2*cumsum(auto[2:end]); M=findfirst([1:length(tau)] .>= c*tau); tau[M])

# ===================================================================
# ===================================================================
# ==================== Functions to plot results ====================
# ===================================================================
# ===================================================================

using Pkg
can_plot = ("PyPlot" in packages_installed)
if can_plot; using PyPlot; end
checkplotting() = (if !can_plot; error("Plotting is disabled since PyPlot is not installed."); end)

# ---------- General plotting functions ----------

function draw_now(number=0)
    checkplotting()
    if number>0; figure(number); end
    pause(0.001) # this should force the figure to update
    show() # bring figure window to the front
end

function open_figure(number;clear_figure=true,figure_size=(5,2.5))
    checkplotting()
    fig = figure(number, figsize=figure_size)
    subplots_adjust(top=0.85, bottom=0.2, left=0.15)
    clear_figure ? clf() : nothing
    # get_current_fig_manager()[:window][:showMaximized]() # maximize figure window
    return fig
end
function labels(title_string,xlabel_string="",ylabel_string="")
    checkplotting()
    title(title_string,fontsize=14)
    xlabel(xlabel_string,fontsize=14)
    ylabel(ylabel_string,fontsize=14)
end
function density_grid(x0,x1,y0,y1,xres,yres)
    xs = range(x0,stop=x1,length=xres)
    ys = range(y0,stop=y1,length=yres)
    X = [x for x=xs, y=ys]
    Y = [y for x=xs, y=ys]
    return xs,ys,X,Y
end

# ---------- Plotting results ----------
color_list = "cgyrbmkw"^1000
marker_list = "ov^Ds<>p*+x1234hHd"

function traceplot(values)
    checkplotting()
    n_values = length(values)
    n_subset = min(2000,n_values)
    if (n_subset < n_values); @warn("Traceplot shows only a subset of points."); end
    subset = round.(Int,range(1,stop=n_values,length=n_subset))
    jitter = (rand(n_subset).-0.5)/2
    PyPlot.plot(subset,values[subset]+jitter, "k.", markersize=1.0)
    draw_now()
end

function traceplot_timewise(result,t_show)
    checkplotting()
    timestep = result.elapsed_time/result.options.n_total
    t_show = min(result.elapsed_time,t_show)
    n_show = round(Int,t_show/timestep)
    n_mark = min(n_show,10000)
    subset = round.(Int,range(round(Int,n_show/n_mark),stop=n_show,length=n_mark))
    PyPlot.plot(subset*timestep,result.t[subset].+rand(n_mark)*0.5 .- 0.25,"k.",markersize=1)
    xlim(0,t_show)
    yticks(1:maximum(result.t[subset].+1))
    ylim(0,ylim()[2])
    labels("Traceplot of t","time (seconds)","t")
    draw_now()
end

function plot_t_running(result)
    checkplotting()
    o = result.options
    cdfs = t_running(result)
    tcolors = "brgycmk"^1000
    n_subset = min(2000,o.n_total)
    subset = round.(Int,range(1,stop=o.n_total,length=n_subset))
    for i = 1:size(cdfs,2)
        PyPlot.plot(subset,cdfs[subset,i],"$(tcolors[i])-")
    end
    yticks(0:0.2:1)
    labels("Running average of CDF(t|x)", "sweep number", "CDF(t|x)")
    draw_now()
end

function plot_autocorrelation(values,t_show,timestep;kwargs...)
    checkplotting()
    a = autocorrelation(values)
    n_plot = min(round(Int,t_show/timestep+1),length(a))
    PyPlot.plot(1000*(0:n_plot-1)*timestep,a[1:n_plot],linewidth=1.5;kwargs...)
    labels("Autocorrelation","time (milliseconds)","")
    xlim(0,1000*t_show)
    ylim(0,1)
    draw_now()
end

function plot_autocorrelation(result,t_show;kwargs...)
    checkplotting()
    o = result.options
    timestep = result.elapsed_time/o.n_total
    plot_autocorrelation(result.t[o.n_burn+1:o.n_total],t_show,timestep;kwargs...)
    labels("Autocorrelation of t","time (milliseconds)","")
    draw_now()
end

function plot_t_posterior(result; kwargs...)
    checkplotting()
    pt = t_posterior(result)
    PyPlot.plot(1:length(pt), pt; kwargs...)
    labels("Posterior on t","t (# of clusters)","p(t|data)")
    xticks(1:length(pt))
    draw_now()
end

function plot_t_posterior_average(results; kwargs...)
    checkplotting()
    t_posteriors = Array{Array{Float64,1},1}()
    for result in results
        push!(t_posteriors, t_posterior(result))
    end
    T = minimum(map(length,t_posteriors))
    avg = mean([p[1:T] for p in t_posteriors])[:]
    PyPlot.plot(1:T,avg; kwargs...)
    xticks(1:T)
    labels("Average posterior on t","t (# of clusters)","p(t|data)")
    draw_now()
end

function plot_k_posterior(result; kwargs...)
    checkplotting()
    pk = k_posterior(result)
    PyPlot.plot(1:length(pk), pk; kwargs...)
    labels("Posterior on k","k (# of components)","p(k|data)")
    xticks(1:length(pk))
    draw_now()
end

function plot_k_posterior_average(results; kwargs...)
    checkplotting()
    k_posteriors = Array{Array{Float64,1},1}()
    for result in results
        push!(k_posteriors, k_posterior(result))
    end
    T = minimum(map(length,k_posteriors))
    avg = mean([p[1:T] for p in k_posteriors])[:]
    PyPlot.plot(1:T,avg; kwargs...)
    xticks(1:T)
    labels("Average posterior on k","k (# of components)","p(k|data)")
    draw_now()
end

# Posterior similarity matrix (probability that i and j are in same cluster)
function plot_similarity_matrix(result; step=10)
    checkplotting()
    n = result.options.n
    title("Posterior similarity matrix")
    imshow(1.0.-similarity_matrix(result),cmap="gray",interpolation="nearest",vmin=0,vmax=1,extent=(1,n,n,1),origin="upper")
    tick_params(axis="x",direction="out",top="off")
    tick_params(axis="y",direction="out",top="off")
    xticks(step:step:n)
    yticks(step:step:n)
    draw_now()
end

function plot_clusters(x,z; colors=color_list, markers=marker_list, markersize=2)
    checkplotting()
    n = length(x)
    d = length(x[1]) # dimension of data
    N = zeros(n); for i = 1:n; N[z[i]] += 1; end
    zcs = sortperm(N,rev=true)
    t = sum(N .> 0)
    if d==1 # univariate data
        for c = 1:t
            xc = x[z.==zcs[c]]
            jitter = (rand(length(xc)).-0.5)/4
            PyPlot.plot(xc,c.+jitter,"bx",markersize=2.0)
        end
        yticks(1:t)
        ylim(0.5,t+0.5)
        grid(true,axis="y")
    elseif d==2 # bivariate data
        T = min(length(colors),length(markers))
        if t>T; @warn("Only plotting the largest $T clusters for $label."); end
        for c = 1:min(t,T)
            xc = x[z.==zcs[c]]
            PyPlot.plot([xi[1] for xi in xc],[xi[2] for xi in xc],
                "$(colors[c])$(markers[c])",markersize=markersize,markeredgewidth=0.1)
        end
    else
        error("Plotting clusters is only enabled for univariate or bivariate data.")
    end
    title("Clusters")
    draw_now()
end

function rug_plot(data) # Rug plot of data
    checkplotting()
    @assert(all(map(length,data).==1),"Rug plot only enabled for univariate data.")
    PyPlot.plot(data,zeros(length(data)),"k+",markersize = 18)
    tick_params(axis="x",direction="out",top="off")
    title("Rug plot")
    draw_now()
end

function plot_density_estimate(result; resolution=-1, kwargs...)
    checkplotting()
    bounds(y) = (mn=minimum(y); mx=maximum(y); s=(mx-mn)/8; (round(mn-s,digits=0),round(mx+s,digits=0)))
    o = result.options
    data,n = o.x,o.n
    d = length(data[1]) # dimension of data
    if d==1
        # Density estimate
        if resolution==-1; resolution = 500; end
        xmin,xmax = bounds(data)
        xs = range(xmin,stop=xmax,length=resolution)
        fs = Float64[density_estimate(xi,result) for xi in xs]
        PyPlot.plot(xs,fs; kwargs...)
    elseif d==2
        if resolution==-1; resolution = 40; end
        n_contours = 10
        power = 0.25
        data_x = Float64[data[i][1] for i=1:n]
        data_y = Float64[data[i][2] for i=1:n]
        xmin,xmax = bounds(data_x)
        ymin,ymax = bounds(data_y)
        xs,ys,X,Y = density_grid(xmin,xmax,ymin,ymax,resolution,resolution)
        subplots_adjust(top=0.845)
        Z = Float64[density_estimate(Float64[x;y],result) for x=xs, y=ys]
        PyPlot.contour(X,Y,Z.^power,n_contours)
    else
        error("Plotting density estimates is only enabled for dimensions 1 and 2.")
    end
    title("Density estimate")
    draw_now()
end

end # module BayesianMixtures

ext: jl


text:
"""
    angularspectrum(Uin, d, λ, lx, ly)

calculate the propagation light field based on the angular spectrum.

## Arguments

- `Uin::AbstractArray{<:Number,2}`: Complex array of input complex amplitude.
- `d::Real`: Distance to propagate in metres.
- `λ::Real`: Wavelength of light to propagate.
- `lx::Real`: The physical size of the input data along the x-axis.
- `ly::Real`: The physical size of the input data along the y-axis.

## Returns

- `::Array{<:Number,2}`: Complex amplitude data after propagation.
"""
function angularspectrum(Uin::AbstractArray{<:Number,2}, d::Real, λ::Real, lx::Real, ly::Real)
    (n, m) = size(Uin)
    ua = (-n/2:n/2-1)/lx*λ
    va = (-m/2:m/2-1)/ly*λ
    ifft(fft(Uin).*ifftshift([1.0-u^2-v^2>=0.0 ? exp(2im*pi*d/λ*sqrt(1.0-u^2-v^2)) : 0.0+0.0im for u in ua, v in va]))
end

"""
    angularspectrum(Uin::MonoLightField2D, d::Unitful.Length) -> MonoLightField2D

calculate the propagation light field for [`MonoLightField2D`](@ref) based on the angular spectrum.
"""
function angularspectrum(Uin::MonoLightField2D, d::Unitful.Length)
    Uout = angularspectrum(Uin.data, auval(d), auval(Uin.wavelength), auval(Uin.size[1]), auval(Uin.size[2]))
    MonoLightField2D(Uin, data=Uout)
end

"""
    angularspectrum(d)

represent an [`angularspectrum`](@ref) propagation functor with the specified distance.
"""
angularspectrum(d) = x->angularspectrum(x, d)


text:
include("SimulationController.jl")
include("TxDAQController.jl")


text:
import accrue.cryptoerase.runtime.Condition;

class C {
    static final Condition{L} c = new Condition();

    int xTESTOUTPUT = 7;
    int yTESTOUTPUT = 7;

    void m() {
	int pubK = (PUBKEY(H){L /this.c H}) 42;
	int privK = (PRIVKEY(H){L /this.c H}) 42;

	this.xTESTOUTPUT = pubK;
	this.yTESTOUTPUT = privK;

    }

    public static void main(String[] args) {
	new C().m();
    }
}

max_stars_repo_name: HarvardPL/cryptoerase

max_issues_repo_name: HarvardPL/cryptoerase

max_forks_repo_name: HarvardPL/cryptoerase


text:
@testset "Momentum" begin
    Random.seed!(SEED)
    @testset "Array" begin
        x = cumsum(randn(N))     # close
        Y = cumsum(randn(N, 2), dims=1)  # high-low
        Z = cumsum(randn(N, 3), dims=1)  # high-low-close
        tmp = aroon(Y)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 3
        @test sum(isnan.(tmp)) != N
        tmp = donch(Y)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 3
        @test sum(isnan.(tmp)) != N
        tmp = momentum(x)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        @test sum(isnan.(tmp)) != N
        tmp = roc(x)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        @test sum(isnan.(tmp)) != N
        tmp = macd(x)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 3
        @test sum(isnan.(tmp)) != N
        tmp = rsi(x)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        @test sum(isnan.(tmp)) != N
        tmp = adx(Z)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 3
        @test sum(isnan.(tmp)) != N
        tmp = adx(Z, wilder=true)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 3
        @test sum(isnan.(tmp)) != N
        tmp = psar(Y)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        @test sum(isnan.(tmp)) != N
        tmp = kst(x)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        @test sum(isnan.(tmp)) != N
        tmp = wpr(Z)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        @test sum(isnan.(tmp)) != N
        tmp = cci(Z)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        @test sum(isnan.(tmp)) != N
        tmp = stoch(Z, kind=:fast)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 2
        @test sum(isnan.(tmp)) != N
        tmp = stoch(Z, kind=:slow)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 2
        @test sum(isnan.(tmp)) != N
        tmp = smi(Z)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 2
        @test sum(isnan.(tmp)) != N
    end
    @testset "Temporal" begin
        x = TS(cumsum(randn(N)))
        Y = TS(cumsum(randn(N, 2), dims=1))
        Z = TS(cumsum(randn(N, 3), dims=1))
        # momentum function
        tmp = aroon(Y)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 3
        tmp = donch(Y)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 3
        tmp = momentum(x)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        tmp = roc(x)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        tmp = macd(x)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 3
        tmp = rsi(x)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        tmp = adx(Z)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 3
        tmp = adx(Z, wilder=true)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 3
        tmp = psar(Y)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        tmp = kst(x)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        tmp = wpr(Z)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        tmp = cci(Z)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        tmp = stoch(Z, kind=:fast)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 2
        tmp = stoch(Z, kind=:slow)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 2
        tmp = smi(Z)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 2
    end
end


text:
module MCDiagnostic
using LinearAlgebra
using DSP: conv
using StatsBase: autocor
using Statistics: var, mean, cov
using Distributions: FDist, quantile

 #=================================================================================
                                IS/SMC Diagnostics
 ==================================================================================#

is_eff_ss(W) = 1/sum((W./sum(W)).^2)

 #=================================================================================
                                MCMC Diagnostics
 ==================================================================================#

# not the most efficient AR calculation, but will do until optimisation reqd
ar_coeffs(X, max_l) = reduce(hcat, [X[i:end-max_l+i-1] for i in 1:max_l]) \ X[max_l+1:end]

function spectrum0_ar(x::Array{T,1}; ncoeffs=nothing) where T <: Real
    nrows = length(x)
    ncoeffs = something(ncoeffs, Int64(ceil(10 * log10(nrows)))) # for AR (not doing AIC)

    # if column is (almost) perfect line, spectrum is 0
    # (avoids pathologies in AR calculation)
    if isapprox(std(diff(x)), 0.)
        spec = 0.
        ncoeffs = 0
    else
        cx = x .- mean(x)
        coeffs = ar_coeffs(cx, ncoeffs)
        x̂ = vcat(cx[1], conv(cx, coeffs)[1:end-ncoeffs])
        res_var = var(cx) - var(x̂)
        spec = res_var/(1 - sum(coeffs))^2
    end

    return convert(T, spec), ncoeffs
end

#= In order for this to function more effectively, the spectrum0_ar
   function needs to choose the order of AR process via AIC. While not
   especially difficult, it is also not trivial to add this in, esp. in
   an efficient manner, so I have left this as WIP for now.
=#
function mcmc_effective_size_ar(x::Array{T,1}) where T <: Real
    @warn "Currently WIP. Need to implement AIC for AR model selection."
    spec, _order = spectrum0_ar(x)
    m = length(x)
    ess = spec != 0. ? m * var(x) / spec : 0.
    return ess
end

function mcmc_effective_size_asymp(x::Array{T,1}) where T <: Real
    # see e.g. BDA3 (sec 11.5) or Thompson 2010: A comparison of ... autocorrelation time.
    n = length(x)
    maxlag = min(n-5, Int64(ceil(30*log10(length(x)))))  # heuristic max
    ρ = autocor(x, 1:maxlag)

    # Geyer, 1992: sums of consecutive ACF values are > 0.
    consec_pairs = conv(ρ, ones(2))
    ncoeffs = something(findfirst(consec_pairs .< 0), length(ρ)+1) - 1
    return  n / (1 + 2*(sum(ρ[1:ncoeffs])))
end

# coda's version of Rhat is really nice. Converted approx. as-is.
function Rhat(xs...; confidence = 0.95)
    #=================================================================
    Unlike BDA3, Plummer et al. in R's coda package use an additional
    adjustment for 1/nchains, as well as add some CIs.
    =================================================================#
    @assert (length(unique(map(length, xs))) == 1) "inconsistent lengths in xs"
    x = reduce(hcat, xs)
    n, nchains = size(x)

    #=================================================================
    Estimate mean within-chain variance (W) and between-chain variance
    (B/Niter), and calculate sampling variances and covariance of the
    estimates (varW, varB, covWB)
    =================================================================#
    ## Univariate (lower case)
    s2 = var(x, dims=1)
    w  = mean(s2)
    x̄  = mean(x, dims=1)
    b  = n * var(x̄)

    muhat  = mean(x̄)
    var_w  = var(s2)/nchains
    var_b  = (2 * b^2)/(nchains - 1)
    cov_wb = (n/nchains) * cov(s2[:], x̄[:].^2) - 2 * muhat * cov(s2[:], x̄[:])

    #=================================================================
    Posterior interval combines all uncertainties in a t interval with
    center muhat, scale sqrt(V), and df_V degrees of freedom.
    =================================================================#
    V     = (n-1)w / n  + (1 + 1/nchains)b / n
    var_V = (var_w * (n-1)^2 +
           (1 + 1/nchains)^2 * var_b + 2 * (n-1) * (1 + 1/nchains) * cov_wb
          )/n^2
    df_V = 2V^2 / var_V

    #=================================================================
                      R = sqrt(V/W) * df_adj
    where df_adj is a degrees of freedom adjustment for the width
    of the t-interval.

    To calculate upper confidence interval we divide R2 = R^2 into two
    parts, fixed and random.  The upper limit of the random part is
    calculated assuming that B/W has an F distribution.
    =================================================================#

    df_adj      = (df_V + 3)/(df_V + 1)
    B_df        = nchains - 1
    W_df        = (2 * w^2)/var_w
    R2_fixed    = (n - 1)/n
    R2_random   = (1 + 1/nchains) * (1/n) * (b/w)
    R2_estimate = R2_fixed + R2_random
    R2_upper    = R2_fixed + quantile(FDist(B_df, W_df), (1 + confidence)/2) * R2_random
    R̂           = sqrt(df_adj * R2_estimate)
    R̂_ci_upper  = sqrt(df_adj * R2_upper)

    return R̂, R̂_ci_upper
end

end


text:
# run all ETL scripts for the disaster database
using Dates

scripts = [
    "./src/load_s_state_geocodes.jl",
    "./src/load_us_county_shp.jl",
    "./src/load_s_earthquake.jl",
    "./src/load_s_storm.jl",
    "./src/load_s_wildfire.jl",
    "./src/make_d_datestamp.jl",
    "./src/make_d_location.jl",
    "./src/make_d_severity.jl",
    "./src/make_f_earthquake.jl",
    "./src/make_f_storm.jl",
    "./src/make_f_wildfire.jl"
];

load_starttime = Dates.now()
elapsed_times = [];
for script in scripts
    starttime = Dates.now()
    run(`julia $(script)`)
    endtime = Dates.now()
    push!(elapsed_times, "$(script): elapsed time is $(Dates.canonicalize(Dates.CompoundPeriod(endtime - starttime)))")
end

println("\n\nload...")
for elapsed_time in elapsed_times
    println(elapsed_time)
end

load_endtime = Dates.now()
println("load: total elapsed time is $(Dates.canonicalize(Dates.CompoundPeriod(load_endtime - load_starttime)))")


text:
"""
    A = get_nrm(pedlist [, check=false, sorthere=true])
    A = get_nrm(T, pedlist [, check=false, sorthere=true])

Calculates the numerator (or additive) relationship matrix, so-called *A-matrix* using the recursive ("tabular") method from the pedigree list `pedlist`.
The pedigree must be sorted as the ancestors precede their progeny unless putting `sorthere=true`.
The resulting matrix is dense. You can specify the type of `A` as `T` (Float64 by default).
By default, this function checks if the pedigree is sorted chronologically or not; use `check=false` to turn it off.
With the option `sorthere=true`, the function reorders the pedigree internally, and returns the correct A-matrix with the original order.
However, it creates a temporary A-matrix with the reordered index, and it needs twice more memory.
"""
function get_nrm(Tv::DataType, pedlist::Matrix{Ti}; check=true, sorthere=false) where Ti<:Integer
   # pedigree checks
   if check
      if sorthere
         # being sorted so no checks for `parentsfirst`
         if !check_ped(ped)
            error("error in pedigree")
         end
      else
         # assuming to be already sorted
         if !check_ped(ped,parentsfirst=true)
            error("error in pedigree, or being not sorted")
         end
      end
   end
   # tabular method
   if sorthere
      # copy the pedigree list
      ped = normalize_pedigree(pedlist)
      perm,invp = find_ped_order(ped)
      permute_ped!(invp,ped)
      B = tabular_method_dense(Tv, ped)
      return B[invp,invp]
   else
      return tabular_method_dense(Tv, pedlist)
   end
end
function get_nrm(pedlist::Matrix{Ti}, Tv=Float64::DataType; check=false, sorthere=false) where Ti<:Integer
   return get_nrm(Tv, pedlist, check=check, sorthere=sorthere)
end

# lower triangular only
function tabular_method_dense(Tv::DataType, pedlist::Matrix{Ti}) where Ti<:Integer
   n::Ti = size(pedlist,2)
   sid::Ti = 0
   did::Ti = 0
   A = zeros(Tv,n,n)
   for j=1:n
      # IDs
      sid = pedlist[1,j]
      did = pedlist[2,j]
      # diagonal
      A[j,j] = 1 + aij(A,n,sid,did)/2
      # off diagonal (i>j i.e. sid<j & did<j)
      for i=j+1:n
         sid = pedlist[1,i]
         did = pedlist[2,i]
         A[i,j] = (aij(A,n,j,sid)+aij(A,n,j,did))/2
         A[j,i] = A[i,j]
      end
   end
   return A
end

# refers to A(i,j)
function aij(A::Matrix{Tv}, n::Ti, i::Ti, j::Ti) where {Tv<:AbstractFloat,Ti<:Integer}
   if i<1 || i>n
      val = zero(Tv)
   elseif j<1 || j>n
      val = zero(Tv)
   else
      val = A[i,j]
   end
   return val
end

"""
    Ainv = get_nrminv(pedlist, f; rank=0)
    Ainv = get_nrminv(pedlist; f, rank=0)

Computes the inverse of numerator (additive) relationship matrix using the "Henderson's" method from the pedigree list `pedlist` and inbreeding coefficients `f`.
Without `f`, this function assumes zero inbreeding for all individuals.
See Henderson (1976) and Quaas (1976) for details.
The output is SparseMatrixCSC with the same type of elements with `f` while SparseMatrixDict (i.e., a hash matrix) is used as a temporary matrix.
If you want to fix the size of A-inverse, specify `rank`.
The function creates the inverse with either a bigger size: `rank` or the maximum ID in the pedigree.
"""
function get_nrminv(pedlist::Matrix{Ti}, f::Vector{Tv}; rank=0) where {Ti<:Integer,Tv<:AbstractFloat}
   # there are UPGs if n<m
   n::Ti = size(pedlist,2)
   m::Ti = max(maximum(pedlist),n,rank)
   AInv = SparseMatrixDict{Tv,Int}(m,m)
   w::Array{Tv,1} = [1.0, -0.5, -0.5]
   triplet = zeros(Ti,3)
   for k=1:n
      sid = pedlist[1,k]
      did = pedlist[2,k]
      triplet .= [k, sid, did]
      if sid<1 || sid>n
         fs = 0.0
         ms = 1.0
      else
         fs = f[sid]
         ms = 0.0
      end
      if did<1 || did>n
         fd = 0.0
         md = 1.0
      else
         fd = f[did]
         md = 0.0
      end
      b = 4.0/((1.0-fs)*(1.0+ms)+(1.0-fd)*(1.0+md))
      for j=1:3
         for i=1:3
            aj = triplet[j]
            ai = triplet[i]
            #if (ai>=1 && ai<=n) && (aj>=1 && aj<=n)
            if ai>=1 && aj>=1
               val = b*w[i]*w[j]
               AInv[ai,aj] = AInv[ai,aj] + val
            end
         end
      end
   end
   return sparse(AInv)
end
function get_nrminv(pedlist::Matrix{Ti}; f::Vector{Tv}=Float64[], rank=0) where {Ti<:Integer,Tv<:AbstractFloat}
   if length(f)==0
      inb = zeros(Float64,maximum(pedlist))
   else
      inb = f
   end
   return get_nrminv(pedlist,inb, rank=rank)
end

"""
    Ainv = get_mgsnrminv(mgslist; rank=0)

Computes the inverse of numerator (additive) relationship matrix using the "Henderson's" method from the pedigree list `mgslist` with sire and maternal-grandsire (MGS).
The inbreeding will be ignored.
See Henderson (1975) for details.
The output is SparseMatrixCSC with Float64 while SparseMatrixDict (i.e., a hash matrix) is used as a temporary matrix.
If you want to fix the size of A-inverse, specify `rank`.
The function creates the inverse with either a bigger size: `rank` or the maximum ID in the pedigree.
"""
function get_mgsnrminv(mgslist::Matrix{Ti}, males=fill(true,size(mgslist,2)); rank=0) where {Ti<:Integer}
   # there are UPGs if n<m
   n::Ti = size(mgslist,2)
   m::Ti = max(maximum(mgslist),n,rank)
   AInv = SparseMatrixDict{Float64,Int}(m,m)
   w::Array{Float64,1} = [1.0, -0.5, -0.25]
   triplet = zeros(Ti,3)
   for k=1:n
      sid = mgslist[1,k]
      mgsid = mgslist[2,k]
      triplet .= [k, sid, mgsid]
      if (sid>=1 && sid<=n) && (mgsid>=1 && mgsid<=n)
         m = 0.0
      elseif (sid>=1 && sid<=n)
         m = 1.0
      elseif (mgsid>=1 && mgsid<=n)
         m = 4.0
      else
         m = 5.0
      end
      b = 16.0/(m + 11.0)
      if males[k]
         for j=1:3
            for i=1:3
               aj = triplet[j]
               ai = triplet[i]
               if ai>=1 && aj>=1
                  val = b*w[i]*w[j]
                  AInv[ai,aj] = AInv[ai,aj] + val
               end
            end
         end
      end
   end
   return sparse(AInv)
end

# change UPG code to 0
function normalize_id_unchanged(pid,maxid)
   return pid
end
function normalize_id_removing_upg(pid,maxid)
   return ifelse(pid>maxid,0,pid)
end

"""
    Q = get_qm(pedlist)
    Q = get_qm(T, pedlist)

Computes the Q-matrix relating unknown parent groups (UPG) to individuals, as suggested by Quaas (1988).
The pedigree list `pedlist` should include a UPG code that is greater than the maximum ID of real animals.
The type of `Q` is `T` (Float64 by default).
"""
function get_qm(Tv::DataType, pedlist::Matrix{Ti}) where Ti<:Integer
   # real animals
   n::Ti = size(pedlist,2)
   # UPGs
   nupg::Ti = maximum(pedlist) - n
   # no UPGs
   if nupg < 1
      error("No groups in the pedigree list")
   end
   # start computations
   Q = zeros(Tv,n,nupg)
   Qp = zeros(Tv,nupg)
   for i=1:n
      Qp .= zeros(Tv,nupg)
      calculate_upg_q!(i,pedlist,n,1.0,Qp)
      Q[i,:] = Qp
   end
   return Q
end
function get_qm(pedlist::Matrix{Ti},Tv::DataType=Float64) where Ti<:Integer
   return get_qm(Tv,pedlist)
end

function calculate_upg_q!(aid::Ti, pedlist::Matrix{Ti}, maxanim::Ti, g::Tv, Qp::Vector{Tv}) where {Ti<:Integer,Tv<:AbstractFloat}
   if aid > maxanim
      group = aid - maxanim
      Qp[group] = Qp[group] + g
   elseif aid > 0
      calculate_upg_q!(pedlist[1,aid],pedlist,maxanim,g/2,Qp)
      calculate_upg_q!(pedlist[2,aid],pedlist,maxanim,g/2,Qp)
   end
end

"""
    P = get_pm(pedlist)
    P = get_pm(T, pedlist)

Computes a P-matrix that links parents with their progeny, as shown in Quaas (1988).
The type of `P` is `T` (Float64 by default).
"""
function get_pm(Tv::DataType, pedlist::Matrix{Ti}; verbose::Bool=true) where Ti<:Integer
   # real animals
   n = size(pedlist,2)
   # UPGs
   nupg = maximum(pedlist) - n
   # zero index
   if minimum(pedlist)<1
      if verbose
         @warn "pedigree with 0 index; not all animals assined to UPG"
      end
      P = zeros(Tv,n,n)
   else
      P = zeros(Tv,n,n+nupg)
   end
   # start computations
   for i=1:n
      sid = pedlist[1,i]
      did = pedlist[2,i]
      if sid>0
         P[i,sid] = 0.5
      end
      if did>0
         P[i,did] = 0.5
      end
   end
   return P
end
function get_pm(pedlist::Matrix{Ti}, Tv::DataType=Float64; verbose::Bool=true) where Ti<:Integer
   return get_pm(Tv,pedlist,verbose=verbose)
end


text:
@theme VisualStudioTheme Dict(
    :name => "Visual Studio",
    :description => "Based on the theme from MS Visual Studio.",
    :comments => "Based on the Visual Studio theme from Pygments.",
    :style => S"bg: ffffff",
    :tokens => Dict(
        COMMENT            => S"fg: 008000",
        COMMENT_PREPROC    => S"fg: 0000ff",
        KEYWORD            => S"fg: 0000ff",
        OPERATOR_WORD      => S"fg: 0000ff",
        KEYWORD_TYPE       => S"fg: 2b91af",
        NAME_CLASS         => S"fg: 2b91af",
        STRING             => S"fg: a31515",
        GENERIC_HEADING    => S"bold",
        GENERIC_SUBHEADING => S"bold",
        GENERIC_EMPH       => S"italic",
        GENERIC_STRONG     => S"bold",
        GENERIC_PROMPT     => S"bold",
        ERROR              => S"bg: ff0000; fg: ffffff",
        TEXT               => S"",
    ),
)


text:
## Analyzing lowered code

function add_docexpr!(docexprs::AbstractDict{Module,V}, mod::Module, ex) where V
    docexs = get(docexprs, mod, nothing)
    if docexs === nothing
        docexs = docexprs[mod] = V()
    end
    push!(docexs, ex)
    return docexprs
end

function lookup_callexpr(frame, stmt)
    fargs = JuliaInterpreter.collect_args(frame, stmt)
    return Expr(:call, fargs...)
end

function assign_this!(frame, value)
    frame.framedata.ssavalues[frame.pc] = value
end

# This defines the API needed to store signatures using methods_by_execution!
# This default version is simple; there's a more involved one in the file Revise.jl
# that interacts with CodeTracking.
const MethodInfo = IdDict{Type,LineNumberNode}
add_signature!(methodinfo::MethodInfo, @nospecialize(sig), ln) = push!(methodinfo, sig=>ln)
push_expr!(methodinfo::MethodInfo, mod::Module, ex::Expr) = methodinfo
pop_expr!(methodinfo::MethodInfo) = methodinfo
add_dependencies!(methodinfo::MethodInfo, be::BackEdges, src, chunks) = methodinfo
add_includes!(methodinfo::MethodInfo, mod::Module, filename) = methodinfo

function minimal_evaluation!(methodinfo, frame)
    src = frame.framecode.src
    be = BackEdges(src)
    chunks = toplevel_chunks(be)
    musteval = falses(length(src.code))
    for chunk in chunks
        if hastrackedexpr(frame.framecode.src, chunk)
            musteval[chunk] .= true
        end
    end
    # Conservatively, we need to step in to each Core.eval in case the expression defines a method.
    hadeval = false
    for id in eachindex(src.code)
        stmt = src.code[id]
        me = false
        if isa(stmt, Expr)
            if stmt.head === :call
                f = stmt.args[1]
                me |= f === :include
                me |= JuliaInterpreter.hasarg(isequal(:eval), stmt.args)
            end
        end
        if me
            chunkid = findfirst(chunk->id∈chunk, chunks)
            musteval[chunks[chunkid]] .= true
        end
    end
    add_dependencies!(methodinfo, be, src, chunks)
    return musteval
end

function methods_by_execution(mod::Module, ex::Expr; kwargs...)
    methodinfo = MethodInfo()
    docexprs = DocExprs()
    value, frame = methods_by_execution!(finish_and_return!, methodinfo, docexprs, mod, ex; kwargs...)
    return methodinfo, docexprs, frame
end

function methods_by_execution!(@nospecialize(recurse), methodinfo, docexprs, mod::Module, ex::Expr; always_rethrow=false, define=true, kwargs...)
    lwr = Meta.lower(mod, ex)
    isa(lwr, Expr) || return nothing, nothing
    frame = prepare_thunk(mod, copy(lwr), true)
    frame === nothing && return nothing, nothing
    define || LoweredCodeUtils.rename_framemethods!(recurse, frame)
    # Determine whether we need interpreted mode
    musteval = minimal_evaluation!(methodinfo, frame)
    if !any(musteval)
        # We can evaluate the entire expression in compiled mode
        if define
            ret = try
                Core.eval(mod, ex) # evaluate in compiled mode if we don't need to interpret
            catch err
                (always_rethrow || isa(err, InterruptException)) && rethrow(err)
                loc = location_string(whereis(frame)...)
                @error "(compiled mode) evaluation error starting at $loc" mod ex exception=(err, trim_toplevel!(catch_backtrace()))
                nothing
            end
        else
            ret = nothing
        end
    else
        # Use the interpreter
        # We have to turn off all active breakpoints, https://github.com/timholy/CodeTracking.jl/issues/27
        bp_refs = JuliaInterpreter.breakpoints()
        if eltype(bp_refs) !== JuliaInterpreter.BreakpointRef
            bp_refs = JuliaInterpreter.BreakpointRef[]
            foreach(bp -> append!(bp_refs, bp.instances), bp_refs)
        end
        active_bp_refs = filter(bp->bp[].isactive, bp_refs)
        foreach(disable, active_bp_refs)
        ret = try
            methods_by_execution!(recurse, methodinfo, docexprs, frame, musteval; define=define, kwargs...)
        catch err
            (always_rethrow || isa(err, InterruptException)) && rethrow(err)
            loc = location_string(whereis(frame)...)
            @error "evaluation error starting at $loc" mod ex exception=(err, trim_toplevel!(catch_backtrace()))
            nothing
        end
        foreach(enable, active_bp_refs)
    end
    return ret, lwr
end

function methods_by_execution!(@nospecialize(recurse), methodinfo, docexprs, frame, musteval; define::Bool=true, skip_include::Bool=true)
    isok(lnn::LineTypes) = !iszero(lnn.line) || lnn.file !== :none   # might fail either one, but accept anything

    mod = moduleof(frame)
    # Hoist this lookup for performance. Don't throw even when `mod` is a baremodule:
    modinclude = isdefined(mod, :include) ? getfield(mod, :include) : nothing
    signatures = []  # temporary for method signature storage
    pc = frame.pc
    while true
        JuliaInterpreter.is_leaf(frame) || (@warn("not a leaf"); break)
        if !musteval[pc] && !define
            pc = next_or_nothing!(frame)
            pc === nothing && break
            continue
        end
        stmt = pc_expr(frame, pc)
        if isa(stmt, Expr)
            head = stmt.head
            if head ∈ structheads
                if define
                    pc = step_expr!(recurse, frame, stmt, true)  # This should check that they are unchanged
                else
                    pc = next_or_nothing!(frame)
                end
            elseif head === :thunk && isanonymous_typedef(stmt.args[1])
                # Anonymous functions should just be defined anew, since there does not seem to be a practical
                # way to "find" them. They may be needed to define later signatures.
                # Note that named inner methods don't require special treatment
                pc = step_expr!(recurse, frame, stmt, true)
            elseif head === :method
                empty!(signatures)
                ret = methoddef!(recurse, signatures, frame, stmt, pc; define=define)
                if ret === nothing
                    # This was just `function foo end` or similar.
                    # However, it might have been followed by a thunk that defined a
                    # method (issue #435), so we still need to check for additions.
                    if !isempty(signatures)
                        file, line = whereis(frame.framecode, pc)
                        lnn = LineNumberNode(Int(line), Symbol(file))
                        for sig in signatures
                            add_signature!(methodinfo, sig, lnn)
                        end
                    end
                    pc = ret
                else
                    pc, pc3 = ret
                    # Get the line number from the body
                    stmt3 = pc_expr(frame, pc3)
                    lnn = nothing
                    if line_is_decl
                        sigcode = @lookup(frame, stmt3.args[2])::Core.SimpleVector
                        lnn = sigcode[end]
                        if !isa(lnn, LineNumberNode)
                            lnn = nothing
                        end
                    end
                    if lnn === nothing
                        bodycode = stmt3.args[end]
                        if !isa(bodycode, CodeInfo)
                            bodycode = @lookup(frame, bodycode)
                        end
                        if isa(bodycode, CodeInfo)
                            lnn = linetable(bodycode, 1)
                            if !isok(lnn)
                                lnn = nothing
                                if length(bodycode.code) > 1
                                    # This may be a kwarg method. Mimic LoweredCodeUtils.bodymethod,
                                    # except without having a method
                                    stmt = bodycode.code[end-1]
                                    if isa(stmt, Expr) && length(stmt.args) > 1
                                        stmt = stmt::Expr
                                        a = stmt.args[1]
                                        nargs = length(stmt.args)
                                        hasself = let stmt = stmt, slotnames::Vector{Symbol} = bodycode.slotnames
                                            any(i->LoweredCodeUtils.is_self_call(stmt, slotnames, i), 2:nargs)
                                        end
                                        if isa(a, Core.SlotNumber)
                                            a = bodycode.slotnames[a.id]
                                        end
                                        if hasself && (isa(a, Symbol) || isa(a, GlobalRef))
                                            thismod, thisname = isa(a, Symbol) ? (mod, a) : (a.mod, a.name)
                                            if isdefined(thismod, thisname)
                                                f = getfield(thismod, thisname)
                                                mths = methods(f)
                                                if length(mths) == 1
                                                    mth = first(mths)
                                                    lnn = LineNumberNode(Int(mth.line), mth.file)
                                                end
                                            end
                                        end
                                    end
                                end
                                if lnn === nothing
                                    # Just try to find *any* line number
                                    for lnntmp in linetable(bodycode)
                                        lnntmp = lnntmp::LineTypes
                                        if isok(lnntmp)
                                            lnn = lnntmp
                                            break
                                        end
                                    end
                                end
                            end
                        elseif isexpr(bodycode, :lambda)
                            bodycode = bodycode::Expr
                            lnntmp = bodycode.args[end][1]::LineTypes
                            if isok(lnntmp)
                                lnn = lnntmp
                            end
                        end
                    end
                    if lnn === nothing
                        i = codelocs(frame, pc3)
                        while i > 0
                            lnntmp = linetable(frame, i)
                            if isok(lnntmp)
                                lnn = lnntmp
                                break
                            end
                            i -= 1
                        end
                    end
                    if lnn !== nothing && isok(lnn)
                        for sig in signatures
                            add_signature!(methodinfo, sig, lnn)
                        end
                    end
                end
            elseif head === :(=) && isa(stmt.args[1], Symbol)
                if define
                    pc = step_expr!(recurse, frame, stmt, true)
                else
                    # FIXME: Code that initializes a global, performs some operations that
                    # depend on the value, and then mutates it will run into serious trouble here.
                    # sym = stmt.args[1]
                    # if isconst(mod, sym)
                    rhs = stmt.args[2]
                    val = isa(rhs, Expr) ? JuliaInterpreter.eval_rhs(recurse, frame, rhs) : @lookup(frame, rhs)
                    assign_this!(frame, val)
                    pc = next_or_nothing!(frame)
                    # else
                    #     pc = step_expr!(recurse, frame, stmt, true)
                    # end
                end
            elseif head === :call
                f = @lookup(frame, stmt.args[1])
                if f === Core.eval
                    # an @eval or eval block: this may contain method definitions, so intercept it.
                    evalmod = @lookup(frame, stmt.args[2])::Module
                    evalex = @lookup(frame, stmt.args[3])
                    thismodexs, thisdocexprs = split_expressions(evalmod, evalex; extract_docexprs=true)
                    for (m, docexs) in thisdocexprs
                        for docex in docexs
                            add_docexpr!(docexprs, m, docex)
                        end
                    end
                    value = nothing
                    for (newmod, newex) in thismodexs
                        newex = unwrap(newex)
                        newframe = prepare_thunk(newmod, newex)
                        newframe === nothing && continue
                        define || LoweredCodeUtils.rename_framemethods!(recurse, newframe)
                        newmusteval = minimal_evaluation!(methodinfo, newframe)
                        push_expr!(methodinfo, newmod, newex)
                        value = methods_by_execution!(recurse, methodinfo, docexprs, newframe, newmusteval; define=define)
                        pop_expr!(methodinfo)
                    end
                    assign_this!(frame, value)
                    pc = next_or_nothing!(frame)
                elseif skip_include && (f === modinclude || f === Base.include || f === Core.include)
                    # Skip include calls, otherwise we load new code
                    add_includes!(methodinfo, mod, @lookup(frame, stmt.args[2]))
                    assign_this!(frame, nothing)  # FIXME: the file might return something different from `nothing`
                    pc = next_or_nothing!(frame)
                elseif !define && f === Base.Docs.doc!
                    fargs = JuliaInterpreter.collect_args(frame, stmt)
                    popfirst!(fargs)
                    length(fargs) == 3 && push!(fargs, Union{})  # add the default sig
                    dmod::Module, b::Base.Docs.Binding, str::Base.Docs.DocStr, sig = fargs
                    m = get!(Base.Docs.meta(dmod), b, Base.Docs.MultiDoc())::Base.Docs.MultiDoc
                    if haskey(m.docs, sig)
                        currentstr = m.docs[sig]::Base.Docs.DocStr
                        redefine = currentstr.text != str.text
                    else
                        push!(m.order, sig)
                        redefine = true
                    end
                    # (Re)assign without the warning
                    if redefine
                        m.docs[sig] = str
                        str.data[:binding] = b
                        str.data[:typesig] = sig
                    end
                    assign_this!(frame, Base.Docs.doc(b, sig))
                    pc = next_or_nothing!(frame)
                else
                    # A :call Expr we don't want to intercept
                    pc = step_expr!(recurse, frame, stmt, true)
                end
            else
                # An Expr we don't want to intercept
                pc = step_expr!(recurse, frame, stmt, true)
            end
        else
            # A statement we don't want to intercept
            pc = step_expr!(recurse, frame, stmt, true)
        end
        pc === nothing && break
    end
    return musteval[frame.pc] ? get_return(frame) : nothing
end


text:
"""

Update data or init R files.

$(SIGNATURES)

# Extended help

### Required arguments
```julia
* `m`                             : CmdStanModels object
* `input`                         : Input data or init values
* `num_chains`                    : Number of chains in model
* `fname_part="data"`             : Data or init R files to be created
```

Not exported.
"""
function update_R_files(m, input, num_chains, fname_part="data")

  m_field = fname_part == "data" ? m.data_file : m.init_file
  if typeof(input) <: NamedTuple || typeof(input) <: Dict
    for i in 1:num_chains
      stan_dump(m.output_base*"_$(fname_part)_$i.R", input, force=true)
      append!(m_field, [m.output_base*"_$(fname_part)_$i.R"])
    end
  elseif  typeof(input) <: Array
    if length(input) == num_chains
      for (i, d) in enumerate(input)
        stan_dump(m.output_base*"_$(fname_part)_$i.R", d, force=true)
        append!(m_field, [m.output_base*"_$(fname_part)_$i.R"])
      end
    else
      @info "Data vector length does not match number of chains,"
      @info "only first element in data vector will be used,"
      for i in 1:num_chains
        stan_dump(m.output_base*"_$(fname_part)_$i.R", input[1], force=true)
        append!(m_field, [m.output_base*"_$(fname_part)_$i.R"])
      end
    end
  elseif typeof(input) <: AbstractString && length(input) > 0
    for i in 1:num_chains
      cp(input, "$(m.output_base)_$(fname_part)_$i.R", force=true)
      append!(m_field, [m.output_base*"_$(fname_part)_$i.R"])
    end
  else
    error("\nUnrecognized input argument: $(typeof(input))\n")
  end

end


text:
using Revise
using ursa_functions
using Plots

ursa_0=URSA(0.95,3,0,2,0.8,10,100)  #parameter of URSA, see definition in ursa_functions
ursa_0.emit_prob_1=[0.4385 0.5615;0.5379 0.4621; 0.6860 0.3140] #opponent reaction probability

pbvi_solver=PBVI_solver(100,ursa_0)  #100 is the number of alpha vectors
solve!(pbvi_solver) #normal PBVI

N_sim=10000
a_v=zeros(Float64,N_sim)
for k=1:Int64(N_sim/2)
a_v[k]=simulate_true_reward_ursa(0.5,pbvi_solver,"non robust", "random", pbvi_solver.ursa)     #0.5 is the initial belief
end

pbvi_solver_rb=PBVI_solver(100,ursa_0)
robust_solve!(pbvi_solver_rb)  #robust PBVI

a_v=zeros(Float64,N_sim)
for k=1:Int64(N_sim/2)
a_v[k]=simulate_true_reward_ursa(0.5,pbvi_solver_rb,"robust", "norminal", pbvi_solver.ursa)
end

pbvi_solver_ch=PBVI_solver(100,ursa_0)
chance_solve!(pbvi_solver_ch,0.05,0.001)  #chance-constraint, 0.05 means 95% confidence

a_v=zeros(Float64,N_sim)
for k=1:Int64(N_sim/2)
a_v[k]=simulate_true_reward_ursa(0.5,pbvi_solver_ch,"chance", "norminal", pbvi_solver.ursa)
end


text:
using JuMP, Gurobi, LinearAlgebra, Suppressor, StatsBase, Compat, MathOptInterface

mutable struct cut
    intercept::Float64
    slope::Array{Float64,2}
end

mutable struct CutIterData #used for plotting performance of method
    time::Float64
    obj::Float64
    bound::Float64
end

# This is a QCQP formulation of the problem
function getMatrixFrobeniusNorm_exact(coords, entries, M, n, numEntriesSampled, k=sqrt(n), gamma=200*sqrt(n)/(numEntriesSampled), imposeConstraintsLazily=true, u_init=rand(n,k), v_init=rand(n,k))

    m=Model(Gurobi.Optimizer)
    set_optimizer_attribute(m, "NonConvex", 2)
    set_optimizer_attribute(m, "TimeLimit", 3600.0)
    set_optimizer_attribute(m, "MIPGap", 1e-2)
    @variable(m, u[1:n, 1:k])
    @variable(m, v[1:n, 1:k])
    @variable(m, X[1:n, 1:n])
    for l=1:(k-1)
        for l2=2:k
            @constraint(m, u[:,l]'*u[:,l2]<=1e-4) #orthogonality constraint, since columns are orthogonal in an svd decomposition.
            @constraint(m, u[:,l]'*u[:,l2]>=-1e-4)
            @constraint(m, v[:,l]'*v[:,l2]<=1e-4)
            @constraint(m, v[:,l]'*v[:,l2]>=-1e-4)
        end
    end
    @constraint(m, X.>=u*v'.-1e-6*ones(n,n))
    @constraint(m, X.<=u*v'.+1e-6*ones(n,n))
    @constraint(m, 0 .<= -X.+n*ones(n,n))
    @constraint(m, 0 .<= X.+n*ones(n,n))
    for i=1:k
        @constraint(m, [1.0; u[:,i]] in SecondOrderCone())
    end
    @constraint(m, 0 .<=u.+1.0*ones(n,k))
    @constraint(m, 0 .<=v.+n*ones(n,k))
    @constraint(m, 0 .<=-u.+1.0*ones(n,k))
    @constraint(m, 0 .<=-v.+n*ones(n,k))

    @objective(m, Min, 1.0/(2.0*gamma)*sum(X[i,j]^2 for i=1:n for j=1:n)+0.5*sum((X[coords[entries[t],1], coords[entries[t],2]]-M[coords[entries[t],1], coords[entries[t],2]])^2 for t in 1:Int(numEntriesSampled)))
    optimize!(m)

    return value.(u)*value.(v)', MOI.get(m, MOI.NodeCount()), objective_value(m)-objective_bound(m)
end

function getMatrixFrobeniusNorm_cuttingplanes_exact_inout(coords, entries, M, n, numEntriesSampled, k=sqrt(n), gamma=200*sqrt(n)/(numEntriesSampled), stabilizationPoint=k/n*Diagonal(ones(n)), X_ws=zeros(n,n), useInOut=true, useSOCs=true, useKelleys=false)
    m=Model(Mosek.Optimizer) #First solve the continuous relaxation
    set_optimizer_attribute(m, "MSK_DPAR_INTPNT_CO_TOL_PFEAS", 1e-10)
    set_optimizer_attribute(m, "MSK_DPAR_INTPNT_CO_TOL_DFEAS", 1e-10)
    @variable(m, Y[1:n, 1:n], PSD)
    @constraint(m, Symmetric(Matrix(1.0I, n, n)-Y) in PSDCone());
    @variable(m, theta>=-1e5)

    @objective(m, Min, theta)
    @constraint(m, sum(Y[i,i] for i=1:n)<=k)

    CurrentTime=time()

    ######### Start defining variables for the exact method ##############
    m_ex=Model(Gurobi.Optimizer)
    set_optimizer_attribute(m_ex, "NonConvex", 2)
    set_optimizer_attribute(m_ex, "TimeLimit", 30000.0)
    set_optimizer_attribute(m_ex, "MIPGap", 1e-4)
    set_optimizer_attribute(m_ex, "FuncPieceError", 1e-6)
    set_optimizer_attribute(m_ex, "FuncPieceLength", 1e-5)
    set_optimizer_attribute(m_ex, "Threads", 12)
    set_optimizer_attribute(m_ex, "Presolve", 0)

    @variable(m_ex, Y_ex[1:n, 1:n], Symmetric)
    @constraint(m_ex, sum(Y_ex[i,i] for i=1:n)<=k)
    @variable(m_ex, U_ex[1:n, 1:k])

    @variable(m_ex, theta_ex>=-1e4)
    @objective(m_ex, Min, theta_ex)

    @constraint(m_ex, defineY[i=1:n, j=i:n], Y_ex[i,j].<=U_ex[i,:]'*U_ex[j,:]+1e-4) #Note that we need only impose this in upper triangle, as Y is symmetric
    @constraint(m_ex, defineY2[i=1:n, j=i:n], Y_ex[i,j].>=U_ex[i,:]'*U_ex[j,:]-1e-4)

    bcdata = CutIterData[]

    for i=1:k
        for j=i:k
            ind=1.0*(i==j)
            @constraint(m_ex, U_ex[:,i]'*U_ex[:,j]<=1.0*ind+1e-6)
            @constraint(m_ex, U_ex[:,i]'*U_ex[:,j]>=1.0*ind-1e-6)
        end
    end
    if useSOCs
    # Strengthen the formulation using second-order cone inequalities on Y
        @constraint(m_ex, imposeSOC[i=1:n, j=1:n], [Y_ex[i,i]+Y_ex[j,j]; Y_ex[i,i]-Y_ex[j,j]; 2.0*Y_ex[i,j]] in SecondOrderCone())
        #@constraint(m_ex, revSOC[i=1:n, j=1:n], Y_ex[i,i]*Y_ex[j,j]<=Y_ex[i,j]^2+1e-4)
    # else # omitted, but could potentially use a diagonally dominant matrix approximation per Ahmadi/Hall
    #     @constraint(m_ex, imposeDD1[i=1:n, j=1:n], Y_ex[i,i]+Y_ex[j,j]+2.0*Y_ex[i,j]>=0.0)
    #     @constraint(m_ex, imposeDD2[i=1:n, j=1:n], Y_ex[i,i]+Y_ex[j,j]-2.0*Y_ex[i,j]>=0.0)
    end

    @constraint(m_ex, imposeDiag[i=1:n], 0.0<=Y_ex[i,i]<=1.0)

    theIndices=zeros(n,n)
    for t=1:Int(numEntriesSampled)
        theIndices[coords[entries[t],1], coords[entries[t],2]]=1.0
    end

    ############ Start the cut loop #######################
    theCutPool=cut[]
    UB=1e12
    LB=-1e12
    rootStabilizationTrick = :inOut
     #Use the solution to the SDP relaxation instead, where appropriate
    rootCutsSense=1
    ε= 1e-10
    λ = (rootStabilizationTrick == :inOut) ? .1 : 1.
    δ = (rootStabilizationTrick == :inOut || rootStabilizationTrick == :twoEps) ? 2*ε : 0.
    rootCutsLim=200
    rootCutCount = 0
    oaRootCutCount=0
    consecutiveNonImprov_1 = 0
    consecutiveNonImprov_2 = 0
    if useInOut
        Y_best=zeros(n,n)
        for epoch in 1:200
          @suppress optimize!(m)
          Ystar=value.(Y)
          stabilizationPoint += Ystar; stabilizationPoint /= 2

          if LB >= objective_value(m) - eps()
            if consecutiveNonImprov_1 == 5
              consecutiveNonImprov_2 += 1
            else
              consecutiveNonImprov_1 += 1
            end
          else
            if consecutiveNonImprov_1 < 5
              consecutiveNonImprov_1 = 0
            end
            consecutiveNonImprov_2 = 0
          end
          LB = max(LB, objective_value(m))

          if consecutiveNonImprov_1 == 5
            λ = 1
          elseif consecutiveNonImprov_2 == 5
            δ = 0.
          end

          Y0 = λ*Ystar + (1-λ)*stabilizationPoint .+ δ*Diagonal(ones(n))
          if useKelleys #override in-out bit of method
              Y0=Ystar
          end

          # Perform an SVD step here to ensure still in convex hull
          u,v,=svd(Y0)
          v=v*k/sum(v)
          v.=min.(ones(n), max.(v, zeros(n)))

          Y0=u*Diagonal(v)*u'

          # Get a cut by solving the dual subproblem
          # Solve subproblem for fixed Y, in alpha
          f_Y, alpha_t=getmatrixcompcut(Y0, theIndices, M, n, gamma)
          if f_Y<UB
              Y_best=Y0
              UB=f_Y
          end
          LB=max(LB, objective_value(m))
          # Impose constraint
          @constraint(m, theta>=f_Y-gamma/2*Compat.dot(alpha_t*alpha_t', Y-Y0))

          @constraint(m_ex, theta_ex>=f_Y-gamma/2*Compat.dot(alpha_t*alpha_t', Y_ex-Y0))
          # Add cut to pool (add in later if it works)
          @show LB, UB, epoch, time()-CurrentTime
          if abs(UB-LB) <= 1e-4 || consecutiveNonImprov_2 >= 10
              break
          end
        end
    end

    # After the cutting-plane method on the continuous relaxation has terminated, apply an outer-approximation method on the exact problem
    # Using lazy callbacks (faster than a multi-tree method)
    cutCount=0
    maxCutCount=1e6
    best_UB=1e6
    function add_lazy_callback(cb_data)
        if cutCount<=maxCutCount
            Y0=zeros(n,n)
            for i=1:n
                for j=1:n
                    Y0[i,j]=callback_value(cb_data, Y_ex[i,j])
                end
            end
            if eigvals(Y0)[1]<=-0.1
                u_t=eigvecs(Y0)[:,1]
                con3 = @build_constraint(Compat.dot(Y_ex,u_t*u_t') >=0.0)
                MOI.submit(m_ex, MOI.LazyConstraint(cb_data), con3)
                noPSDCuts+=1
            end
            U0=zeros(n,k)
            for i=1:n
                for j=1:k
                    U0[i,j]=callback_value(cb_data, U_ex[i,j])
                end
            end

            for j=1:k
                if norm(U0[:,j])>1e-4
                    U0[:,j]=U0[:,j]/norm(U0[:,j])
                else
                    U0[1,j]=1.0
                    U0[:,j]=U0[:,j]/norm(U0[:,j])
                end
            end
            #@show U0'*U0
            Y0=U0*U0'

            f_Y, alpha_t=getmatrixcompcut(Y0, theIndices, M, n, gamma)
            if f_Y<best_UB
                @show best_UB=f_Y
            end
            b_t=f_Y+gamma/2*Compat.dot(alpha_t*alpha_t', Y0)
            # f_0=b_t-floor(b_t)

            A_t=(gamma/2)*alpha_t*alpha_t'#+1e-3*Diagonal(ones(n))

            con = @build_constraint(theta_ex>=b_t-Compat.dot(A_t, Y_ex))
            MOI.submit(m_ex, MOI.LazyConstraint(cb_data), con)
            #MOI.submit(m_ex, MOI.LazyConstraint(cb_data), con2)
            cutCount+=1
        end
    end
    MOI.set(m_ex, MOI.LazyConstraintCallback(), add_lazy_callback)

    noPSDCuts=0
    maxPSDCuts=0
    function addPSDCuts(cb_data)
        if noPSDCuts<=maxPSDCuts
            Y0=zeros(n,n)
            for i=1:n
                for j=1:n
                    Y0[i,j]=callback_value(cb_data, Y_ex[i,j])
                end
            end
            if eigvals(Y0)[1]<=-0.1
                u_t=eigvecs(Y0)[:,1]
                con = @build_constraint(Compat.dot(Y_ex,u_t*u_t') >=0.0)
                MOI.submit(m_ex, MOI.UserCut(cb_data), con)
                noPSDCuts+=1
            end
        end
    end
    #MOI.set(m_ex, MOI.UserCutCallback(), addPSDCuts)

    U_ws,=svd(X_ws) #Compute Y_ws given X_ws, then use SVD to make sure it is indeed of rank at most k
    U_ws=U_ws[:, 1:k] # Julia automatically orders the columns of U in terms of the size of their singular values

    Y_ws=U_ws*U_ws'

    f_Y, alpha_t=getmatrixcompcut(Y_ws, theIndices, M, n, gamma)

    @constraint(m_ex, theta_ex>=f_Y-gamma/2*Compat.dot(alpha_t*alpha_t', Y_ex-Y_ws))

    f_Y_ws=f_Y

    # set_start_value(theta_ex, f_Y_ws)
    set_start_value.(U_ex, U_ws)
    set_start_value.(Y_ex, Y_ws)

    optimize!(m_ex)

    @show cutCount, best_UB

    Y_ex0=zeros(n,n)
    status = MOI.get(m_ex, MOI.TerminationStatus())
    primalstatus=MOI.get(m_ex, MOI.PrimalStatus())
    nodecount= MOI.get(m_ex, MOI.NodeCount())
    objval= -1.0
    objbound= 0.0

    if status == MOI.OPTIMAL || primalstatus==MOI.FEASIBLE_POINT
        # Set Y_ex0 equal to incumbent
        Y_ex0.=value.(Y_ex)
        objval=objective_value(m_ex)
        objbound=objective_bound(m_ex)
    else
        # No incumbent, set Y_ex0 equal to warm start (Gurobi hasn't found a solution due to numerical resolution issues), keep obj val and obj bound such that gap is -1.0
        Y_ex0.=Y_ws
    end

    # Recover X_ex at the very end (alla rounding)

    # Apply greedy rounding mechanism to purge any accumulated numerical errors
    u,v,w=LinearAlgebra.svd(Y_ex0)

    diag_v=zeros(n)
    diag_v[sortperm(v, rev=true)[1:k]].=1.0
    Y_rounded=u*Diagonal(diag_v)*u'

    # solve for X_ex

    m2=Model(Mosek.Optimizer)
    @variable(m2, X2[1:n, 1:n])
    @objective(m2, Min, 1.0/(2.0*gamma)*sum(X2[i,j]^2 for i=1:n for j=1:n)+sum((X2[coords[entries[t],1], coords[entries[t],2]]-M[coords[entries[t],1], coords[entries[t],2]])^2 for t in 1:Int(numEntriesSampled)))
    @constraint(m2, X2.==Y_rounded*X2)
    @suppress optimize!(m2)
    X_best=value.(X2)

   return X_best, nodecount, cutCount, objval-objbound
   # Note: LB from the in-out method, which may not always converge to the sdo bound; UB from the greedy rounding which may sometimes differ from Gurobi's bound
end

function getMatrixFrobeniusNorm_cuttingplanesmultitree_exact_inout(coords, entries, M, n, numEntriesSampled, k=sqrt(n), gamma=200*sqrt(n)/(numEntriesSampled), stabilizationPoint=k/n*Diagonal(ones(n)), X_ws=zeros(n,n), useInOut=true, fileName="")
    m=Model(Mosek.Optimizer) #First solve the continuous relaxation
    set_optimizer_attribute(m, "MSK_DPAR_INTPNT_CO_TOL_PFEAS", 1e-10)
    set_optimizer_attribute(m, "MSK_DPAR_INTPNT_CO_TOL_DFEAS", 1e-10)
    @variable(m, Y[1:n, 1:n], PSD)
    @constraint(m, Symmetric(Matrix(1.0I, n, n)-Y) in PSDCone());
    @variable(m, theta>=-1e4)
    @objective(m, Min, theta)
    @constraint(m, sum(Y[i,i] for i=1:n)<=k)

    m_ex=Model(Gurobi.Optimizer) #Second solve the exact problem
    set_optimizer_attribute(m_ex, "NonConvex", 2)
    set_optimizer_attribute(m_ex, "TimeLimit", 30.0)
    set_optimizer_attribute(m_ex, "MIPGap", 1e-2)
    set_optimizer_attribute(m_ex, "FuncPieceError", 1e-6)
    set_optimizer_attribute(m_ex, "FuncPieceLength", 1e-5)
    set_optimizer_attribute(m_ex, "Threads", 12)
    set_optimizer_attribute(m_ex, "Presolve", 0)
    #set_optimizer_attribute(m_ex, "Heuristics", 0)
    set_optimizer_attribute(m_ex, "FuncMaxVal", 1e4)
    # set_optimizer_attribute(m_ex, "ResultFile", "theResultFile.sol")

    @variable(m_ex, Y_ex[1:n, 1:n], Symmetric)
    @constraint(m_ex, sum(Y_ex[i,i] for i=1:n)==k)
    @variable(m_ex, U_ex[1:n, 1:k])

    @constraint(m_ex, strengthenY0[i=1:n], Y_ex[i,i]>=sum(U_ex[i,t]^2 for t=1:k)) #Strengthen the relaxation
    @constraint(m_ex, strengthenY1[i=1:n, j=1:n], [Y_ex[i,i]+Y_ex[j,j]-2.0*Y_ex[i,j]+1.0;Y_ex[i,i]+Y_ex[j,j]-2.0*Y_ex[i,j]-1.0; U_ex[i,:]-U_ex[j,:]] in SecondOrderCone())
    @constraint(m_ex, strengthenY2[i=1:n, j=1:n], [Y_ex[i,i]+Y_ex[j,j]+2.0*Y_ex[i,j]+1.0;Y_ex[i,i]+Y_ex[j,j]+2.0*Y_ex[i,j]-1.0; U_ex[i,:]+U_ex[j,:]] in SecondOrderCone())

    @constraint(m_ex, defineY[i=1:n, j=i:n], Y_ex[i,j].<=U_ex[i,:]'*U_ex[j,:]+1e-4)
    @constraint(m_ex, defineY2[i=1:n, j=i:n], Y_ex[i,j].>=U_ex[i,:]'*U_ex[j,:]-1e-4)

    for i=1:k
        for j=i:k
            ind=1.0*(i==j)
            @constraint(m_ex, U_ex[:,i]'*U_ex[:,j]>=1.0*ind-1e-4)
            @constraint(m_ex, U_ex[:,i]'*U_ex[:,j]<=1.0*ind+1e-4)
        end
    end
    # # Strengthen the formulation using second-order cone inequalities on Y
    @constraint(m_ex, imposeSOC[i=1:n, j=1:n], [Y_ex[i,i]+Y_ex[j,j]; Y_ex[i,i]-Y_ex[j,j]; 2.0*Y_ex[i,j]] in SecondOrderCone())
    @constraint(m_ex, imposeDiag[i=1:n], 0.0<=Y_ex[i,i]<=1.0)
    @variable(m_ex, theta_ex>=-1e4)
    @objective(m_ex, Min, theta_ex)

    theIndices=zeros(n,n)
    for t=1:Int(numEntriesSampled)
        theIndices[coords[entries[t],1], coords[entries[t],2]]=1.0
    end

    UB=1e12
    LB=-1e12
    rootStabilizationTrick = :inOut
     #Use the solution to the SDP relaxation instead, where appropriate
    rootCutsSense=1
    ε= 1e-10
    λ = (rootStabilizationTrick == :inOut) ? .1 : 1.
    δ = (rootStabilizationTrick == :inOut || rootStabilizationTrick == :twoEps) ? 2*ε : 0.
    rootCutsLim=200
    rootCutCount = 0
    oaRootCutCount=0
    consecutiveNonImprov_1 = 0
    consecutiveNonImprov_2 = 0
    # theCutPool=Cut[]
    if useInOut
        Y_best=zeros(n,n)
        for epoch in 1:200
          @suppress optimize!(m)
          Ystar=value.(Y)
          stabilizationPoint += Ystar; stabilizationPoint /= 2

          if LB >= objective_value(m) - eps()
            if consecutiveNonImprov_1 == 5
              consecutiveNonImprov_2 += 1
            else
              consecutiveNonImprov_1 += 1
            end
          else
            if consecutiveNonImprov_1 < 5
              consecutiveNonImprov_1 = 0
            end
            consecutiveNonImprov_2 = 0
          end
          LB = max(LB, objective_value(m))

          if consecutiveNonImprov_1 == 5
            λ = 1
          elseif consecutiveNonImprov_2 == 5
            δ = 0.
          end

          Y0 = λ*Ystar + (1-λ)*stabilizationPoint .+ δ*Diagonal(ones(n))

          # Perform an SVD step here to ensure still in convex hull
          u,v,=svd(Y0)
          v=v*k/sum(v)
          v.=min.(ones(n), max.(v, zeros(n)))

          Y0=u*Diagonal(v)*u'

          f_Y, alpha_t=getmatrixcompcut(Y0, theIndices, M, n, gamma)
          if f_Y<UB
              Y_best=Y0
              UB=f_Y
          end
          LB=max(LB, objective_value(m))
          # Impose constraint
          @constraint(m, theta>=f_Y-gamma/2*Compat.dot(alpha_t*alpha_t', Y-Y0))
          # Also impose constraint in exact master problem (this is the whole point of the cut loop)
          @constraint(m_ex, theta_ex>=f_Y-gamma/2*Compat.dot(alpha_t*alpha_t', Y_ex-Y0))
          # Add cut to pool (add in later if it works)
          @show LB, UB, epoch
          if abs(UB-LB) <= 1e-4 || consecutiveNonImprov_2 >= 10
              break
          end
        end
    end

    U_ws,=svd(X_ws) #Compute Y_ws given X_ws, then use SVD to make sure it is indeed of rank at most k
    U_ws=U_ws[:, 1:k] # Julia automatically orders the columns of U in terms of the size of their singular values
    Y_ws=U_ws*U_ws'
    # @show Y_ws[1,1]
    # @show typeof(Y_ws[1,1])

    set_start_value.(Y_ex, Y_ws)
    set_start_value.(U_ex, U_ws)

    f_Y, alpha_t=getmatrixcompcut(Y_ws, theIndices, M, n, gamma)
    @constraint(m_ex, theta_ex>=f_Y-gamma/2*Compat.dot(alpha_t*alpha_t', Y_ex-Y_ws))

    noPSDCuts=0
    maxPSDCuts=100
    function addPSDCuts(cb_data)
        if noPSDCuts<=maxPSDCuts
            Y0=zeros(n,n)
            for i=1:n
                for j=1:n
                    Y0[i,j]=callback_value(cb_data, Y_ex[i,j])
                end
            end

            if eigvals(Y0)[1]<=-0.1
                u_t=eigvecs(Y0)[:,1]
                con = @build_constraint(Compat.dot(Y_ex,u_t*u_t') >=0.0)
                MOI.submit(m_ex, MOI.UserCut(cb_data), con)
                noPSDCuts+=1
            end
        end
    end
    MOI.set(m_ex, MOI.UserCutCallback(), addPSDCuts)

    UB_master=1e12
    LB_master=-1e12
    cutCount=0
    nc_agg=0
    Y_best=zeros(n,n)
    maxOAIters=20
    epochData = DataFrame(n=Int64[], gamma=Float64[], iter=Int64[], nodes=Float64[])

    theNodeData = similar(epochData, 0)

    consecNonImprove=0 #Use to decide if we need to perform a more expensive iteration, since the lower bound isn't changing
    for t=1:maxOAIters
        if t>1
            set_optimizer_attribute(m_ex, "MIPGap", 1e-1/t)
            set_optimizer_attribute(m_ex, "Threads", 12)
            set_optimizer_attribute(m_ex, "Presolve", 0)
            noPSDCuts=0
            set_optimizer_attribute(m_ex, "TimeLimit", min(10.0*t, 300.0))
            set_optimizer_attribute(m_ex, "PoolSolutions", 2) #No point in keeping a large number of solutions; just slows down the next iteration
            if t<=2
                set_start_value.(Y_ex, value.(Y_ex))
                set_start_value.(U_ex, value.(U_ex))
            end
            if consecNonImprove>=2
                set_optimizer_attribute(m_ex, "TimeLimit", 300.0)  #Expensive iteration, to move the bound
                consecNonImprove=0
            end

        end

        @time optimize!(m_ex)
        U_iter=value.(U_ex)
        Y_iter=value.(Y_ex)
        # Impose soc callback cut
        if eigvals(Y_iter-U_iter*U_iter')[1]<=-0.01
            u_t=eigvecs(Y_iter-U_iter*U_iter')[:,1]
            @constraint(m_ex, [u_t'*Y_ex*u_t+1.0, u_t'*Y_ex*u_t-1.0, 2.0*U_iter'*u_t] in SecondOrderCone())
        end
        # Impose OA cut at the exact value identified by Gurobi
        f_Y, alpha_t=getmatrixcompcut(Y_iter, theIndices, M, n, gamma)
        @constraint(m_ex, theta_ex>=f_Y-gamma/2*Compat.dot(alpha_t*alpha_t', Y_ex-Y_iter))

        # Now round the matrix to ensure that it is infact an orthogonal projection matrix, and impose a cut there too.
        u,v,=svd(Y_iter)
        Y_iter=u[:,1:k]*u[:,1:k]'
        LB_diff=objective_bound(m_ex)-LB_master
        LB_master=max(LB_master, objective_bound(m_ex))
        f_Y, alpha_t=getmatrixcompcut(Y_iter, theIndices, M, n, gamma)
        @constraint(m_ex, theta_ex>=f_Y-gamma/2*Compat.dot(alpha_t*alpha_t', Y_ex-Y_iter))
        if f_Y<UB_master
            UB_master=f_Y
            Y_best=Y_iter
        elseif LB_diff<1e-3 #Neither bound moved, so need a more expensive iteration
            consecNonImprove+=1
        end
        @show typeof(n)
        @show typeof(gamma)
        @show typeof(cutCount)
        @show typeof(MOI.get(m_ex, MOI.NodeCount()))
        @show size(MOI.get(m_ex, MOI.NodeCount()))
        nc_agg+=MOI.get(m_ex, MOI.NodeCount())
        push!(theNodeData, [n, gamma, cutCount, MOI.get(m_ex, MOI.NodeCount())])
        cutCount+=1
        @show LB_master, UB_master, t
        @show MOI.get(m_ex, MOI.NodeCount())
        if abs(UB_master-LB_master) <= 1e-2
            break
        end
    end

    CSV.write("nodesvsiters"*fileName*".csv", theNodeData)

    @show cutCount
    @show MOI.get(m_ex, MOI.NodeCount())
    # Recover X_ex at the very end (alla rounding)
    Y_ex0=Y_best#value.(Y_ex)
    # Apply greedy rounding mechanism to purge any accumulated numerical errors
    u,v,w=LinearAlgebra.svd(Y_ex0)

    diag_v=zeros(n)
    diag_v[sortperm(v, rev=true)[1:k]].=1.0
    Y_rounded=u*Diagonal(diag_v)*u'

    # solve for X_ex

    m2=Model(Mosek.Optimizer)
    @variable(m2, X2[1:n, 1:n])
    @objective(m2, Min, 1.0/(2.0*gamma)*sum(X2[i,j]^2 for i=1:n for j=1:n)+sum((X2[coords[entries[t],1], coords[entries[t],2]]-M[coords[entries[t],1], coords[entries[t],2]])^2 for t in 1:Int(numEntriesSampled)))
    @constraint(m2, X2.==Y_rounded*X2)
    @suppress optimize!(m2)
    X_best=value.(X2)

   return X_best, MOI.get(m_ex, MOI.NodeCount()), cutCount, objective_value(m_ex)-max(objective_bound(m_ex), objective_bound(m)), nc_agg
   # Note: LB from the in-out method, which may not always converge to the sdo bound; UB from the greedy rounding which may sometimes differ from Gurobi's bound
end

function getmatrixcompcut(Y0::Array{Float64, 2}, theIndices::Array{Float64, 2}, M::Array{Float64, 2}, n::Int64, gamma::Float64)
    u,v=svd(0.5*Y0+0.5*Y0')

    u0=u*Diagonal(sqrt.(v))
    m2=Model(Mosek.Optimizer)
    set_optimizer_attribute(m2, "MSK_DPAR_INTPNT_CO_TOL_PFEAS", 1e-8)
    set_optimizer_attribute(m2, "MSK_DPAR_INTPNT_CO_TOL_DFEAS", 1e-8)
    @variable(m2, alpha[1:n, 1:n])
    @variable(m2, r[1:n, 1:n])
    @constraint(m2, r.==alpha'*u0)
    @constraint(m2, alpha.*(ones(n,n)-theIndices).==0.0) # If not in set of indices then eliminate variable
    @objective(m2, Max, -0.5*sum((alpha[i,j]+M[i,j])^2*theIndices[i,j] for i=1:n for j=1:n)-(gamma/2.0)*Compat.dot(r,r))
    @suppress optimize!(m2)
    alpha_t=value.(alpha).*theIndices
    # Get UB and cut from second-stage optimization problem
    f_Y=objective_value(m2)+0.5*sum((M.*theIndices).^2)
    return f_Y, alpha_t
end
# Pareto-optimal version, assuming we have access to k
function getmatrixcompcut_pareto(Y0::Array{Float64, 2}, theIndices::Array{Float64, 2}, M::Array{Float64, 2}, n::Int64, gamma::Float64, k::Int64)
    u,v=svd(Y0)
    u_inner=(sqrt(k)/sqrt(n))*Diagonal(ones(n)) #For the Pareto-optimal cut

    u0=u*Diagonal(sqrt.(v))
    m2=Model(Mosek.Optimizer)
    set_optimizer_attribute(m2, "MSK_DPAR_INTPNT_CO_TOL_PFEAS", 1e-8)
    set_optimizer_attribute(m2, "MSK_DPAR_INTPNT_CO_TOL_DFEAS", 1e-8)
    @variable(m2, alpha[1:n, 1:n])
    @variable(m2, r[1:n, 1:n])
    @constraint(m2, r.==alpha'*u0)
    @constraint(m2, alpha.*(ones(n,n)-theIndices).==0.0) # If not in set of indices then eliminate variable
    @objective(m2, Max, -0.5*sum((alpha[i,j]+M[i,j])^2*theIndices[i,j] for i=1:n for j=1:n)-(gamma/2.0)*Compat.dot(r,r))
    @suppress optimize!(m2)
    # Get UB and cut from second-stage optimization problem
    f_Y=objective_value(m2)+0.5*sum((M.*theIndices).^2)

    m3=Model(Mosek.Optimizer)
    set_optimizer_attribute(m3, "MSK_DPAR_INTPNT_CO_TOL_PFEAS", 1e-8)
    set_optimizer_attribute(m3, "MSK_DPAR_INTPNT_CO_TOL_DFEAS", 1e-8)
    @variable(m3, alpha2[1:n, 1:n])
    @variable(m3, r2[1:n, 1:n])
    @constraint(m3, r2.==alpha2'*u0)
    @variable(m3, r3[1:n, 1:n])
    @constraint(m3, r3.==alpha2'*u_inner)
    @constraint(m3, alpha2.*(ones(n,n)-theIndices).==0.0) # If not in set of indices then eliminate variable
    @constraint(m3, -0.5*sum((alpha2[i,j]+M[i,j])^2*theIndices[i,j] for i=1:n for j=1:n)-(gamma/2.0)*Compat.dot(r2,r2)>=objective_value(m2)) #Ensures that the cut supplies the right OFV
    @objective(m3, Max, -0.5*sum((alpha2[i,j]+M[i,j])^2*theIndices[i,j] for i=1:n for j=1:n)-(gamma/2.0)*Compat.dot(r3,r3))
    @suppress optimize!(m3)
    alpha_t=value.(alpha2).*theIndices
    return f_Y, alpha_t
end

max_stars_repo_path: Matrix Completion/exact_approaches.jl

max_stars_repo_name: ryancorywright/MixedProjectionSoftware

max_issues_repo_path: Matrix Completion/exact_approaches.jl

max_issues_repo_name: ryancorywright/MixedProjectionSoftware

max_forks_repo_path: Matrix Completion/exact_approaches.jl

max_forks_repo_name: ryancorywright/MixedProjectionSoftware