text:
"""
    SPAMethod

Abstract type for nesting the various methods that can be used to perform an SPA test. Subtypes are: \n
    SPABoot \n
The subtypes have entries in the help (?) menu.
"""
abstract type SPAMethod end

"""
    SPABoot(data ; alpha, kernelfunction, bandwidth, bootinput_kwargs...)

Method type for doing an SPA test using a dependent bootstrap. A constructor
that takes the data and several keyword arguments is provided. The data refers
to the loss differences for each forecast relative to the basecase. Use ?spa for more detail.
Relevant keyword arguments are: \n
    alpha=0.05 <- The confidence level to use in the test \n
    kernelfunction=KernelEpanechnikov() <- The kernel function to use with HAC variance estimator. See ?hacvariance for more detail. \n
	bandwidth=-1 <- The bandwidth for HAC variance estimator. If bandwidth <= -1 then bandwidth is estimated using Politis (2003) "Adaptive Bandwidth Choice" \n
    blocklength=0.0 <- Block length to use with the bootstrap. Default of 0.0 implies
        the block length will be optimally estimated from the data use the method deemed
        most appropriate by the DependentBootstrap package (typically the selection procedure
        of Patton, Politis, and White (2009)). \n
    numresample=1000 <- Number of resamples to use when bootstrapping. \n
    bootmethod=:stationary <- Bootstrap methodology to use, where the default is the
        stationary bootstrap of Politis and Romano (1993) \n
Note, see the DependentBootstrap docs for more information on valid keyword arguments
for the BootInput constructor.
"""
struct SPABoot{Tk<:KernelFunction} <: SPAMethod
	alpha::Float64
	bootinput::BootInput
	kernelfunction::Tk
	bandwidth::Int
	function SPABoot(alpha::Float64, bootinput::BootInput, kernelfunction::Tk, bandwidth::Int) where {Tk<:KernelFunction}
		!(0.0 < alpha < 0.5) && error("Confidence level set to $(alpha) which is not on the (0, 0.5) interval")
		new{Tk}(alpha, bootinput, kernelfunction, bandwidth)
	end
end
function SPABoot(data ; alpha::Float64=0.05, kernelfunction=KernelEpanechnikov(), bandwidth::Int=-1, kwargs...)
    return SPABoot(alpha, BootInput(data ; flevel1=mean, kwargs...), get_kernel_type(kernelfunction), bandwidth)
end

"""
	SPATest(rejH0::Int, mu_u::Vector{Float64}, mu_c::Vector{Float64}, mu_l::Vector{Float64}, pvalue_u::Float64, pvalue_c::Float64, pvalue_l::Float64, pvalue::Float64, teststat::Float64)

Output type for an SPA test proposed in Hansen (2005) "A Test for Superior Predictive Ability". The names of the fields
for this type follow those in the referenced paper, so see that paper for more detail. The fields of this type follow: \n
	rejH0 <- true if the null is rejected, false otherwise
	mu_u <- mu associated with upper bound on threshold rate
	mu_c <- Recommended mu
	mu_l <- mu associated with lower bound on threshold rate
	pvalue_u <- pvalue associated with mu_u
	pvalue_c <- pvalue assoicated with mu_c
	pvalue_l <- pvalue associated with mu_l
	pvalueauto <- Recommended p-value (usually pvalue_c except in rare situations where this method can fail)
	teststat::Float64 <- SPA test statistic (look for tSPA in source paper) \n
Note the different mu variables, {u, c, l} are described in the referenced article on p370-371. \n
Most users will just be interested in the field pvalueauto.
"""
struct SPATest
	rejH0::Bool
	mu_u::Vector{Float64} #mu associated with upper bound on threshold rate
	mu_c::Vector{Float64} #Recommended mu
	mu_l::Vector{Float64} #mu associated with lower bound on threshold rate
	pvalue_u::Float64 #pvalue associated with mu_u
	pvalue_c::Float64 #pvalue assoicated with mu_c
	pvalue_l::Float64 #pvalue associated with mu_l
	pvalueauto::Float64 #Recommended p-value (usually pvalue_c except in rare situations)
	teststat::Float64 #SPA test statistic
end
function Base.show(io::IO, a::SPATest)
    println("SPA test results:")
    println("    test stat: $(a.teststat)")
    println("    pvalue (u): $(a.pvalue_u)")
    println("    pvalue (c): $(a.pvalue_c)")
    println("    pvalue (l): $(a.pvalue_l)")
end

"""
	spa{T<:Number}(lossDiff::Matrix{T}, method ; kwargs...)::SPATest
	spa{T<:Number}(lossDiff::Matrix{T} ; kwargs...)

This function implements the SPA test proposed in Hansen (2005) "A Test for Superior Predictive Ability". \n
Let x_0 denote a base-case forecast, x_k, k = 1, ..., K, denote K alternative forecasts, and y denote the forecast target.
Let L(., .) denote a loss function. The first argument of spa is a matrix where the kth column of the matrix is created by the operation: \n
L(x_k, y) - L(x_0, y) \n
Note that the forecast loss comes first and the base case loss comes second. This is the opposite to what is described in Hansen's paper. \n
The second method argument determines which methodology to use. Currently, only SPABoot is available and if
this input type is provided, the keyword arguments are not needed.
Alternatively, the user can omit the second argument, and then any keyword arguments will be passed to the
SPABoot constructor. See ?SPABoot for more detail. The most relevant keyword arguments are:
    alpha=0.05 <- The confidence level to use in the test \n
    kernelfunction=KernelEpanechnikov() <- The kernel function to use with HAC variance estimator. See ?hacvariance for more detail. \n
	bandwidth=-1 <- The bandwidth for HAC variance estimator. If bandwidth <= -1 then bandwidth is estimated using Politis (2003) "Adaptive Bandwidth Choice" \n
    blocklength=0.0 <- Block length to use with the bootstrap. Default of 0.0 implies
        the block length will be optimally estimated from the data use the method deemed
        most appropriate by the DependentBootstrap package (typically the selection procedure
        of Patton, Politis, and White (2009)). \n
    numresample=1000 <- Number of resamples to use when bootstrapping. \n
    bootmethod=:stationary <- Bootstrap methodology to use, where the default is the
        stationary bootstrap of Politis and Romano (1993) \n
The output of an SPA test is of type SPATest. Use ?SPATest for more information. \n
Note, Hansen suggests using the Stationary Bootstrap implied HAC variance estimator, which is not currently supported in this package.
However, note that any consistent HAC estimator is valid and in many cases may be preferred.
"""
function spa(lD::Matrix{<:Number}, method::SPABoot)::SPATest
	lD = -1.0*lD	#Hansen's loss differentials have base case first
	numObs = size(lD, 1)
	numModel = size(lD, 2)
	numResample = method.bootinput.numresample
	numObs < 2 && error("Number of observations = $(numObs) which is not enough to perform a reality check")
	numModel < 1 && error("Input dataset is empty")
	inds = dbootinds(Float64[], method.bootinput) #Get bootstrap indices
	#Get hac variance estimators
	wSqVec = Array{Float64}(undef, numModel)
	for k = 1:numModel
		(wSqVec[k], _) = hacvariance(lD[:, k], kf=method.kernelfunction, bw=method.bandwidth)
		wSqVec[k] <= 0.0 && (wSqVec[k] = nextfloat(0.0))
	end
	wInv = Float64[ 1 / sqrt(wSqVec[k]) for k = 1:numModel ]
	#Get bootstrapped mean loss differentials
	mldBoot = Array{Float64}(undef, numModel, numResample)
	for j = 1:numResample
		for k = 1:numModel
			mldBoot[k, j] = mean(view(lD, 1:numObs, k)[inds[j]])
		end
	end
	#Get mu definitions
	mu_u = Float64[ mean(view(lD, 1:numObs, k)) for k = 1:numModel ]
	mu_l = Float64[ max(mu_u[k], 0) for k = 1:numModel ]
	multTerm_mu_c = (1/numObs) * 2 * log(log(numObs))
	mu_c = Float64[ (mu_u[k] >= -1 * sqrt(multTerm_mu_c * wSqVec[k])) * mu_u[k] for k = 1:numModel ]
	#Get test statistic
	tSPA = maximum([ max(sqrt(numObs) * mu_u[k] * wInv[k], 0) for k = 1:numModel ])
	#Get p-values for each mu term
	muVecVec = Vector{Float64}[mu_u, mu_c, mu_l] #Order must be u, c, l
	pValVec = Array{Float64}(undef, length(muVecVec))
	for q = 1:length(muVecVec)
		z = sqrt(numObs) * (wInv .* (mldBoot .- muVecVec[q]))
		tSPAmu = Float64[ max(0, maximum(view(z, 1:numModel, i))) for i = 1:numResample ]
		pValVec[q] = (1 / numResample) * sum(tSPAmu .> tSPA)
	end
	any(mu_c .!= 0.0) ? (pvalueAuto = pValVec[2]) : (pvalueAuto = pValVec[3]) #Ad hoc rule as to which p-value is likely to be the most reliable
	pvalueAuto < method.alpha ? (rejH0 = true) : (rejH0 = false)
	spaTest = SPATest(rejH0, mu_u, mu_c, mu_l, pValVec[1], pValVec[2], pValVec[3], pvalueAuto, tSPA)
	return(spaTest)
end
#Keyword wrapper
spa(lossDiff::Matrix{<:Number} ; kwargs...)::SPATest = spa(lossDiff, SPABoot(lossDiff ; kwargs...))


text:
using LinearAlgebra
using GaussianDistributions
using SARA
using SARA: random_sampling, uncertainty_sampling, outer_sampling,
        TemperatureProfile, stripe_to_global
using ForwardDiff
using Base.Threads
using SARA: oSARA

function optical_gradient_map(x::AbstractMatrix, y::AbstractMatrix, σ²::AbstractMatrix)
    subsampling = 1
    stride = 1
    subsample(x, k::Int) = @view x[end+1-stride*k:stride:end, 1:subsampling:end]
    T_length = [30] # range(20, 35, step = 5) # different T length scales
    nsamples = [16] # vcat(1, 2, 4, 6, 8, 12, 16)
    k_scale = 1/4 # kernel scaling hyper-parameter
    σ² ./= k_scale^2
    τ_length = [.4] # τ length scales
    n = 512
    T = range(200, 1200, length = n + 1) # plotting ranges
    τ = range(2.5, 4, length = n)
    gradients = zeros(length(T), length(τ), length(nsamples), length(T_length), length(τ_length))
    for i in eachindex(T_length)
        T_l = T_length[i]
        for (j, τ_l) in enumerate(τ_length)
            @threads for ik in eachindex(nsamples) # number of gradient samples per stripe
                println(i, j, ik)
                k = nsamples[ik]
                xk = subsample(x, k)
                yk = subsample(y, k)
                σ²k = subsample(σ², k)

                # kernel = k_scale^2 * oSARA([T_l, τ_l])
                kernel = oSARA([T_l, τ_l])
                G = Gaussian(kernel)
                xk = vec(xk)
                yk = vec(yk)
                # μ_yk = mean(yk)
                # yk = yk .- μ_yk
                σ²k = vec(σ²k)
                @time C = conditional(G, xk, Gaussian(yk, Diagonal(σ²k)), tol = 1e-6)
                # C = C + μ_yk
                gij = @view gradients[:, :, ik, i, j]
                # gij = vec(gij)
                # @time gij .= mean(C, tuple.(T, τ'))
                @time gij .= mean(C).(tuple.(T, τ'))
            end
        end
    end
    if true
        f = h5open("outer_optical_gradient_map_stride_1.h5", "w")
        f["temperatures"] = collect(T)
        f["dwelltimes"] = collect(τ)
        f["gradients"] = gradients
        f["nsample"] = collect(nsamples)
        f["temperature lengthscales"] = collect(T_length)
        f["dwelltime lengthscales"] = collect(τ_length)
        f["kernel scale"] = k_scale
        f["dimensions"] = ["temperature", "dwelltime", "nsample",
                            "temperature lengthscale", "dwelltime lengthscale"]
        f["subsampling"] = subsampling
        f["stride"] = stride
        close(f)
    end
    return T, τ, gradients
end

# load outer loop data
using HDF5
# optical data
# path = "SARA/NatCom2020/outer/data/"
# path = "NatCom2020/"
# file = "Bi2O3_19F44_01_outer_loop_data.h5" # without input noise
# file = "Bi2O3_19F44_01_outer_loop_input_noise_data.h5"
# file = "Bi2O3_19F44_01_outer_loop_input_noise_constant_offset.h5"
# xrd data
path = "SARA/NatCom2020/outer/data/"
# path = "NatCom2020/"
# file = "Bi2O3_19F44_01_outer_xrd_gradients_input_noise.h5"
file = "Bi2O3_19F44_01_outer_optical_gradients_l_2_input_noise_50.h5"

f = h5open(path * file, "r")
temperatures = read(f, "temperatures")
dwelltimes = read(f, "dwelltimes")
gradients = read(f, "gradients")
var_gradients = read(f, "var(gradients)")
close(f)

# normalize global data
y_normalization = std(gradients)
gradients ./= y_normalization
var_gradients ./= y_normalization^2 # noise variance

x = tuple.(temperatures, dwelltimes)
y = gradients
σ² = var_gradients

T, τ, gradients = optical_gradient_map(x, y, σ²)

################################################################################
# plot it
# using HDF5
# # f = h5open("SARA/NatCom2020/results/outer_optical_gradient_map.h5")
# f = h5open("outer_optical_gradient_map.h5")
#
# T = read(f, "temperatures")
# τ = read(f, "dwelltimes")
# gradients = read(f, "gradients")
# A = gradients[:, :, 1, 1, 1]
#
# using Plots
# pyplot()
# heatmap(T, τ, A')
# gui()
# close(f)

################################################################################

# calculate single map
# subsample(x, k::Int) = @view x[end+1-k:end, :]
# n = 8
# x = subsample(x, n)
# y = subsample(y, n)
# σ² = subsample(σ², n)

# using Kernel
# # from marginal likelihood opt:
# T_l = 26.629076190570203
# τ_l = .3 # 15.924969663921164
# NN_const = 1111.2566343353083
# # T_l, τ_l = 10., .2
# # kernel = oSARA([T_l, τ_l])
# kernel = Kernel.NN(NN_const)
# scaling = [T_l, τ_l]
# scale(k) = (x, y) -> k(x./scaling, y./scaling)
# kernel = scale(kernel)
# G = Gaussian(kernel)
# x = vec(x)
# y = vec(y)
# # μ_yk = mean(yk)
# # yk = yk .- μ_yk
# σ² = vec(σ²)
# @time C = conditional(G, x, Gaussian(y, Diagonal(σ²)), tol = 1e-6)
# # C = C + μ_yk
# n = 64
# T = range(200, 1200, length = n) # plotting ranges
# τ = range(2.5, 4, length = n + 1)
# gradients = mean(C).(tuple.(T, τ'))
# heatmap(T, τ, gradients')

# using Plots
# pyplot()
# n = 64
# T = range(200, 1200, length = n)
# τ = range(2.5, 4, length = n)
# heatmap(T, τ, (x, y)->mean(C)((x, y)))
# gui()

# marginal likelihood optimization of nn kernel
# function nn_kernel(θ)
#     l1, l2, σ = θ
#     scaling = [l1, l2]
#     scale(k) = (x, y) -> k(x./scaling, y./scaling)
#     return scale(Kernel.NN(σ))
# end
# using GaussianDistributions: optimize
# θ = [40, 0.3, 1000]
#
# subsample_spec(x, k::Int) = @view x[:, 1:k]
# d = 8
# xk, yk, σ²k = subsample_spec(x, d), subsample_spec(y, d), subsample_spec(σ², d)
# xk, yk, σ²k = vec.((xk, yk, σ²k))
# optimize(nn_kernel, θ, xk, Gaussian(yk, σ²k))

ext: jl


text:
using Conda

Conda.add("ase", channel="conda-forge")
Conda.add("rdkit", channel="conda-forge")
Conda.add("pymatgen", channel="conda-forge")


text:
using Base.Test
using ReferenceTests

# check for ambiguities
refambs = detect_ambiguities(Base, Core)
using SmartGameFormat
ambs = detect_ambiguities(Base, Core)
@test length(setdiff(ambs, refambs)) == 0

tests = [
    "tst_lexer.jl",
    "tst_types.jl",
    "tst_parser.jl",
]

for t in tests
    @testset "$t" begin
        include(t)
    end
end


text:
abstract type Message end

abstract type Report <: Message end

struct Abort <: Message end

struct DumpState <: Message end

struct Initialize{R <: AbstractPenaltyParameter} <: Message
    create_subproblem::Function
    create_subproblem_args::Tuple
    create_subproblem_kwargs::NamedTuple
    r::Type{R}
    scenarios::Set{ScenarioID}
    scenario_tree::ScenarioTree
    warm_start::Bool
    subproblem_callbacks::Vector{SubproblemCallback}
end

struct InitializeLowerBound <: Message
    create_subproblem::Function
    create_subproblem_args::Tuple
    create_subproblem_kwargs::NamedTuple
    scenarios::Set{ScenarioID}
    scenario_tree::ScenarioTree
    warm_start::Bool
end

struct PenaltyInfo <: Message
    scen::ScenarioID
    penalty::Union{Float64,Dict{VariableID,Float64}}
end

struct Ping <: Message end

struct ReportBranch <: Report
    scen::ScenarioID
    sts::MOI.TerminationStatusCode
    obj::Float64
    time::Float64
    vals::Dict{VariableID,Float64}
end

struct ReportLeaf <: Report
    scen::ScenarioID
    vals::Dict{VariableID,Float64}
end

struct ReportLowerBound <: Report
    scen::ScenarioID
    sts::MOI.TerminationStatusCode
    obj::Float64
    time::Float64
end

struct ShutDown <: Message end

struct Solve <: Message
    scen::ScenarioID
    w_vals::Dict{VariableID,Float64}
    xhat_vals::Dict{VariableID,Float64}
    niter::Int
end

struct SolveLowerBound <: Message
    scen::ScenarioID
    w_vals::Dict{VariableID,Float64}
end

struct SubproblemAction <: Message
    scen::ScenarioID
    action::Function
    args::Tuple
    kwargs::NamedTuple
end

struct VariableMap <: Message
    scen::ScenarioID
    var_info::Dict{VariableID,VariableInfo} # VariableInfo definition in subproblem.jl
end

struct WorkerState{T} <: Message
    id::Int
    state::T
end


text:
module _Machin

using Brainstorm._Math._Series._Arctan: arctanhsum

# http://numbers.computation.free.fr/Constants/Log2/log2.html
ln2_N2a(digits::Integer) = arctanhsum(digits, [5,  7], [2, 2])
ln2_N2b(digits::Integer) = arctanhsum(digits, [6, 99], [4, 2])
ln2_N2c(digits::Integer) = arctanhsum(digits, [7, 17], [4, 2])

ln2_N3a(digits::Integer) = arctanhsum(digits, [11,  111, 19601], [ 8, -4, -2])
ln2_N3b(digits::Integer) = arctanhsum(digits, [17,   79,  1351], [10,  8,  4])
ln2_N3c(digits::Integer) = arctanhsum(digits, [23,   65,   485], [14,  6, -4])
ln2_N3d(digits::Integer) = arctanhsum(digits, [26, 4801,  8749], [18, -2,  8])
ln2_N3e(digits::Integer) = arctanhsum(digits, [31,   49,   161], [14, 10,  6])

ln2_N4a(digits::Integer) = arctanhsum(digits, Int64[ 52, 4801, 8749, 70226], [ 36, -2,   8,  18])
ln2_N4b(digits::Integer) = arctanhsum(digits,      [127,  449, 4801,  8749], [ 72, 54,  34, -10])
ln2_N4c(digits::Integer) = arctanhsum(digits,      [251,  449, 4801,  8749], [144, 54, -38,  62])

ln2_N5a(digits::Integer) = arctanhsum(digits, Int64[575, 3361, 8749, 13121, 56251], [342, 198, 222, 160, 106])

end


text:
struct DecisionIndex <: MOI.AbstractConstraintAttribute end
MOI.is_set_by_optimize(::DecisionIndex) = true

struct DecisionCoefficientChange{T} <: MOI.AbstractFunctionModification
    decision::MOI.VariableIndex
    new_coefficient::T
end

struct DecisionMultirowChange{T} <: MOI.AbstractFunctionModification
    decision::MOI.VariableIndex
    new_coefficients::Vector{Tuple{Int64, T}}
end

struct DecisionStateChange <: MOI.AbstractFunctionModification
    index::Int
    new_state::DecisionState
end

struct KnownValuesChange <: MOI.AbstractFunctionModification end

const DecisionModification = Union{DecisionCoefficientChange, DecisionStateChange, KnownValuesChange}
const VectorDecisionModification = Union{DecisionMultirowChange, DecisionStateChange, KnownValuesChange}


text:
module Utils

compile_enabled() = Base.JLOptions().compile_enabled == 1
# allow `== 2` (`--compile=all`) as well?

end  # module


text:
Base.print(io::IO, loss::SupervisedLoss, args...) = print(io, typeof(loss).name.name, args...)
Base.print(io::IO, loss::L1DistLoss, args...) = print(io, "L1DistLoss", args...)
Base.print(io::IO, loss::L2DistLoss, args...) = print(io, "L2DistLoss", args...)
Base.print(io::IO, loss::LPDistLoss{P}, args...) where {P} = print(io, typeof(loss).name.name, " with P = $(P)", args...)
Base.print(io::IO, loss::L1EpsilonInsLoss, args...) = print(io, typeof(loss).name.name, " with \$\\epsilon\$ = $(loss.ε)", args...)
Base.print(io::IO, loss::L2EpsilonInsLoss, args...) = print(io, typeof(loss).name.name, " with \$\\epsilon\$ = $(loss.ε)", args...)
Base.print(io::IO, loss::QuantileLoss, args...) = print(io, typeof(loss).name.name, " with \$\\tau\$ = $(loss.τ)", args...)
Base.print(io::IO, loss::SmoothedL1HingeLoss, args...) = print(io, typeof(loss).name.name, " with \$\\gamma\$ = $(loss.gamma)", args...)
Base.print(io::IO, loss::HuberLoss, args...) = print(io, typeof(loss).name.name, " with \$\\alpha\$ = $(loss.d)", args...)
Base.print(io::IO, loss::DWDMarginLoss, args...) = print(io, typeof(loss).name.name, " with q = $(loss.q)", args...)
Base.print(io::IO, loss::PeriodicLoss, args...) = print(io, typeof(loss).name.name, " with c = $(round(2π / loss.k, digits=1))", args...)
Base.print(io::IO, loss::ScaledLoss{T,K}, args...) where {T,K} = print(io, "$(K) * ($(loss.loss))", args...)

_round(num) = round(num) == round(num, digits=1) ? round(Int, num) : round(num, digits=1)
function _relation(num)
    if num <= 0
        "negative only"
    elseif num >= 1
        "positive only"
    elseif num < 0.5
        "1:$(_round((1-num)/num)) weighted"
    else
        "$(_round(num/(1-num))):1 weighted"
    end
end
Base.print(io::IO, loss::WeightedMarginLoss{T,W}, args...) where {T,W} = print(io, "$(_relation(W)) $(loss.loss)", args...)


text:
module TestBuildFields
using Compat
using Compat.Test
using CategoricalArrays
using CategoricalArrays: DefaultRefType

@testset "buildindex(), buildinvindex(), buildorder() for b a c" begin
    index = ["b", "a", "c"]

    invindex = Dict(
        "b" => DefaultRefType(1),
        "a" => DefaultRefType(2),
        "c" => DefaultRefType(3),
    )

    order = [
        DefaultRefType(2),
        DefaultRefType(1),
        DefaultRefType(3),
    ]

    pool = CategoricalPool(index, invindex)

    levels = ["c", "a", "b"]

    built_index = CategoricalArrays.buildindex(invindex)
    @test isa(index, Vector)
    @test built_index == index

    built_invindex = CategoricalArrays.buildinvindex(index)
    @test isa(invindex, Dict)
    @test built_invindex == invindex

    neworder = [
        DefaultRefType(3),
        DefaultRefType(2),
        DefaultRefType(1),
    ]

    built_order = CategoricalArrays.buildorder(pool.invindex, levels)
    @test isa(order, Vector{DefaultRefType})
    @test built_order == neworder
end

end


text:
# MIT License
#
# Copyright (c) 2020 Joshua E Gentry

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.
using Dates

# ========================================================================================+
"""
    mutable struct GitHubFile
        filename::String
        timestamp::Float64
        localsha1::String
        remotesha1::String
        url::String
    end

This object holds details about a file that is downloaded from GitHub.  It stores the SHA1
(as calculated by GitHub) for both the server version and the local version.  Along with
with the timestamp of when the local file was last modified.

__Values:__

| Name        | Description
| :---------- | :----------
| filename    | The full name (without the path) of the training data file.
| timestamp   | The last time we modifed the file.
| localsha1   | The local SHA1 for the file.
| remotesha1  | The remote SHA1 for the same file.
| url         | The URL that can be used to download the file.
"""
mutable struct GitHubFile
    filename::String
    timestamp::Float64
    localsha1::String
    remotesha1::String
    url::String
    GitHubFile() = new("", 0.0, "", "", "")
    GitHubFile(filename, timestamp, localsha, remotesha, url) = new(filename, timestamp, localsha, remotesha, url)
end

# =================================================================================================
"""
    parse_github_file(
        record
    )::Union{DfFGitHubFileile, Nothing}

Construct a GitHubFile object from a JSON record.  If the JSON is invalid returns
`nothing`.

__Arguments:__

| T | Name   | Default | Description
|:--| :----- | :------ | :----------
| R | record |         | The JSON object to pull the values from.
"""
function parse_github_file(
            record
        )::Union{GitHubFile, Nothing}

    if isa(record, AbstractDict) == false
        return nothing
    end

    # -------------------------------------------------------------------------------------
    # Extract the values.
    local filename    = get_json_safe(AbstractString, record, "filename")
    local timestamp   = get_json_safe(Float64, record, "timestamp")
    local localsha1   = get_json_safe(AbstractString, record, "localsha1")
    local remotesha1  = get_json_safe(AbstractString, record, "remotesha1")
    local url         = get_json_safe(AbstractString, record, "url")

    # -------------------------------------------------------------------------------------
    # Make sure the values exist.
    filename === nothing && return nothing
    timestamp === nothing && return nothing
    localsha1 === nothing && return nothing
    remotesha1 === nothing && return nothing
    url === nothing && return nothing

    # -------------------------------------------------------------------------------------
    # Create the GitHubFile object.
    return GitHubFile(filename, timestamp, localsha1, remotesha1, url)
end

# =================================================================================================
"""
file_out_of_date(
        file::GitHubFile,
        path::AbstractString
    )::Bool

Check if the file is out of date with the server.

__Arguments:__

| T | Name | Default | Description
|:--| :--- | :------ | :----------
| R | file |         | The file to check if it's out of date.
| R | path |         | The directory to look for the file in.

__Exceptions:__

| Name         | Reason
| :----------- | :-----
| FileError    | There was an issue reading from the file.

__Details:__

Primarily this method checks if the local sha1 matches the server's sha1.  If they match no
update is needed.  However if the local training data file is missing, then we need to
update it.  Also if the timestamp of the file is different than the timestamp we have we
need to re-calculate the sha1 because it might have changed.

This method will update the localsh1 and timestamp properties if they are out of sync with
the local file system.
"""
function file_out_of_date(
            file::GitHubFile,
            path::AbstractString
        )::Bool
    local fullpath = joinpath(path, file.filename)

    # ---------------------------------------------------------------------------------------------
    # If the file doesn't exist, then it needs to be updated.
    if isfile(fullpath) == false
        return true
    end

    # ---------------------------------------------------------------------------------------------
    # If the timestamp doesn't match our timestamp, recalculate the sha1.
    local timestamp = mtime(fullpath)
    if timestamp != file.timestamp
        file.localsha1 = calc_sha1(fullpath)
        file.timestamp = timestamp
    end

    # ---------------------------------------------------------------------------------------------
    # Compare our sha1 with the server's.
    return file.localsha1 != file.remotesha1
end


text:
include("ModAlgAss/ModAlgAss.jl")
include("ModAlgAss/MeatAxe.jl")


text:
module RecSys

println("You have decided to learn Recommender Systems using Julia. Please refer to the document Julia_RecSys.jl")

end


text:
module SerialCOBS
    using LibSerialPort
    using StaticArrays
    using ProtoBuf

    include("proto_helper.jl")
    include("Arduino.jl")
    include("COBS.jl")

    export Arduino, messageBytesAvailable
    export message, recieve!
    export encode, decode

    function message(ard::Arduino, msg::MSG_TYPE) where {MSG_TYPE <: ProtoBuf.ProtoType}
    # function message(ard::Arduino, msg::AbstractVector{UInt8})
        iob = IOBuffer()
        msg_size = writeproto(iob, msg)
        msg = @view iob.data[1:msg_size]

        encoded_msg = encode(ard, msg)

        status = write(ard, encoded_msg)

        return status
    end

    # function recieve(ard::Arduino{MSG_TYPE}, msg::MSG_TYPE) where {MSG_TYPE <: ProtoBuf.ProtoType}
    function recieve!(ard::Arduino, msg::MSG_TYPE) where {MSG_TYPE <: ProtoBuf.ProtoType}
        encoded_msg = (SerialCOBS.retrieve_encoded_msg(ard))
        if encoded_msg !== nothing
            decoded_msg = decode(ard, encoded_msg)
            readproto(IOBuffer(decoded_msg), msg)

            return true
        end

        return false
    end
end  # module SerialCOBS


text:
module PipsNlpInterface

#TODO: update pips-nlp interface to work on optigraphs in addition to optinodes for subproblems

using SparseArrays
using LinearAlgebra
using DataStructures
using PipsNLP
import MPI

import JuMP
import MathOptInterface
const MOI = MathOptInterface
using Plasmo

include("pipsnlp_c_interface.jl")
using .PipsNlpSolver

include("pips_utils.jl")
include("nlp_evaluator.jl")

export pipsnlp_solve

function pipsnlp_solve(graph::OptiGraph) #Assume graph variables and constraints are first stage
    #Check structure: 3 possibilities
    #1.) optigraph without subgraphs (only link constraints, no first stage -> second stage constraints)
    #2.) optigraph with single subgraph (first stage and linkconstraints between optinodes)
    #3.) Not yet supported: optigraph with multiple subgraphs (first stage and linkconstraints between subgraphs)

    comm = MPI.COMM_WORLD
    if MPI.Comm_rank(comm) == 0
        println("Building Model for PIPS-NLP")
    end

    #NOTE: either the interface figures out the linking structure, or the user provided it (e.g. using distribute_optigraph)
    if haskey(graph.ext,:user_pips_data)
        worker_data = graph.ext[:user_pips_data]
    else
        worker_data = PipsNLP._setup_pips_nlp_data!(graph)
    end

    first_stage = worker_data.first_stage
    submodels = setdiff(all_nodes(graph),[first_stage])
    n_sub_models = length(submodels)
    model_list = [first_stage; submodels]

    #unpack data
    first_stage_data = PipsNLP._get_pips_data(first_stage)
    nlinkeq = worker_data.n_linkeq_cons
    nlinkineq = worker_data.n_linkineq_cons
    ineqlink_lb = worker_data.link_ineq_lower
    ineqlink_ub = worker_data.link_ineq_upper
    eqlink_lb = worker_data.link_eq_lower
    eqlink_ub = worker_data.link_eq_upper

    #initialize values
    function str_init_x0(nodeid, x0)
        node = model_list[nodeid+1]
        jump_initval = copy(JuMP.start_value.(JuMP.all_variables(node)))
        local_initval = Vector{Float64}(undef,JuMP.num_variables(node))

        nothing_values = isa.(jump_initval,Nothing)
        float_values = .!(nothing_values)
        local_initval[nothing_values] .= 0  #set to 0 by default
        local_initval[float_values] .= jump_initval[float_values]
        original_copy(local_initval,x0)
    end

    #get problem info
    function str_prob_info(nodeid,flag,mode,col_lb,col_ub,row_lb,row_ub)
        if flag != 1
            node = model_list[nodeid+1]
            local_data = PipsNLP._get_pips_data(node)
            #do a warm start on the node
            if !(local_data.loaded)
                local_data.loaded = true
                if (nodeid > 0)
                    if haskey(node.ext, :warmStart)
    	                if node.ext[:warmStart] == true
              	      	   optimize!(node)
                   	   	end
                    end
            	end

                nlp_lb, nlp_ub = constraintbounds(node.ext[:constraint_data])       #This is every constraint in the model
         		local_data.local_m  = length(nlp_lb)          #number of local constraints (rows)
    			newRowId = Array{Int}(undef,local_data.local_m)
    			eqId = 1
    			ineqId = 1

                #go through all constraints, check if they are equality of inequality
                for c in 1:local_data.local_m
                    if nlp_lb[c] == nlp_ub[c]
                        push!(local_data.eq_idx, c)
                        newRowId[c] = eqId
                        eqId +=  1
                    else
        				push!(local_data.ineq_idx, c)
        				newRowId[c] = ineqId
        				ineqId += 1
                    end
         		end

                #Local node data
                local_data.m  = local_data.local_m + local_data.num_eqconnect + local_data.num_ineqconnect
                local_data.n = JuMP.num_variables(node)
                local_data.x_sol = zeros(Float64,local_data.n)
                local_data.firstJeqmat = sparse(local_data.firstIeq, local_data.firstJeq, local_data.firstVeq, local_data.num_eqconnect, first_stage_data.n)
                local_data.secondJeqmat = sparse(local_data.secondIeq, local_data.secondJeq, local_data.secondVeq, local_data.num_eqconnect, local_data.n)
                local_data.firstJineqmat = sparse(local_data.firstIineq, local_data.firstJineq, local_data.firstVineq, local_data.num_ineqconnect, first_stage_data.n)
                local_data.secondJineqmat = sparse(local_data.secondIineq, local_data.secondJineq, local_data.secondVineq, local_data.num_ineqconnect, local_data.n)

                if node.nlp_data == nothing
                    JuMP._init_NLP(node)
                end
                local_data.d = JuMP.NLPEvaluator(node)
                MOI.initialize(local_data.d, [:Grad,:Jac, :Hess])
                jac_structure = pips_jacobian_structure(local_data.d)
                Ijac = [jac_structure[i][1] for i = 1:length(jac_structure)]
                Jjac = [jac_structure[j][2] for j = 1:length(jac_structure)]
                Ijaceq = Int[]
                Jjaceq = Int[]
                Ijacineq = Int[]
                Jjacineq = Int[]
                jac_eq_index = Int[]
                jac_ineq_index = Int[]
                for i in 1:length(Ijac)
                    c = Ijac[i]
                    if nlp_lb[c] == nlp_ub[c]
                        modifiedrow = newRowId[c]
                        push!(Ijaceq, modifiedrow)
                        push!(Jjaceq, Jjac[i])
                        push!(jac_eq_index, i)
                    else
                        modifiedrow = newRowId[c]
                        push!(Ijacineq, modifiedrow)
                        push!(Jjacineq, Jjac[i])
                        push!(jac_ineq_index,i)
                    end
                end

        		node.ext[:Ijaceq] = Ijaceq
         		node.ext[:Jjaceq] = Jjaceq
                node.ext[:Ijacineq] = Ijacineq
         		node.ext[:Jjacineq] = Jjacineq
         		node.ext[:jac_eq_index] = jac_eq_index
         		node.ext[:jac_ineq_index] = jac_ineq_index
                hess_structure = pips_hessian_lagrangian_structure(local_data.d)
                Ihess = [hess_structure[i][1] for i = 1:length(hess_structure)]
                Jhess = [hess_structure[j][2] for j = 1:length(hess_structure)]

         		Hmap = Bool[]
         		node_Hrows = Int[]
         		node_Hcols = Int[]
         		for i in 1:length(Ihess)
            	    if Jhess[i] <= Ihess[i]
               	        push!(node_Hrows, Ihess[i])
                   		push!(node_Hcols, Jhess[i])
                   		push!(Hmap, true)
            	    else
                        push!(Hmap, false)
                    end
                end
         		val = ones(Float64,length(node_Hrows))
         		mat = sparse(node_Hrows,node_Hcols,val, local_data.n, local_data.n)
         		node.ext[:Hrows] = node_Hrows
         		node.ext[:Hcols] = node_Hcols
         		node.ext[:Hmap] = Hmap
         		local_hessnnz = length(mat.rowval)
         		local_data.local_unsym_hessnnz = length(Ihess)
         		local_data.hessnnz = local_hessnnz
		    end

	 	    if mode == :Values
                nlp_lb, nlp_ub = constraintbounds(node.ext[:constraint_data])
    			eq_lb=Float64[]
    			eq_ub=Float64[]
    			ineq_lb=Float64[]
    			ineq_ub=Float64[]
    			for i in 1: length(nlp_lb)
    			    if nlp_lb[i] == nlp_ub[i]
    			       push!(eq_lb, nlp_lb[i])
    			       push!(eq_ub, nlp_ub[i])
    			    else
    			       push!(ineq_lb, nlp_lb[i])
    			       push!(ineq_ub, nlp_ub[i])
    			    end
    			end

    			if nodeid !=  0 #0 is the master
    			   eq_lb = [eq_lb; local_data.eqconnect_lb]
    			   eq_ub = [eq_ub; local_data.eqconnect_ub]
    			   ineq_lb = [ineq_lb; local_data.ineqconnect_lb]
    			   ineq_ub = [ineq_ub; local_data.ineqconnect_ub]
    			end

    			original_copy([eq_lb;ineq_lb], row_lb)
    			original_copy([eq_ub;ineq_ub], row_ub)
                colLower = variablelowerbounds(node)
                colUpper = variableupperbounds(node)
                original_copy(colLower, col_lb)
    			original_copy(colUpper, col_ub)
            end
		    return (local_data.n,local_data.m)
		else
		    if mode == :Values
                original_copy([eqlink_lb; ineqlink_lb], row_lb)
                original_copy([eqlink_ub; ineqlink_ub], row_ub)
            end
            return (0, nlinkeq + nlinkineq)
		end
    end

    #evaluate objective
    #x0 is first stage variable values, x1 is local values
    function str_eval_f(nodeid,x0,x1)
    	node = model_list[nodeid+1] #Julia doesn't start index at 0
        local_data = PipsNLP._get_pips_data(node)
        local_d = PipsNLP._get_pips_data(node).d
        if nodeid ==  0
            local_x = x0
        else
            local_x = x1
        end
        #check objective sign
        local_scl = (JuMP.objective_sense(node) == MOI.MAX_SENSE) ? -1.0 : 1.0
        f = local_scl*pips_eval_objective(local_d,local_x)
        return f
    end

    #evaluate constraints
    function str_eval_g(nodeid,x0,x1,new_eq_g, new_inq_g)
        node = model_list[nodeid+1]
        local_data = PipsNLP._get_pips_data(node)
        local_d = PipsNLP._get_pips_data(node).d
        if nodeid ==  0
            local_x = x0
        else
            local_x = x1
        end
        local_g = Array{Float64}(undef,local_data.local_m)
        pips_eval_constraint(local_d, local_g, local_x)

        #evaluates local constraints and connecting constraints with first stage
        new_eq_g[1:end] = [local_g[local_data.eq_idx]; local_data.firstJeqmat*x0+local_data.secondJeqmat*x1]
        new_inq_g[1:end] = [local_g[local_data.ineq_idx]; local_data.firstJineqmat*x0+local_data.secondJineqmat*x1]
	    return Int32(1)
    end

    #get solution
    function str_write_solution(id::Integer, x::Vector{Float64}, y_eq::Vector{Float64}, y_ieq::Vector{Float64})
        node = model_list[id+1]
        local_data = PipsNLP._get_pips_data(node)
        local_data.x_sol = copy(x)

        #TODO: Grab dual values
        rank = MPI.Comm_rank(comm)
        local_data.coreid = rank
    end

    #evaluate gradient
    function str_eval_grad_f(rowid,colid,x0,x1,new_grad_f)
        node = model_list[rowid+1]
        if rowid == colid
            local_data = PipsNLP._get_pips_data(node)
            local_d = PipsNLP._get_pips_data(node).d
            if colid ==  0
                local_x = x0
            else
                local_x = x1
            end
            local_grad_f = Array{Float64}(undef,local_data.n)
            pips_eval_objective_gradient(local_d, local_grad_f, local_x)
            #eval_grad_f(local_d, local_grad_f, local_x)
            #local_scl = (node.objSense == :Min) ? 1.0 : -1.0
            local_scl = (JuMP.objective_sense(node) == MOI.MAX_SENSE) ? -1.0 : 1.0
            #scale!(local_grad_f,local_scl)
            rmul!(local_grad_f,local_scl)
            original_copy(local_grad_f, new_grad_f)
        elseif colid == 0
            new_grad_f[1:end] .= 0
        else
            @assert(false)
        end
        return Int32(1)
    end

    #copy row/column indices to C pointer
    function array_copy(src,dest)
        @assert(length(src)==length(dest))
        for i in 1:length(src)
            dest[i] = src[i]-1
        end
    end

    #copy array elements
    function original_copy(src,dest)
        @assert(length(src)==length(dest))
        for i in 1:length(src)
            dest[i] = src[i]
        end
    end

    #evaluate jacobian
    function str_eval_jac_g(rowid,colid,flag, x0,x1,mode,e_rowidx,e_colptr,e_values,i_rowidx,i_colptr,i_values)
        if flag != 1 #populate parent child structure
            node = model_list[rowid+1]
            local_data = PipsNLP._get_pips_data(node)
            local_m_eq = length(local_data.eq_idx)
            local_m_ineq = length(local_data.ineq_idx)
            if mode == :Structure
                if (rowid == colid)
                    Ieq=[node.ext[:Ijaceq];local_data.secondIeq .+ local_m_eq]
                    Jeq=[node.ext[:Jjaceq];local_data.secondJeq]
                    Veq= ones(Float64, length(Ieq))
                    Iineq=[node.ext[:Ijacineq];local_data.secondIineq .+ local_m_ineq]
                    Jineq=[node.ext[:Jjacineq];local_data.secondJineq]
                    Vineq=ones(Float64, length(Iineq))
                    eqmat = sparse(Ieq, Jeq, Veq, local_m_eq + local_data.num_eqconnect, local_data.n)
                    ineqmat = sparse(Iineq, Jineq, Vineq, local_m_ineq + local_data.num_ineqconnect, local_data.n)
                else
                    eqmat = sparse(local_m_eq .+ local_data.firstIeq, local_data.firstJeq, local_data.firstVeq, local_m_eq .+ local_data.num_eqconnect, first_stage_data.n)
                    ineqmat = sparse(local_m_ineq .+ local_data.firstIineq, local_data.firstJineq, local_data.firstVineq, local_m_ineq .+ local_data.num_ineqconnect, first_stage_data.n)
                end
                return(length(eqmat.rowval), length(ineqmat.rowval))
            else #mode = :Values
                if rowid == colid    #evaluate Jacobian for block
                    if colid ==  0   #if it's the root node
                        local_x = x0
                    else
                        local_x = x1
                    end
                    local_values = Array{Float64}(undef,length(node.ext[:Ijaceq])+length(node.ext[:Ijacineq]))
                    #eval_jac_g(local_data.d, local_values, local_x)
                    pips_eval_constraint_jacobian(local_data.d, local_values, local_x)
                    jac_eq_index = node.ext[:jac_eq_index]
                    jac_ineq_index = node.ext[:jac_ineq_index]
                    Ieq=[node.ext[:Ijaceq];local_data.secondIeq .+ local_m_eq]
                    Jeq=[node.ext[:Jjaceq];local_data.secondJeq]
                    Veq=[local_values[jac_eq_index];local_data.secondVeq]
                    Iineq=[node.ext[:Ijacineq];local_data.secondIineq .+ local_m_ineq]
                    Jineq=[node.ext[:Jjacineq];local_data.secondJineq]
                    Vineq=[local_values[jac_ineq_index];local_data.secondVineq]
                    eqmat = _sparse_keep_zero(Ieq, Jeq, Veq, local_m_eq + local_data.num_eqconnect, local_data.n)
                    ineqmat = _sparse_keep_zero(Iineq, Jineq, Vineq, local_m_ineq + local_data.num_ineqconnect, local_data.n)
                else #evaluate
                    eqmat = _sparse_keep_zero(local_m_eq .+ local_data.firstIeq, local_data.firstJeq, local_data.firstVeq, local_m_eq .+ local_data.num_eqconnect, first_stage_data.n)
                    ineqmat = _sparse_keep_zero(local_m_ineq .+ local_data.firstIineq, local_data.firstJineq, local_data.firstVineq, local_m_ineq .+ local_data.num_ineqconnect, first_stage_data.n)
                end
                if length(eqmat.nzval) > 0
                    array_copy(eqmat.rowval,e_rowidx)
                    array_copy(eqmat.colptr,e_colptr)
                    original_copy(eqmat.nzval,e_values)
                end
                if length(ineqmat.nzval) > 0
                    array_copy(ineqmat.rowval,i_rowidx)
                    array_copy(ineqmat.colptr,i_colptr)
                    original_copy(ineqmat.nzval, i_values)
                end
            end
        else #populate linkconstraint structure
            node = model_list[rowid+1]
            local_data = PipsNLP._get_pips_data(node)
            linkIeq = local_data.linkIeq
            linkJeq = local_data.linkJeq
            linkVeq = local_data.linkVeq
            linkIineq = local_data.linkIineq
            linkJineq = local_data.linkJineq
            linkVineq = local_data.linkVineq

            if mode == :Structure
                return(length(linkVeq), length(linkVineq))
            else
                eqmat = sparse(linkIeq, linkJeq, linkVeq, nlinkeq, local_data.n)
                ineqmat = sparse(linkIineq, linkJineq, linkVineq, nlinkineq, local_data.n)
                if length(eqmat.nzval) > 0
                    array_copy(eqmat.rowval,e_rowidx)
                    array_copy(eqmat.colptr,e_colptr)
                    original_copy(eqmat.nzval,e_values)
                end
                if length(ineqmat.nzval) > 0
                    array_copy(ineqmat.rowval,i_rowidx)
                    array_copy(ineqmat.colptr,i_colptr)
                    original_copy(ineqmat.nzval, i_values)
                end
            end
        end
        return Int32(1)
    end

    #evaluate hessian
    function str_eval_h(rowid,colid,x0,x1,obj_factor,lambda,mode,rowidx,colptr,values)
        node = model_list[colid+1]
        local_data = PipsNLP._get_pips_data(node)
        if mode == :Structure
            if rowid == colid
                node_Hrows  = node.ext[:Hrows]
                node_Hcols = node.ext[:Hcols]
                node_val = ones(Float64,length(node_Hrows))
                mat = _sparse_keep_zero(node_Hrows,node_Hcols,node_val, local_data.n, local_data.n)
                return length(mat.rowval)
            elseif colid == 0
                node_Hrows  = node.ext[:Hrows]
                node_Hcols = node.ext[:Hcols]
                node_Hrows,node_Hcols = _exchange(node_Hrows,node_Hcols)
                node_val = ones(Float64,length(node_Hrows))
                mat = _sparse_keep_zero(node_Hrows,node_Hcols,node_val, local_data.n, local_data.n)
                return length(mat.rowval)
            else
                return 0
            end
        else #fill in values
            if rowid == colid
                node_Hmap = node.ext[:Hmap]
                node_Hrows  = node.ext[:Hrows]
                node_Hcols = node.ext[:Hcols]
                if colid ==  0
                    local_x = x0
                else
                    local_x = x1
                end
                local_unsym_values = Array{Float64}(undef,local_data.local_unsym_hessnnz)
                node_val = ones(Float64,length(node_Hrows))
                local_scl = (JuMP.objective_sense(node) == MOI.MAX_SENSE) ? -1.0 : 1.0
                local_m_eq = length(local_data.eq_idx)
                local_m_ineq = length(local_data.ineq_idx)
                local_lambda = zeros(Float64, local_data.local_m)
                for i in 1:local_m_eq
                    local_lambda[local_data.eq_idx[i]] = lambda[i]
                end
                for i in 1:local_m_ineq
                    local_lambda[local_data.ineq_idx[i]] = lambda[i+local_m_eq]
                end

                #eval_hesslag(local_data.d, local_unsym_values, local_x, obj_factor*local_scl, local_lambda)
                pips_eval_hessian_lagrangian(local_data.d, local_unsym_values, local_x, obj_factor*local_scl, local_lambda)
                local_sym_index=1
                for i in 1:local_data.local_unsym_hessnnz
                    if node_Hmap[i]
                        node_val[local_sym_index] = local_unsym_values[i]
                        local_sym_index +=1
                    end
                end
                node_Hrows,node_Hcols = _exchange(node_Hrows,node_Hcols)
                mat = _sparse_keep_zero(node_Hrows,node_Hcols,node_val, local_data.n, local_data.n)
                array_copy(mat.rowval,rowidx)
                array_copy(mat.colptr,colptr)
                original_copy(mat.nzval,values)
            elseif colid ==	0
                node_Hrows  = node.ext[:Hrows]
                node_Hcols = node.ext[:Hcols]
                node_val = zeros(Float64,length(node_Hrows))
                node_Hrows,node_Hcols = _exchange(node_Hrows,node_Hcols)
                mat = _sparse_keep_zero(node_Hrows,node_Hcols,node_val, local_data.n, local_data.n)
                array_copy(mat.rowval,rowidx)
                array_copy(mat.colptr,colptr)
                original_copy(mat.nzval,values)
            else
            end
        end
        return Int32(1)
    end

    #initalize MPI if not already
    if !MPI.Initialized()
        MPI.Init()
    end

    # comm = MPI.COMM_WORLD
    if (MPI.Comm_rank(comm) == 0)
        t1 = time()
    end

    root = 0
    rank = MPI.Comm_rank(comm)

    ###################################
    #PIPS-NLP FUNCTIONS
    ###################################
    #Create PipsModel (The PIPS interface model) and pass all the functions it requires
    model = PipsModel(:Min, 0,
    n_sub_models,
    str_init_x0,
    str_prob_info,
    str_eval_f,
    str_eval_g,
    str_eval_grad_f,
    str_eval_jac_g,
    str_eval_h,
    str_write_solution)

    #create the C++ problem struct
    prob = createProblemStruct(comm, model, true)
    if rank == 0
        println("Created PIPS-NLP Problem Struct")
    end

    #probably don't need this barrier
    MPI.Barrier(comm)

    #solve the problem struct.  MPI will do communication in the solver.
    ret = solveProblemStruct(prob)
    ###################################

    if rank == 0
        println("Solved PIPS-NLP Problem Struct")
    end

    if rank == 0
        println("Timing Results:")
        println("init_x0  ",prob.t_jl_init_x0)
        println("str_init_x0  ", prob.t_jl_str_prob_info)
        println("eval_f  ", prob.t_jl_eval_f)
        println("eval_g0  ",prob.t_jl_eval_g)
        println("eval_grad_f  ", prob.t_jl_eval_grad_f)
        println("eval_jac  ", prob.t_jl_eval_jac_g)
        println("str_eval_jac  ", prob.t_jl_str_eval_jac_g)
        println("eval_h  ",  prob.t_jl_eval_h)
        println("str_eval_h ",  prob.t_jl_str_eval_h)
        println("eval_write_solution  ",  prob.t_jl_write_solution)
        println("PIPS-NLP time:   ",  time() - t1, " (s)")
    end

    status_dict = Dict(
    0 => MOI.LOCALLY_SOLVED,        #:Optimal
    1 => MOI.OTHER_LIMIT,           #:Not_Finished
    2 => MOI.ITERATION_LIMIT,       #:Iteration_Limit
    3 => MOI.LOCALLY_INFEASIBLE,    #:Infeasible_Problem_Detected
    4 => MOI.NUMERICAL_ERROR)       #:Restoration_needed

    if ret in keys(status_dict)
        status = status_dict[ret]
    else
        status = MOI.OTHER_ERROR
    end

    #set solution
    for (idx,node) in enumerate(model_list)  #set solution values for each model
        local_data = PipsNLP._get_pips_data(node)
        if idx != 1 #all cores have the first stage
            coreid = zeros(Int, 1)
            sc = MPI.Reduce(local_data.coreid, MPI.SUM, root, comm)

            if rank == root
                coreid[1] = sc
            end

            MPI.Bcast!(coreid, length(coreid), root, comm)

            n = zeros(Int, 1)
            n[1] = local_data.n
            MPI.Bcast!(n,length(n),coreid[1], comm)

            if rank != coreid[1]
                local_data.n = n[1]
                local_data.x_sol = zeros(Float64,local_data.n)
            end
            MPI.Bcast!(local_data.x_sol, local_data.n, coreid[1], comm)
        end

        vars = Plasmo.all_variables(node)
        vals = local_data.x_sol
        Plasmo.set_node_primals(node,vars,vals)
        Plasmo.set_node_status(node,status)
        #TODO: set duals
        #TODO: load results into a graph backend optimizer
    end

    if rank == 0
        println(status)
    end

    return status
end  #end pips_nlp_solve

end #end module


text:
module test_spineruntimes_C_InterfaceTestFixture

using Test
using SpineRuntimes: Lib
using .Lib: spAtlas_createFromFile, spAtlas_dispose, spSkeletonJson_create, spSkeletonJson_readSkeletonDataFile, spSkeletonJson_dispose,
            spAnimationStateData_create, spSkeleton_create, spAnimationState_create, spSkeleton_setToSetupPose,
            spAnimationState_setAnimationByName, spAnimationState_addAnimationByName,
            spSkeleton_update, spAnimationState_update, spAnimationState_apply,
            spSkeleton_dispose, spAnimationState_dispose, spAnimationStateData_dispose, spSkeletonData_dispose
using .Lib: spAtlas, spSkeletonJson, spSkeletonData, spAnimationStateData, spAnimation

# https://github.com/EsotericSoftware/spine-runtimes/blob/4.0/spine-c/spine-c-unit-tests/tests/C_InterfaceTestFixture.cpp

function readSkeletonJsonData(jsonName::String, patlas::Ptr{spAtlas})
    pjson::Ptr{spSkeletonJson} = spSkeletonJson_create(patlas)
    pskeletonData::Ptr{spSkeletonData} = spSkeletonJson_readSkeletonDataFile(pjson, Base.unsafe_convert(Ptr{Cchar}, jsonName))
    spSkeletonJson_dispose(pjson)
    return pskeletonData
end

function isnull(x)
    x == C_NULL
end

function enumerateAnimations(pskeletonData)::Vector{spAnimation}
    if !isnull(pskeletonData)
        skeletonData = unsafe_load(pskeletonData)
        animations = unsafe_wrap(Array, skeletonData.animations, skeletonData.animationsCount)
        return unsafe_load.(animations)
    end
    return []
end

macro isnotnull(p)
    name = string(p)
    quote
        ptr = $(esc(p))
        if ptr == C_NULL
            throw(AssertionError(string($name, " != C_NULL")))
        end
    end
end

function testRunner(block, jsonName::String, atlasName::String)
    path::Ptr{Cchar} = Base.unsafe_convert(Ptr{Cchar}, atlasName)
    patlas = spAtlas_createFromFile(path, C_NULL)
    @isnotnull patlas
    atlas::spAtlas = unsafe_load(patlas)
    pages = unsafe_load(atlas.pages)
    pskeletonData = readSkeletonJsonData(jsonName, patlas)
    pstateData::Ptr{spAnimationStateData} = spAnimationStateData_create(pskeletonData)
    stateData = unsafe_load(pstateData)
    stateData.defaultMix = 0.2
    pskeleton = spSkeleton_create(pskeletonData)
    pstate = spAnimationState_create(pstateData)
    spSkeleton_setToSetupPose(pskeleton)
    anims = enumerateAnimations(pskeletonData)
    # spAnimationState_setAnimationByName(pstate, 0, first(anims).name, false)
    # for anim in anims[2:end]
    #     spAnimationState_addAnimationByName(pstate, 0, anim.name, false, 0.0f0)
    # end
    # timeSlice = 1.0f0 / 60.0f0
    # spSkeleton_update(pskeleton, timeSlice)
    # spAnimationState_update(pstate, timeSlice)
    # spAnimationState_apply(pstate, pskeleton)

    block(anims)

    # Dispose Instance
    spSkeleton_dispose(pskeleton)
    spAnimationState_dispose(pstate)

    # Dispose Global
    spAnimationStateData_dispose(pstateData)
    spSkeletonData_dispose(pskeletonData)
    spAtlas_dispose(patlas)
end

const SPINEBOY_JSON  = normpath(@__DIR__, "testdata/spineboy/spineboy-ess.json")
const SPINEBOY_ATLAS = normpath(@__DIR__, "testdata/spineboy/spineboy.atlas")
const RAPTOR_JSON    = normpath(@__DIR__, "testdata/raptor/raptor-pro.json")
const RAPTOR_ATLAS   = normpath(@__DIR__, "testdata/raptor/raptor.atlas")
const GOBLINS_JSON   = normpath(@__DIR__, "testdata/goblins/goblins-pro.json")
const GOBLINS_ATLAS  = normpath(@__DIR__, "testdata/goblins/goblins.atlas")

get_anim_name(x) = unsafe_string(x.name)

testRunner(SPINEBOY_JSON, SPINEBOY_ATLAS) do anims
    @test length(anims) == 7
    @test get_anim_name.(anims) == ["death", "hit", "idle", "jump", "run", "shoot", "walk"]
end
testRunner(RAPTOR_JSON,   RAPTOR_ATLAS) do anims
    @test length(anims) == 5
    @test get_anim_name.(anims) == ["gun-grab", "gun-holster", "jump", "roar", "walk"]
end
testRunner(GOBLINS_JSON,  GOBLINS_ATLAS) do anims
    @test length(anims) == 1
    @test get_anim_name.(anims) == ["walk"]
end

end # module test_spineruntimes_C_InterfaceTestFixture


text:
function http_get_call_with_url(url::String)::PSResult

    # check: is the URL string empty?
    # check: is it a legit URL string?

    # ok, so we are going to make a HTTP GET call with the URL that was passed in -
    response = HTTP.request("GET",url)

    # ok, so let's check if we are getting a 200 back -
    if (response.status == 200)
        return PSResult{String}(String(response.body))
    else
        # create an error, and throw it back to the caller -
        error_message = "http status flag $(response.status) was returned from url $(url)"
        return PSResult{PSError}(PSError(error_message))
    end
end


text:
const examples = [
    (
        "Empty_Example",
        raw"""Project(
            (
                start=0,
                finish=0,
            ),
            [
                :start => :finish,
            ]
        )
        """,
        8,
    ),

    (
        "One_Task_Example",
        raw"""Project(
            (
                start=0,
                a = PertBeta(0,1,2),
                finish=0,
            ),
            [
                :start => :a,
                :a => :finish,
            ]
        )
        """,
        4,
    ),

    (
        "Shortlong_Example",
        raw"""Project(
            (
                start=0,
                short = PertBeta(0,1,2),
                long = PertBeta(10,20,200),
                finish=0,
            ),
            [
                :start .=> [:short, :long];
                [:short, :long] .=> :finish;
            ]
        )
        """,
        4,
    ),

    (
        "Chain_Example",
        raw"""Project(
            (
                start=0,
                a=PertBeta(1,2,3),
                b=PertBeta(10,20,30),
                c=PertBeta(100,200,600),
                finish=0,
            ),
            [
                :start => :a,
                :a => :b,
                :b => :c,
                :c => :finish,
            ]
        )
        """,
        2.5,
    ),

    (
        "Realistic_Example",
        raw"""Project(
            (
                start=0,
                a=PertBeta(2,3,4),
                b=PertBeta(1,2,6),
                c=PertBeta(0,1,2),
                d=PertBeta(0,1,2),
                e=PertBeta(0,1,2),
                f=PertBeta(3.0, 6.5, 13.0),
                g=PertBeta(3.5, 7.5, 14.0),
                h=PertBeta(0, 1, 2),
                j=PertBeta(0, 1, 2),
                finish=0,
            ),
            [
                :start .=> [:a, :b, :c, :d];
                :a => :f;
                :b .=> [:f, :g];
                [:c, :d] .=> :e;
                :e .=> [:f, :g, :h];
                [:f, :g, :h] .=> :j;
                :j => :finish;
            ],
        )
        """,
        2.5,
    ),
]

# Taken from https://github.com/GiovineItalia/Compose.jl/blob/7fde5c28ceb46dae0927f8fdb6b347680eb86387/docs/make.jl#L4-L7
struct SVGJSWritable{T}
    x :: T
end
Base.show(io::IO, m::MIME"text/html", x::SVGJSWritable) = show(io, m, x.x)

function create_demofiles()
    demo_dir = mkpath(joinpath(@__DIR__, "src", "demos",))
    template = String(read(joinpath(@__DIR__, "demo_template.md")))
    for (name, construction, fontsize) in examples
        filled_template = replace(template, "{{NAME}}" => name)
        filled_template = replace(filled_template, "{{CONSTRUCTION}}" => construction)
        filled_template = replace(filled_template, "{{FONTSIZE}}" => fontsize)
        open(joinpath(demo_dir, "$name.md"), "w") do fh
            write(fh, filled_template)
        end
    end
end


text:
include("../graph_utils.jl")

function cell_state_update_graph(cfg::Dict)
    nin = Neurodevo.input_lengths(cfg)[4]
    nout = Neurodevo.output_lengths(cfg)[4]
    mg = get_graph(nin, nout, 7)
    add_node!(mg, nin+1, CGP.Config.f_const, p=0.0)
    add_node!(mg, nin+2, CGP.Config.f_max, x=nin+1, y=nin)
    add_node!(mg, nin+3, CGP.Config.f_const, p=0.9)
    add_node!(mg, nin+4, CGP.Config.f_const, p=0.1)
    add_node!(mg, nin+5, CGP.Config.f_mult, x=nin+3, y=cfg["n_cell_params"]+1)
    add_node!(mg, nin+6, CGP.Config.f_mult, x=nin+4, y=nin+2)
    add_node!(mg, nin+7, CGP.Config.f_sum, x=nin+5, y=nin+6)
    set_outputs!(mg, nin, nout, 7, vcat(nin+7, cfg["n_cell_params"] .+ 2:cfg["n_cell_state"]))
    mg
end


text:
module LatticeWalks

import ..LatticeVars.init!, ..LatticeVars.get_num_sites_visited
using ..Lattices
using ..WalksBase
import ..WalksBase: get_position, get_time, get_nsteps, step!
using ..Actors

export AbstractLatticeWalk, LatticeWalk, get_num_sites_visited

abstract type AbstractLatticeWalk <: AbstractWalkGeneral
end

function init!(lattice_walk::AbstractLatticeWalk)
    init!(lattice_walk.lattice)
    WalksBase.init!(lattice_walk.walk)
    return nothing
end

###
### LatticeWalk
###

"""
    struct LatticeWalk

Represents a walk on a lattice.
"""
struct LatticeWalk{L <: AbstractLattice, W <: AbstractWalk} <: AbstractLatticeWalk
    lattice::L
    walk::W
end

for f in (:get_position, :get_time, :get_nsteps)
    @eval ($f)(lw::LatticeWalk) = ($f)(lw.walk)
end

for f in (:get_num_sites_visited,)
    @eval ($f)(lw::LatticeWalk) = ($f)(lw.lattice)
end

Base.getindex(lattice_walk::LatticeWalk, inds...) = getindex(lattice_walk.lattice, inds...)

LatticeWalk(lat, walk) = LatticeWalk{typeof(lat), typeof(walk)}(lat, walk)

function Base.show(io::IO, lw::LatticeWalk)
    println(io, "LatticeWalk")
    show(io, lw.lattice)
    println(io)
    show(io, lw.walk)
    println(io)
end

function step!(lattice_walk::LatticeWalk{<:Any, <:Any})
    result = step!(lattice_walk.walk, lattice_walk.lattice)
    result && return result
    return result
end

function step!(walk::WalkB, lattice::TrapsLattice)
    pos = get_position(walk)
    jump_time = lattice[1, pos]
    addto_time!(walk, jump_time)
    step!(walk)
    return true
end

function step!(lattice_walk::LatticeWalk{<:Any, <:MortalWalk})
    result = step!(lattice_walk.walk, lattice_walk.lattice)
    result && return result
    set_status!(lattice_walk.walk, false)
    return result
end

function step!(walk::MortalWalk, lattice::TrapsLattice)
    pos = get_position(walk)
    jump_time = lattice[1, pos]
    decay_time = lattice[2, pos]
    if decay_time < jump_time
        addto_time!(walk, decay_time)
        return false
    else
        addto_time!(walk, jump_time)
        step!(walk)
        return true
    end
end

#function step!(walk::WalkF{X,Y,Z}, lattice::SitePercolationLattice) where {X,Y,Z}
function step!(walk::WalkB, lattice::SitePercolationLattice)
    disp = try_step_increment(walk)
    @show disp
    newpos = get_position(walk) + disp
    @show newpos
    site_state = lattice[1, newpos]
    @show site_state
    if site_state
        set_position!(walk, newpos)
    end
    return true
end

end  # module LatticeWalks


text:
using Pkg; Pkg.activate("~/.julia/dev/Tullio/benchmarks/02/")
pkg"add Tullio, LoopVectorization, LinearAlgebra, BenchmarkTools, Einsum, TensorOperations, MKL_jll, https://github.com/haampie/FastTranspose.jl, https://github.com/mcabbott/ArrayMeta.jl, StructArrays"
pkg"up"

Threads.nthreads() == 6

using Tullio, LoopVectorization, Einsum, TensorOperations, FastTranspose, ArrayMeta,Libdl, Test, StructArrays, BenchmarkTools, Plots

TVERSION = VersionNumber(Pkg.TOML.parsefile(joinpath(pkgdir(Tullio), "Project.toml"))["version"])

# adapted from https://github.com/haampie/FastTranspose.jl/blob/master/benchmark/mkl.jl
using MKL_jll
const libMKL = Libdl.dlopen(MKL_jll.libmkl_rt)
const MKL_domatcopy = Libdl.dlsym(libMKL, :mkl_domatcopy)
const MKL_somatcopy = Libdl.dlsym(libMKL, :mkl_somatcopy)
function mkl_matcopy!(B::Matrix{Float64}, A::Matrix{Float64}, alpha = 1.0)
    m, n = size(A)
    ordering = 'C'
    trans = 'T'
    ccall(MKL_domatcopy, Cvoid, (Ref{Cchar}, Ref{Cchar}, Ref{Csize_t}, Ref{Csize_t}, Ref{Cdouble}, Ptr{Float64}, Ref{Csize_t}, Ptr{Float64}, Ref{Csize_t}),
           ordering, trans, m, n, alpha, A, m, B, n)
    return nothing
end
function mkl_matcopy!(B::Matrix{Float32}, A::Matrix{Float32}, alpha = 1.0)
    m, n = size(A)
    ordering = 'C'
    trans = 'T'
    ccall(MKL_somatcopy, Cvoid, (Ref{Cchar}, Ref{Cchar}, Ref{Csize_t}, Ref{Csize_t}, Ref{Cdouble}, Ptr{Float32}, Ref{Csize_t}, Ptr{Float32}, Ref{Csize_t}),
           ordering, trans, m, n, alpha, A, m, B, n)
    return nothing
end

for T in [Float64, Float32]
    local A,B = rand(T, 3,3), rand(T, 3,3)
    mkl_matcopy!(B,A)
    # @test B == A' # fails for Float32?
end

# matrix transpose

base_transpose!(y,x) = permutedims!(y, x, (2,1))
base_lazy_transpose!(y,x) = copyto!(y, transpose(x))
einsum_transpose!(y,x) = @einsum y[i,j] = x[j,i]
arraymeta_transpose!(y,x) = @arrayop y[i,j] = x[j,i]
tensor_transpose!(y,x) = @tensor y[i,j] = x[j,i]
avx_transpose!(y,x) = @tullio y[i,j] = x[j,i] threads=false tensor=false
tullio_transpose!(y,x) = @tullio y[i,j] = x[j,i] tensor=false

functions2 = (base_transpose!, base_lazy_transpose!,
    einsum_transpose!, arraymeta_transpose!, tensor_transpose!,
    avx_transpose!, tullio_transpose!,
    recursive_transpose!, mkl_matcopy!)
sizes2 = sort(vcat(vec((2 .^ (4:12))' .+ [-1,0,1]), [25, 50, 100, 150, 200, 1500]))

# 3-array permutedims

base_312!(y, x) = permutedims!(y, x, (3,1,2))
base_lazy_312!(y, x) = copyto!(y, PermutedDimsArray(x, (3,1,2)))
einsum_312!(y, x) = @einsum y[c,a,b] = x[a,b,c]
arraymeta_312!(y, x) = @arrayop  y[c,a,b] = x[a,b,c]
tensor_312!(y, x) = @tensor y[c,a,b] = x[a,b,c]
avx_312!(y, x) = @tullio y[c,a,b] = x[a,b,c] threads=false tensor=false
tullio_312!(y, x) = @tullio y[c,a,b] = x[a,b,c] tensor=false

functions3 = (base_312!, base_lazy_312!,
    einsum_312!, arraymeta_312!, tensor_312!,
    avx_312!, tullio_312!)
sizes3 = sort(vcat(vec((2 .^ (3:9))' .+ [-1,0,1] ), [25, 50, 100, 150, 200]))

# 4-array permutedims

base_4321!(y, x) = permutedims!(y, x, (4,3,2,1))
base_lazy_4321!(y, x) = copyto!(y, PermutedDimsArray(x, (4,3,2,1)))
einsum_4321!(y, x) = @einsum y[d,c,b,a] = x[a,b,c,d]
arraymeta_4321!(y, x) = @arrayop y[d,c,b,a] = x[a,b,c,d]
tensor_4321!(y, x) = @tensor y[d,c,b,a] = x[a,b,c,d]
avx_4321!(y, x) = @tullio y[d,c,b,a] = x[a,b,c,d] threads=false tensor=false
tullio_4321!(y, x) = @tullio y[d,c,b,a] = x[a,b,c,d] tensor=false

functions4 = (base_4321!, base_lazy_4321!,
    einsum_4321!, arraymeta_4321!, tensor_4321!,
    avx_4321!, tullio_4321!)
sizes4 = sort(vcat(vec((2 .^ (2:6))' .+ [-1,0,1]), [25, 50, 99, 100, 101])) # crashes on 129!

# based on https://github.com/chriselrod/PaddedMatrices.jl/blob/master/benchmark/blasbench.jl

function benchmark_fun!(f!, C, A, force_belapsed = false, reference = nothing)
    tmin = @elapsed f!(C, A)
    if force_belapsed || 2tmin < BenchmarkTools.DEFAULT_PARAMETERS.seconds
        tmin = min(tmin, @belapsed $f!($C, $A))
    elseif tmin < BenchmarkTools.DEFAULT_PARAMETERS.seconds
        tmin = min(tmin, @elapsed f!(C, A))
        if tmin < 2BenchmarkTools.DEFAULT_PARAMETERS.seconds
            tmin = min(tmin, @elapsed f!(C, A))
        end
    end
    isnothing(reference) || @assert C ≈ reference
    tmin
end

function runbench(funs::Tuple, dims::Int=2, sizes=[10,30,100])
    (StructVector ∘ map)(sizes) do n
        A = rand(ntuple(_->n, dims)...)
        C0, C = similar(A), similar(A)
        funs[1](C0, A)
        times = map(funs) do f!
            t = benchmark_fun!(f!, C, A)
            C ≈ C0 || @error "disagreement for $f! at size $n"
            t
        end
        nt = NamedTuple{map(Symbol, funs)}(times)
        res = (size=n, length=length(C), nt...)
        @show res
    end
end

r0 = runbench((base_transpose!, mkl_matcopy!, tullio_transpose!), 2, [10,30,100])

function makeplot(res, title="")
    plt = plot()
    for lab in propertynames(res)[3:end]
        times = getproperty(res, lab)
        flops = 1e-9 * res.length ./ times # was wrong in 0.2.0 plots

        if startswith(string(lab), "tullio")
            plot!(plt, res.size, flops, lab=string(lab), m=:circle)
        else
            plot!(plt, res.size, flops, lab=string(lab))
        end
    end
    plot!(plt, yaxis=("numbers / ns", ([0.5,1,2,4,8], ["1/2","1",2,4,8]), :log10), xaxis=("size per dimension", :log10), legend=:bottomleft)
    plot!(plt, title = title * "Julia " * string(VERSION) * ", Intel " * split(Sys.cpu_info()[1].model)[3])
end

makeplot(r0)
# savefig(joinpath("~/.julia/dev/Tullio", "benchmarks/02/trash.png"))

res2 = runbench(functions2, 2, sizes2)
plot2 = makeplot(res2, "Tullio $TVERSION, ")
savefig(plot2, joinpath("~/.julia/dev/Tullio", "benchmarks/02/transpose-$TVERSION-$VERSION.png"))

res3 = runbench(functions3, 3, sizes3)
plot3 = makeplot(res3, "Tullio $TVERSION, ")
savefig(plot3, joinpath("~/.julia/dev/Tullio", "benchmarks/02/permute3-$TVERSION-$VERSION.png"))

res4 = runbench(functions4, 4, sizes4)
plot4 = makeplot(res4, "Tullio $TVERSION, ")
savefig(plot4, joinpath("~/.julia/dev/Tullio", "benchmarks/02/permute4-$TVERSION-$VERSION.png"))


text:
struct CuDeviceArray{T,N} <: AbstractArray{T,N}
    shape::NTuple{N,Int}
    ptr::Ptr{T}
end

@inline function Base.getindex(A::CuDeviceArray{T}, index::Integer) where {T}
    @boundscheck checkbounds(A, index)
    Base.unsafe_load(A.ptr, index)::T
end

Base.size(g::CuDeviceArray) = g.shape

Base.checkbounds(::CuDeviceArray, I...) = nothing

function kernel(a, i)
    return a[i]
end

code_llvm(stdout, kernel, Tuple{CuDeviceArray{Float32,1}, Int})

using TinyCassette
struct GPUctx end

code_llvm(stdout, TinyCassette.execute,
          Tuple{GPUctx, typeof(kernel), CuDeviceArray{Float32,1}, Int})


text:
using Distributions
using Test

@testset "Log of Beta-binomial distribution" begin
    d = BetaBinomial(50, 0.2,0.6)

    for k in Base.OneTo(50)
        p  = pdf(d, k)
        lp = logpdf(d, k)
        @test_approx_eq lp log(p)
        @test insupport(d, k)
    end
    @test !insupport(d, 51)
end


text:
function leitor(P)

end

function geraP1(P)
    n = length(P)
    saida = zeros(n)
    for i = 1:n
        saida[j] = P[n - j + 1]
    end
    return saida
end

function geraP2(P)
    n = length(P)
    saida = zeros(n)
    for i = 1:n
        saida[j] = P[n - j + 1]
    end
    return saida

end

function Lagrange(n, k, B, an)
    return 1.0 + (B/an)^(1.0/(n - k))
end

function encontrarK(P)
    n = length(P)
    j = 1
    k = 0
    acabou = false
    while ~acabou
        if P[j] < 0
            k = n - j
            acabou = true
        else
            k = 0
        end

        j = j + 1

        if j > n
            acabou  = true
        end

    end
    return k
end

function acharB(P)

    b = minimum(P)

    if b > 0
        b = 0
    else
        b = -b
    end
end

function Limite(P)
    n = length(P) - 1
    a = P[1]
    k = encontrarK(P)
    B = acharB(P)
    return Lagrange(n, k, B, a)
end

function final(P)
    L = Limite(P)
    L1 = 1.0/Limite(gerarP1(P))
    L2 = Limite(gerarP2(P))
    L3 = -1.0/Limite(geraP2(P1))
    return L2, L3, L1, L
end


text:
using FinancialMonteCarlo
import FinancialMonteCarlo.pricer, FinancialMonteCarlo.AbstractMethod
abstract type AbstractIntegrationMethod <: FinancialMonteCarlo.AbstractMethod end
abstract type AbstractFFTMethod <: AbstractIntegrationMethod end
abstract type AbstractIntegralMethod <: AbstractIntegrationMethod end


text:
include("FEM1D.jl")

using LinearAlgebra, SparseArrays

function BPX_half_matrix(L::Integer, l::Integer)
    ψ = Array{Float64, 1}(undef, 2^(L-l))
    χ = Array{Float64, 1}(undef, 2^(L-l))
    ψ .= range(1, 2^(L-l), step=1)/2^(L-l)
    χ .= 1 .- ψ
    n = 2^l
    Id = spdiagm(0 => ones(n))
    S = spdiagm(-1 => ones(n-1))
    B = (kron(Id, ψ) + kron(S, χ))[1:(end-1), 1:(end-1)]
    return B
end

function Schwarz_matrix_3D(A::SparseMatrixCSC{Float64, Int64}, L::Integer)
    B = spdiagm(0 => diag(A).^-1)
    for l=1:(L-1)
        P_x = BPX_half_matrix(L, l)
        P = kron(P_x, P_x, P_x)
        if l == 1
            A_P = sparse(inv(Array(transpose(P)*A*P)))
        else
            A_P = spdiagm(0 => diag(transpose(P)*A*P).^-1)
        end
        B += P*A_P*transpose(P)
    end
    return B
end

function Laplace_3D(L::Integer)
    N = 2^L - 1
    A = spdiagm(0 => -2*ones(N), 1 => ones(N-1), -1 => ones(N-1))
    D = spdiagm(0 => ones(N))
    T = [[D, D, A], [A, D, D], [D, A, D]]
    M = sum([kron(t[1], kron(t[2], t[3])) for t = T])
    return M
end

function CG_solver!(A::SparseMatrixCSC{Float64, Int64}, b::Array{<:Real, 1}, x::Array{<:Real, 1}, tol::Real)
    r::Array{Float64, 1} = b - A*x
    p::Array{Float64, 1} = r .+ 0.0
    p_A::Array{Float64, 1} = r .+ 0.0
    β::Float64 = 1.0
    α::Float64 = 1.0
    γ_0::Float64 = norm(r, Inf)
    error::Float64 = 1.0
    γ::Float64 = γ_0
    N_sweeps::Integer = 0
    while error/γ_0 > tol^2
        γ = dot(r, r)
        p_A .= A*p
        α = γ/dot(p, p_A)
        x .+= α*p
        r .-= α*(p_A)
        β = dot(r, r)/γ
        p .= r + β*p
        error = norm(b - A*x, Inf)
        N_sweeps += 1
    end
    return N_sweeps
end

function PCG_solver!(A::SparseMatrixCSC{Float64, Int64}, B::SparseMatrixCSC{Float64, Int64}, b::Array{<:Real, 1}, x::Array{<:Real, 1}, tol::Real)
    r::Array{Float64, 1} = b - A*x
    p::Array{Float64, 1} = r .+ 0.0
    p_A::Array{Float64, 1} = r .+ 0.0
    z::Array{Float64, 1} = r .+ 0.0
    β::Float64 = 1.0
    α::Float64 = 1.0
    γ_0::Float64 = norm(r, Inf)
    error::Float64 = 1.0
    γ::Float64 = γ_0
    N_sweeps::Integer = 0
    while error/γ_0 > tol^2
        γ = dot(r, z)
        p_A .= A*p
        α = γ/dot(p, p_A)
        x .+= α*p
        r .-= α*(p_A)
        z .= B*r
        β = dot(r, z)/γ
        p .= z + β*p
        error = norm(b - A*x, Inf)
        N_sweeps += 1
    end
    return N_sweeps
end

function test_CG(L::Integer)
    A = Laplace_3D(L)
    N = 2^L - 1
    h = 1/2^L
    x = randn(N^3)
    b = randn(N^3)
    N_sweeps = CG_solver!(A, b, x, 1e-6)
    println(N_sweeps)
end

function test_PCG(L::Integer)
    A = Laplace_3D(L)
    B = Schwarz_matrix_3D(A, L)
    N = 2^L - 1
    h = 1/2^L
    x = randn(N^3)
    b = randn(N^3)
    N_sweeps = PCG_solver!(A, B, b, x, 1e-6)
    println(N_sweeps)
end

function test_sparse(L::Integer)
    A = Laplace_3D(L)
    N = 2^L - 1
    b = randn(N^3)
    F = lu(A)
    x = F \ b
    return nothing
end

max_stars_repo_name: VLSF/neuralBPX

max_issues_repo_name: VLSF/neuralBPX

max_forks_repo_name: VLSF/neuralBPX