text:
using JuliaBlackBoard
using Test

@testset "JuliaBlackBoard.jl" begin
    d = tempdir()
    cdir = @__DIR__
    cp(joinpath(cdir, "test-setup.jl"), joinpath(d,"test-setup.jl"), force=true)
    include(joinpath(d, "test-setup.jl"))
end


text:
using Miletus.DayCounts
using Dates

@test 0.497724380567 ≈ yearfraction(ISDAActualActual(), Date(2003,November,1), Date(2004, May, 1))
@test 0.50           ≈ yearfraction(ISMAActualActual(), Date(2003,November,1), Date(2004, May, 1))
@test 0.497267759563 ≈ yearfraction(AFBActualActual(), Date(2003,November,1), Date(2004, May, 1))

@test 0.410958904110 ≈ yearfraction(ISDAActualActual(), Date(1999,February,1), Date(1999, July, 1))
@test 0.410958904110 ≈ yearfraction(ISMAActualActual(), Date(1999,February,1), Date(1999, July, 1), Date(1998, July, 1), Date(1999, July, 1))
@test 0.410958904110 ≈ yearfraction(AFBActualActual(), Date(1999,February,1), Date(1999, July, 1))

@test 1.001377348600 ≈ yearfraction(ISDAActualActual(), Date(1999,July,1), Date(2000, July, 1))
@test 1.0            ≈ yearfraction(ISMAActualActual(), Date(1999,July,1), Date(2000, July, 1), Date(1999, July, 1), Date(2000, July, 1))
@test 1.0            ≈ yearfraction(AFBActualActual(), Date(1999,July,1), Date(2000, July, 1))

@test 0.915068493151 ≈ yearfraction(ISDAActualActual(), Date(2002,August,15), Date(2003, July, 15))
@test 0.915760869565 ≈ yearfraction(ISMAActualActual(), Date(2002,August,15), Date(2003, July, 15), Date(2003, January, 15), Date(2003, July, 15))
@test 0.915068493151 ≈ yearfraction(AFBActualActual(), Date(2002,August,15), Date(2003, July, 15))

@test 0.504004790778 ≈ yearfraction(ISDAActualActual(), Date(2003,July,15), Date(2004, January, 15))
@test 0.5            ≈ yearfraction(ISMAActualActual(), Date(2003,July,15), Date(2004, January, 15), Date(2003, July, 15), Date(2004, January, 15))
@test 0.504109589041 ≈ yearfraction(AFBActualActual(), Date(2003,July,15), Date(2004, January, 15))

@test 0.503892506924 ≈ yearfraction(ISDAActualActual(), Date(1999,July,30), Date(2000, January, 30))
@test 0.5            ≈ yearfraction(ISMAActualActual(), Date(1999,July,30), Date(2000, January, 30), Date(1999,July,30), Date(2000, January, 30))
@test 0.504109589041 ≈ yearfraction(AFBActualActual(), Date(1999,July,30), Date(2000, January, 30))

@test 0.415300546448 ≈ yearfraction(ISDAActualActual(), Date(2000,January,30), Date(2000, June, 30))
@test 0.417582417582 ≈ yearfraction(ISMAActualActual(), Date(2000,January,30), Date(2000, June, 30), Date(2000,January,30), Date(2000, July, 30))
@test 0.415300546448 ≈ yearfraction(AFBActualActual(), Date(2000,January,30), Date(2000, June, 30))

using BusinessDays
c=BusinessDays.UnitedKingdom()
adjust(ModPreceding(), c, Date(2016,5,2)) == Date(2016, 5, 3)

ext: jl


text:
using SafeTestsets
@safetestset "My f tests" begin include("my_f_test.jl") end
@safetestset "fVAI() tests" begin include("fVAI_test.jl") end


text:
# This file was generated, do not modify it. # hide
@load OneHotEncoder
MultinomialClassifier = @load MultinomialClassifier pkg="MLJLinearModels"

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
center = 5
sigma = 2

mf = GaussianMF(center, sigma)

# Evaluation tests
@assert Fuzzy.eval(mf, center + sigma) == exp(-0.5)
@assert Fuzzy.eval(mf, center) == 1
@assert Fuzzy.eval(mf, center + sigma) == Fuzzy.eval(mf, center - sigma)

# Mean finding tests
@assert Fuzzy.mean_at(mf, 1) == Fuzzy.mean_at(mf, 0.6) == Fuzzy.mean_at(mf, 0.3) == Fuzzy.mean_at(mf, 0) == center


text:
module Util

using Cxx

function cstrings(jlstrings)
    x = icxx"vector<string>();"
    for _ in jlstrings
        icxx"$x.push_back($(pointer(bytestring(_))));"
    end
    x
end

function cstring(jlstring)
    b = bytestring(jlstring)
    icxx"string($(pointer(b)), $(length(b)));"
end

cstring(jlstring::Symbol) = cstring(string(jlstring))

function cstring_map(jlmap)
    m = icxx"map<string, string>();"
    for (key, val) in jlmap
        icxx"$m[$(cstring(key))] = $(cstring(val));"
    end
    m
end

function cstring_vector(strings)
    x = icxx"vector<string>();"
    for _ in strings
        icxx"$x.push_back($(cstring(_)));"
    end
    x
end

function jlstring(cstr)
    icxx"$cstr.data();" |> bytestring
end

function jl_vector(cvec, T)
    b = icxx"$cvec->begin();"
    e = icxx"$cvec->end();"
    x = T[]
    while icxx"$b != $e"
        push!(x, T(icxx"auto q=*$b; return q;"))
        b = icxx"$b+1;"
    end
    x
end

end


text:
# This file was generated, do not modify it. # hide
rmsl(ŷ, y[test])


text:
using Statistics
using Random
import ChainRulesCore: rrule, NoTangent

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

# TODO: migrate to a common broadcasting rrule

# from Zygote:
# https://github.com/FluxML/Zygote.jl/blob/d5be4d5ca80e79278d714eaac15ca71904a262e3/src/lib/array.jl#L177-L185
struct StaticGetter{i} end
(::StaticGetter{i})(v) where {i} = v[i]

@generated function _unzip(tuples, ::Val{N}) where {N}
  Expr(:tuple, (:(map($(StaticGetter{i}()), tuples)) for i ∈ 1:N)...)
end

function unzip(tuples)
  N = length(first(tuples))
  _unzip(tuples, Val(N))
end

function rrule(::typeof(Broadcast.broadcasted), f::F, args...) where F
    ys, pbs = unzip(rrule.(f, args...))
    function pullback(Δ)
        dxs = map((pb, Δ) -> pb(Δ), pbs, Δ)
        return NoTangent(), unzip(dxs)...
    end
    return ys, pullback
end

Yota.update_chainrules_primitives!()

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

obj(Y, X, b) = mean((Y .- X * b) .^ 2.0) # objective to minimize

# modified example from GradDescent.jl documentation
@testset "linreg" begin

    Random.seed!(1) # set seed
    n = 1000 # number of observations
    d = 10   # number of covariates
    X = randn(n, d) # simulated covariates
    b = randn(d)    # generated coefficients
    ϵ = randn(n) * 0.1 # noise
    Y = X * b + ϵ # observed outcome

    epochs = 100 # number of epochs

    θ = randn(d) # initialize model parameters

    loss_first, g = grad(obj, Y, X, θ)
    loss_last = loss_first
    for i in 1:epochs
        loss_last, g = grad(obj, Y, X, θ)
        # println("Epoch: $i; loss = $loss_last")
        δ = 0.01 * g[4]
        θ -= δ
    end
    @test loss_last < loss_first / 10   # loss reduced at least 10 times

end

mutable struct Linear2{T}
    W::AbstractMatrix{T}
    b::AbstractVector{T}
end

linear_loss(m::Linear2, X, y) = mean((m.W * X .+ m.b .- y) .^ 2.0)

@testset "linreg 2" begin

    n_vars = 10
    n_out = 2
    n_obs = 1000
    X = randn(n_vars, n_obs)
    W_true = rand(n_out, n_vars)
    b_true = rand(n_out)
    ϵ = randn(n_out, n_obs) * 0.1 # noise
    Y = W_true * X .+ b_true .+ ϵ

    m = Linear2(rand(n_out, n_vars), rand(n_out))

    epochs = 100

    for i in 1:epochs
        val, g = grad(linear_loss, m, X, Y)
        # println("Epoch: $i; loss = $val")
        update!(m, g[2], (x, gx) -> x - 0.1 * gx)
    end

    @test isapprox(m.W, W_true; atol=0.1)
    @test isapprox(m.b, b_true; atol=0.1)

end


text:
load("kernels.jl")

dest1 = Array(Int, (2,3))
dest2 = Array(Int, (2,3))
B = [1 2 3
     4 5 6]
C = [1 0 1
     0 1 0]
D = [ 0  0  0
     10 10 10]

# code = quote
#     A = B.*C + D
#     dest1[...] = A
#     dest2[...] = A + C
# end

@kernel let
    A = B.*C + D
    dest1[...] = A
    dest2[...] = A + C
end

function f(B,C,D)
    A = B.*C + D
    dest1 = A
    dest2 = A + C
    (dest1, dest2)
end
dest1f, dest2f = f(B,C,D)

println("kernel:     dest1 =\n$dest1")
println("Array ops:  dest1 =\n$dest1f")
println()
println("kernel:     dest2 =\n$dest2")
println("Array ops:  dest2 =\n$dest2f")
println()

function timeit(n::Int)
    A = Array(Float, n, n)
    temp = Array(Float, n, n)
    B = rand(n, n)
    C = rand(n, n)
    D = rand(n, n)

    print("Array routines:     "); tic()
    B.*C + D
    ta = toc()

    print("For loops:          "); tic()
    for j=1:n,i=1:n
        temp[i,j] = B[i,j].*C[i,j]
    end
    for j=1:n,i=1:n
        A[i,j] = temp[i,j] + D[i,j]
    end
    tfor = toc()

    print("Handwritten kernel: "); tic()
    for j=1:n,i=1:n
        A[i,j] = B[i,j].*C[i,j] + D[i,j]
    end
    tmk = toc()

    print("Kernel:             "); tic()
    @kernel let
        A[...] = B.*C + D
    end
    tk = toc()
    (ta, tmk, tfor, tk)
end

timeit(1000)

max_stars_repo_name: toivoh/julia-kernels

max_issues_repo_name: toivoh/julia-kernels

max_forks_repo_name: toivoh/julia-kernels


text:
using Documenter, BernsteinExpansions

DocMeta.setdocmeta!(BernsteinExpansions, :DocTestSetup,
                   :(using BernsteinExpansions); recursive=true)

makedocs(
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true",
        assets = ["assets/aligned.css"]),
    sitename = "BernsteinExpansions.jl",
    doctest = true,
    strict = true,
    pages = [
        "Home" => "index.md",
        "Manual" => Any[
            "Introduction" => "man/introduction.md",
            "Range Enclosure" => "man/range.md",
        "Benchmarks" => "man/benchmarks.md"],
        "Library" => Any[
            "Types" => "lib/types.md",
            "Methods" => "lib/methods.md"],
        "About" => "about.md"
    ]
)

deploydocs(
    repo = "github.com/JuliaReach/BernsteinExpansions.jl.git",
    push_preview=true
)

max_stars_repo_name: JuliaReach/BernsteinExpansions

max_issues_repo_name: JuliaReach/BernsteinExpansions

max_forks_repo_name: JuliaReach/BernsteinExpansions


text:
### bootstrap.jl
#
# Copyright (C) 2016, 2017 Mosè Giordano.
#
# Maintainer: Mosè Giordano <mose AT gnu DOT org>
# Keywords: periodogram, lomb scargle, bootstrapping
#
# This file is a part of LombScargle.jl.
#
# License is BSD 3-clause "New" or "Revised".
#
### Commentary
#
# This file contains facilities to perform bootstrapping and to calculate
# false-alarm probability and its inverse.
#
### Code:

using Random

struct Bootstrap{T<:AbstractFloat}
    p::Vector{T} # Vector of highest peaks
end

const _default_rng = if VERSION < v"1.3"
    () -> Random.GLOBAL_RNG
else
    Random.default_rng
end

bootstrap(rng::AbstractRNG, N::Integer, p::PeriodogramPlan) =
    Bootstrap(sort!([maximum(normalize!(_periodogram!(shuffle(rng, p.times), p), p)) for _ in 1:N], rev = true))

bootstrap(rng::AbstractRNG, N::Integer, t::AbstractVector{<:Real}, rest...; kwargs...) =
    bootstrap(rng, N, plan(t, rest...; kwargs...))

bootstrap(N::Integer, t::AbstractVector{<:Real}, rest...; kwargs...) =
    bootstrap(_default_rng(), N, plan(t, rest...; kwargs...))

"""
    LombScargle.bootstrap(N::Integer,
                          times::AbstractVector{Real},
                          signal::AbstractVector{Real},
                          errors::AbstractVector{Real}=ones(signal); ...)

Create `N` bootstrap samples, perform the Lomb–Scargle analysis on them, and
store all the highest peaks for each one in a `LombScargle.Bootstrap` object.
All the arguments after `N` are passed around to [`lombscargle`](@ref).
"""
bootstrap(::Integer, ::AbstractVector{<:Real})

"""
    LombScargle.bootstrap(N::Integer, plan::PeriodogramPlan)

Create `N` bootstrap samples, perform the Lomb–Scargle analysis on them for the given
`plan`, and store all the highest peaks for each one in a `LombScargle.Bootstrap` object.

See documentation of [`LombScargle.plan`](@ref) for how to plan a Lomb–Scargle periodogram.
"""
bootstrap(::Integer, ::PeriodogramPlan)

"""
    fap(b::Bootstrap, power::Real)

Return the false-alarm probability for `power` in the bootstrap sample `b`.

Its inverse is the [`fapinv`](@ref) function.
"""
fap(b::Bootstrap{<:AbstractFloat}, power::Real) =
    length(findall(x -> x >= power, b.p))/length(b.p)

"""
    fapinv(b::Bootstrap, prob::Real)

Return the power value whose false-alarm probability is `prob` in the bootstrap
sample `b`.

It returns `NaN` if the requested probability is too low and the power cannot be
determined with the bootstrap sample `b`.  In this case, you should enlarge your
bootstrap sample so that `N*fap` can be rounded to an integer larger than or
equal to 1.

This is the inverse of [`fap`](@ref) function.
"""
fapinv(b::Bootstrap{<:AbstractFloat}, prob::Real) =
    get(b.p, round(Int, length(b.p)*prob), NaN)


text:
# This implements an efficient sampler for the univariate truncated
# standard normal distribution from:

# Y. Li and S. Ghosh (2015)
# "Efficient Sampling Methods for Truncated Multivariate Normal and Student-t
#  Distributions Subject to Linear Inequality Constraints"
# Journal of Statistical Theory & Practice
# doi: 10.1080/15598608.2014.996690
# http://dx.doi.org/10.1080/15598608.2014.996690

immutable TruncatedNormalSampler <: Sampleable{Univariate, Continuous}
  a::Float64       # Lower bound
  b::Float64       # Upper bound
  m::ASCIIString   # Sampling method to use:
                   #   NR: Normal rejection sampling
                   #   HR: Half-Normal rejection sampling
                   #   UR: Uniform rejection sampling
                   #   ER: Translated exponential rejection sampling
  sgn::Int64       # Change sign of result to be returned?
end

function TruncatedNormalSampler(a::Float64, b::Float64)
  m::ASCIIString
  sgn = 1
  a0 = 0.2570

  # Needs to work for untruncated regions
  if a == -Inf && b == Inf
    return TruncatedNormalSampler(a, b, "NR", sgn)
  end

  # Symmetric case, return negative of result
  if b <= 0  || a == -Inf
    sgn = -1
    olda = a
    a = -b
    b = -olda
  end

  b1 = a + sqrt(π/2) * exp(a*a/2)
  b2 = a + 2 / (a + sqrt(a*a + 4)) * exp((a*a - a*sqrt(a*a +4))/4 + .5)

  # Determine method to use
  a < 0 ?
    b == Inf ? m = "NR" :
      b - a <= sqrt(2π) ? m = "UR" : m = "NR" :
    b == Inf ?
      a < a0 ? m = "HR" : m = "ER" :
      a < a0 ?
        b <= b1 ? m = "UR" : m = "HR" :
        b <= b2 ? m = "UR" : m = "ER"

  return TruncatedNormalSampler(a, b, m, sgn)
end

function rand(s::TruncatedNormalSampler)
#   a = s.a
#   b = s.b
#   sgn = s.sgn
  if s.m == "NR"
    while true
      x = randn()
      x >= s.a && x <= s.b && return x*s.sgn
    end
  elseif s.m == "HR"
    while true
      x = abs(randn())
      if x >= s.a && x <= s.b
        return x*s.sgn
      end
    end
  elseif s.m == "UR"
    Zd = Normal(0, 1)
    Uab = Uniform(s.a, s.b)
    pc = (cdf(Zd, s.b) - cdf(Zd, s.a))
    if s.a >= 0
      pc *= exp(s.a*s.a/2)
    elseif s.b <= 0
      pc *= exp(s.b*s.b/2)
    end
    u = 0.0
    x = 0.0
    while true
      u = rand()
      x = rand(Uab)
      if u < exp(-x*x/2) * pc
        return x*s.sgn
      end
    end
  elseif s.m == "ER"
    λ = (s.a + sqrt(s.a*s.a + 4)) / 2
    u = 0.0
    x = 0.0
    p = 0.0
    while true
      x = randexp() / λ + s.a
      u = rand()
      p = exp(-(x-λ)^2 / 2)
      u < p && x <= s.b && return x*s.sgn
    end
  end
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 5b182574-5656-4035-85c1-89d1c92ff719
begin
	import Pkg
	Pkg.develop("PoreMatMod")
end

# ╔═╡ 1a73bdc4-9394-4212-9ae8-3b8654a496c0
using PoreMatMod

# ╔═╡ f05cdf05-1661-41c7-863d-15a436791ac4
using PoreMatMod.ExampleHelpers

# ╔═╡ 0b644231-c25a-4d9f-895d-16fc612613ec
md"# Azepine"

# ╔═╡ 916b9bc9-c6c9-444b-9281-e3709c248b75
begin
	benzodiazepine_parent = moiety("benzodiazepine.xyz")
	remove_bonds!(benzodiazepine_parent)
	local new_box = replicate(unit_cube(), (20, 20, 20))
	local new_atoms = Frac(Cart(benzodiazepine_parent.atoms, benzodiazepine_parent.box), new_box)
	local new_charges = Frac(Cart(benzodiazepine_parent.charges, benzodiazepine_parent.box), new_box)
	benzodiazepine_parent = Crystal(benzodiazepine_parent.name, new_box, new_atoms, new_charges)
	infer_bonds!(benzodiazepine_parent, false)
	benzo_fragment = moiety("benzo_fragment.xyz")
	chloro_benzo_fragment = moiety("chloro_benzo_fragment.xyz")
	diazepine_fragment = moiety("diazepine_fragment.xyz")
	N_methyl_diazepine = moiety("N_methyl_diazepine.xyz")
end

# ╔═╡ da1f22fb-8e16-439b-934a-280ed31635df
begin
	intermediate = replace(benzodiazepine_parent, diazepine_fragment => N_methyl_diazepine)
	diazepam = replace(intermediate, benzo_fragment => chloro_benzo_fragment)
end

# ╔═╡ 7d73fa2e-9ee5-4002-97a6-496f4d186512
begin
	local temp = deepcopy(intermediate)
	translate_by!(temp.atoms.coords, Frac([0.5, 0.5, 0.5]))
	wrap!(temp)
	write_cif(temp, "intermediate.cif")

	local temp = deepcopy(diazepam)
	translate_by!(temp.atoms.coords, Frac([0.5, 0.5, 0.5]))
	wrap!(temp)
	write_cif(temp, "diazepam.cif")
end

# ╔═╡ 6fe41bf5-cf36-4313-8886-cfca2825c6cd
md"# Slab"

# ╔═╡ faff7bea-0d3a-45fe-8a88-b288b390e35b
function read_poscar(filename::String)::Crystal
	filedata = readlines(filename)
	box_matrix = Matrix(reduce(hcat, [parse.(Float64, row) for row in split.(filedata[3:5])])')
	species = Symbol.(split(filedata[6]))
	atom_counts = parse.(Int, split(filedata[7]))
	nb_atoms = sum(atom_counts)
	species_vec = reduce(vcat, [[element for _ in 1:atom_counts[i]] for (i, element) in enumerate(species)])

	box = Box(box_matrix)
	coords = Frac(reduce(hcat, [parse.(Float64, row) for row in split.(filedata[9:nb_atoms+8])]))
	atoms = Atoms(nb_atoms, species_vec, coords)
	charges = Charges(nb_atoms, zeros(nb_atoms), coords)

	return Crystal(filename, box, atoms, charges)
end

# ╔═╡ db771098-8097-4748-85c7-ece4faed4e43
begin
	poscar = read_poscar("data/crystals/POSCAR")
	infer_bonds!(poscar, true)
end

# ╔═╡ 35c1c1d5-ea01-4d15-80c3-63330744b037
write_cif(poscar, "poscar.cif")

# ╔═╡ b09ea1c5-67b5-4953-8fef-c5144f09186b
hydrated_Pd2 = moiety("hydrated_Pd2.xyz")

# ╔═╡ b7a6a162-39cc-4137-b803-94c2a5326b6e
OA_Pd2 = moiety("OA_Pd2.xyz")

# ╔═╡ e232cfc8-a954-4089-be59-a9fa5b3a2736
oxidative_addition = replace(poscar, hydrated_Pd2 => OA_Pd2)

# ╔═╡ 908ab709-c12a-455b-994b-89e600d47b06
write_cif(oxidative_addition, "oxidative_addition.cif")

# ╔═╡ Cell order:
# ╟─5b182574-5656-4035-85c1-89d1c92ff719
# ╠═1a73bdc4-9394-4212-9ae8-3b8654a496c0
# ╠═f05cdf05-1661-41c7-863d-15a436791ac4
# ╟─0b644231-c25a-4d9f-895d-16fc612613ec
# ╠═916b9bc9-c6c9-444b-9281-e3709c248b75
# ╠═da1f22fb-8e16-439b-934a-280ed31635df
# ╠═7d73fa2e-9ee5-4002-97a6-496f4d186512
# ╟─6fe41bf5-cf36-4313-8886-cfca2825c6cd
# ╠═faff7bea-0d3a-45fe-8a88-b288b390e35b
# ╠═db771098-8097-4748-85c7-ece4faed4e43
# ╠═35c1c1d5-ea01-4d15-80c3-63330744b037
# ╠═b09ea1c5-67b5-4953-8fef-c5144f09186b
# ╠═b7a6a162-39cc-4137-b803-94c2a5326b6e
# ╠═e232cfc8-a954-4089-be59-a9fa5b3a2736
# ╠═908ab709-c12a-455b-994b-89e600d47b06


text:
function η_karnik(u::Matrix{FuzzyNumber}, d2::Matrix{FuzzyNumber}; m::Real=2.0)
	m > 1 || error("fuzzifier m ∈ (1, ∞), got $m")
	levels = u[1, 1].levels
	num_levels = length(levels)
	N, c = size(u)

	η = Vector{FuzzyNumber}(undef, c)
	for i = 1:c
		grades = Vector{Interval}(undef, num_levels)
		for (lvl, α) in enumerate(levels)
			u_cut = cut(FuzzyVector(u[:, i]), α)
			d2_cut = cut(FuzzyVector(d2[:, i]), α)
			c_left = km_iwa(d2_cut, u_cut; bound="lower", m=m)
			c_right = km_iwa(d2_cut, u_cut; bound="upper", m=m)
			grades[lvl] = Interval(c_left, c_right)
		end
		η[i] = FuzzyNumber(levels, grades)
	end
	η
end

function u_lpcm(X⃗::FuzzyVector, C⃗::FuzzyVector, η::Real; m::Real=2.0)
	m > 1 || error("fuzzifier m ∈ (1, ∞), got $m")
	η > 0 || error("η must be suitable positive number")
	levels = X⃗[1].levels
	num_levels = length(levels)

	grades = Vector{Interval}(undef, num_levels)
	D2 = FuzzySets.d_interval(X⃗, C⃗, squared=true)
	for (lvl, α) in enumerate(levels)
		D2_cut = cut(D2, α)
		grades[lvl] = u_lpcm(D2_cut, η; m=m)
	end
	u = FuzzyNumber(levels, grades)
	u
end

function u_lpcm(d2::Interval, η::Real; m::Real=2.0)
	m > 1 || error("fuzzifier m ∈ (1, ∞)")
	η > 0 || error("η must be suitable positive number")
	power = 1 / (m - 1)

	u = Interval(1) / (Interval(1) + (d2 / η)^power)
	u
end


text:
using Base.Test
using UnitQuaternions

####################################################################################################
# Constructors
####################################################################################################

print("Testing constructors...")

q = UnitQuaternion([1, 1, 1], 1)
for i = 1:4
    @test_approx_eq q[i] 0.5
end

@test_throws ErrorException UnitQuaternion([1, 1, 1, 1], 1)

@test UnitQuaternion() == UnitQuaternion([0, 0, 0], 1)

for i = 1:100
    a, b, c, d = 5randn(4)
    p = UnitQuaternion(a, b, c, d)
    q = UnitQuaternion([a, b, c, d])
    @test p == UnitQuaternion([a, b, c], d)
    @test q == UnitQuaternion([a, b, c], d)
    @test p == q
    @test_approx_eq sqrt(q.η^2 + norm(q.ϵ)^2) 1
end

@test UnitQuaternion([0, 0, 0]) == UnitQuaternion()
@test UnitQuaternion([π, 0, 0]) == UnitQuaternion(1, 0, 0, 0)
@test UnitQuaternion([2π, 0, 0]) == UnitQuaternion(0, 0, 0, -1)
@test UnitQuaternion([0, π/2, 0]) == UnitQuaternion(0, sqrt(2)/2, 0, sqrt(2)/2)
@test UnitQuaternion([0, 0, π/3]) == UnitQuaternion(0, 0, 0.5, sqrt(3)/2)

println(" done.")

####################################################################################################
# Operators
####################################################################################################

print("Testing operators...")

@test_approx_eq +(UnitQuaternion()) eye(4)
@test_approx_eq +(UnitQuaternion(1, 0, 0, 0)) [0 0 0 1; 0 0 1 0; 0 -1 0 0; -1 0 0 0]
@test_approx_eq +(UnitQuaternion(1, 1, 1, 1)) [ 0.5  0.5 -0.5  0.5;
                                               -0.5  0.5  0.5  0.5;
                                                0.5 -0.5  0.5  0.5;
                                               -0.5 -0.5 -0.5  0.5;]

@test UnitQuaternion(1, 2, 3, 4) + UnitQuaternion() == UnitQuaternion(1, 2, 3, 4)
@test UnitQuaternion() + UnitQuaternion(1, 2, 3, 4) == UnitQuaternion(1, 2, 3, 4)
@test UnitQuaternion(1, 2, 3, 4) ⊕ UnitQuaternion() == UnitQuaternion(1, 2, 3, 4)
@test UnitQuaternion() ⊕ UnitQuaternion(1, 2, 3, 4) == UnitQuaternion(1, 2, 3, 4)

q = UnitQuaternion()
@test q^2 == q
@test q^3 == q
q = UnitQuaternion(sind(90/2), 0, 0, cosd(90/2))
@test q^2 == UnitQuaternion(sind(180/2), 0, 0, cosd(180/2))
@test q^0.5 == UnitQuaternion(sind(45/2), 0, 0, cosd(45/2))

p = UnitQuaternion(sqrt(2)/2, 0, 0, sqrt(2)/2)
q = UnitQuaternion(0, 0, sqrt(2)/2, sqrt(2)/2)
@test p + q == UnitQuaternion(0.5, 0.5, 0.5, 0.5)
@test q ⊕ p == UnitQuaternion(0.5, 0.5, 0.5, 0.5)
@test q + p == UnitQuaternion(0.5, -0.5, 0.5, 0.5)
@test p ⊕ q == UnitQuaternion(0.5, -0.5, 0.5, 0.5)

@test q ⊞ [0, 0, 0] == q
@test p ⊞ (q ⊟ p) == q
@test_approx_eq (q ⊞ [0.1, 0.2, 0.3]) ⊟ q [0.1, 0.2, 0.3]

r = inv(p) + q
@test r == UnitQuaternion(-0.5, -0.5, 0.5, 0.5)
@test_approx_eq q ⊟ p 2atan2(norm(r.ϵ), r.η) * r.ϵ / norm(r.ϵ)

@test p ⊞ [0, 0, π/2] == UnitQuaternion(0.5, 0.5, 0.5, 0.5)

@test_throws BoundsError q[5]
@test_approx_eq q[1] 0
@test_approx_eq q[4] sqrt(2)/2

println(" done.")

####################################################################################################
# Methods
####################################################################################################

print("Testing methods...")

q = UnitQuaternion(0, 0, sqrt(2)/2, sqrt(2)/2)

@test_approx_eq angle(q) π/2

@test_approx_eq axis(q) [0, 0, 1]

qmean = UnitQuaternion(sind(90/2), 0, 0, cosd(90/2))
v = [UnitQuaternion(sind(60/2), 0, 0, cosd(60/2)), UnitQuaternion(sind(120/2), 0, 0, cosd(120/2))]
cov = covariance(v, qmean, [1, 1])
@test_approx_eq rad2deg(sqrt(cov[1])) 30

@test inv(q) == UnitQuaternion(0, 0, -sqrt(2)/2, sqrt(2)/2)

qs = [UnitQuaternion(), UnitQuaternion()]
weights = [rand(), rand()]
@test mean(qs, weights) == UnitQuaternion()
@test mean(qs) == UnitQuaternion()
qs = [UnitQuaternion(1, 0, 0, 0), UnitQuaternion(1, 0, 0, 0)]
weights = [rand(), rand()]
@test mean(qs, weights) == UnitQuaternion(1, 0, 0, 0)
@test mean(qs) == UnitQuaternion(1, 0, 0, 0)
qs = [UnitQuaternion(1.1, 0, 0, 0), UnitQuaternion(0.9, 0, 0, 0)]
weights = [1, 1]
@test mean(qs, weights) == UnitQuaternion(1, 0, 0, 0)
@test mean(qs) == UnitQuaternion(1, 0, 0, 0)
qs = [UnitQuaternion(0, 1.2, 0, 0), UnitQuaternion(0, 0.8, 0, 0)]
weights = [1, 1]
@test mean(qs, weights) == UnitQuaternion(0, 1, 0, 0)
@test mean(qs) == UnitQuaternion(0, 1, 0, 0)
qs = [UnitQuaternion(0, 0, 1.3, 0), UnitQuaternion(0, 0, 0.7, 0)]
weights = [1, 1]
@test mean(qs, weights) == UnitQuaternion(0, 0, 1, 0)
@test mean(qs) == UnitQuaternion(0, 0, 1, 0)
qs = [UnitQuaternion(0, 0, 0, 1.4), UnitQuaternion(0, 0, 0, 0.6)]
weights = [1, 1]
@test mean(qs, weights) == UnitQuaternion(0, 0, 0, 1)
@test mean(qs) == UnitQuaternion(0, 0, 0, 1)

@test_approx_eq log(UnitQuaternion()) [0, 0, 0]
@test_approx_eq log(UnitQuaternion(1, 1, 1, 0)) π/2 * [sqrt(3)/3, sqrt(3)/3, sqrt(3)/3]
@test_approx_eq log(q) [0, 0, π/4]

@test_approx_eq rotatevector(q, [1, 0, 0]) [0, 1, 0]
p = UnitQuaternion(0.5, 0.5, 0.5, 0.5)
@test_approx_eq rotatevector(p, [1, 0, 0]) [0, 1, 0]

for i = 1:10
    r = rand(3)
    p = UnitQuaternion(rand(4))
    q = UnitQuaternion(rand(4))
    @test_approx_eq rotateframe(q, r) rotationmatrix(q) * r
    @test_approx_eq rotateframe(p, r) rotationmatrix(p) * r
    @test_approx_eq rotationmatrix(p + q) rotationmatrix(p) * rotationmatrix(q)
end

p = UnitQuaternion(sind(30/2), 0, 0, cosd(30/2))
q = UnitQuaternion(sind(60/2), 0, 0, cosd(60/2))
@test slerp(p, q, 0.5) == UnitQuaternion(sind(45/2), 0, 0, cosd(45/2))
p = UnitQuaternion(0, sind(10/2), 0, cosd(10/2))
q = UnitQuaternion(0, sind(90/2), 0, cosd(90/2))
@test slerp(p, q, 0.1) == UnitQuaternion(0, sind(18/2), 0, cosd(18/2))
@test slerp(p, q, 0.4) == UnitQuaternion(0, sind(42/2), 0, cosd(42/2))
@test slerp(p, q, 1.1) == UnitQuaternion(0, sind(98/2), 0, cosd(98/2))

q = UnitQuaternion(0, 0, sqrt(2)/2, sqrt(2)/2)
@test_approx_eq vector(UnitQuaternion()) [0, 0, 0, 1]
@test_approx_eq vector(q) [0, 0, sqrt(2)/2, sqrt(2)/2]

println(" done.")


text:
function serve_update_page_from_dir(folder)

    folder = realpath(folder)
    @assert isdir(folder) "$folder is not a valid directory."
    refimages_name = last(splitpath(realpath(folder)))
    @info "Refimage set name is $refimages_name"

    router = HTTP.Router()

    function receive_update(req)
        data = JSON3.read(req.body)
        images = data["images"]
        tag = data["tag"]

        tempdir = tempname()
        recorded_folder = joinpath(folder, "recorded")
        reference_folder = joinpath(folder, "reference")

        @info "Copying reference folder to \"$tempdir\""
        cp(reference_folder, tempdir)

        for image in images
            @info "Overwriting \"$image\" in new reference folder"
            copy_filepath = joinpath(tempdir, image)
            copy_dir = splitdir(copy_filepath)[1]
            # make the path in case a new refimage is in a not yet existing folder
            mkpath(copy_dir)
            cp(joinpath(recorded_folder, image), copy_filepath, force = true)
        end

        @info "Uploading updated reference images under tag \"$tag\""
        try
            upload_reference_images(tempdir, tag; name = refimages_name)

            HTTP.Response(200, "Upload successful")
        catch e
            showerror(stdout, e, catch_backtrace())
            HTTP.Response(404)
        end
    end

    function serve_local_file(req)
        if req.target == "/"
            s = read(normpath(joinpath(dirname(pathof(ReferenceTests)), "reference_images.html")), String)
            s = replace(s, "DEFAULT_TAG" => "'$(last_major_version())'")
            return HTTP.Response(200, s)
        end
        file = HTTP.unescapeuri(req.target[2:end])
        filepath = normpath(joinpath(folder, file))
        # check that we don't go outside of the artifact folder
        if !startswith(filepath, folder)
            @info "$file leads to $filepath which is outside of the artifact folder."
            return HTTP.Response(404)
        end

        if !isfile(filepath)
            return HTTP.Response(404)
        else
            return HTTP.Response(200, read(filepath))
        end
    end

    HTTP.@register(router, "POST", "/", receive_update)
    HTTP.@register(router, "GET", "/", serve_local_file)

    @info "Starting server. Open localhost:8000 in your browser to view."
    HTTP.serve(router, HTTP.Sockets.localhost, 8000)
end

function serve_update_page(; commit = nothing, pr = nothing)
    authget(url) = HTTP.get(url, Dict("Authorization" => "token $(ENV["GITHUB_TOKEN"])"))

    commit !== nothing && pr !== nothing && error("Keyword arguments `commit` and `pr` can't be set at once.")
    if pr !== nothing
        prinfo = JSON3.read(authget("https://api.github.com/repos/JuliaPlots/Makie.jl/pulls/$pr").body)
        headsha = prinfo["head"]["sha"]
        @info "PR is $pr, using last commit hash $headsha"
    elseif commit !== nothing
        headsha = commit
    else
        error("You have to specify either the keyword argument `commit` or `pr`.")
    end

    checksinfo = JSON3.read(authget("https://api.github.com/repos/JuliaPlots/Makie.jl/commits/$headsha/check-runs").body)

    checkruns = filter(checksinfo["check_runs"]) do x
        right_combination = any(["GLMakie", "WGLMakie", "CairoMakie"]) do package
            occursin(package, x["name"]) && occursin("1.6", x["name"])
        end
        if right_combination
            job = JSON3.read(authget("https://api.github.com/repos/JuliaPlots/Makie.jl/actions/jobs/$(x["id"])").body)
            run = JSON3.read(authget(job["run_url"]).body)
            if run["status"] != "completed"
                @info "$(x["name"])'s run hasn't completed yet, no artifacts will be available."
                false
            else
                true
            end
        else
            false
        end
    end

    if isempty(checkruns)
        error("No check runs fit the criteria, check if something about names or versions might have changed.")
    end

    menu = REPL.TerminalMenus.RadioMenu([x["name"] for x in checkruns])
    choice = REPL.TerminalMenus.request("Choose a workflow run:", menu)

    if choice == -1
        error("Cancelled")
    end
    check = checkruns[choice]

    job = JSON3.read(authget("https://api.github.com/repos/JuliaPlots/Makie.jl/actions/jobs/$(check["id"])").body)
    run = JSON3.read(authget(job["run_url"]).body)

    artifacts = JSON3.read(authget(run["artifacts_url"]).body)["artifacts"]
    for a in artifacts
        if endswith(a["name"], "1.6")
            @info "Choosing artifact $(a["name"])"
            download_url = a["archive_download_url"]
            @info "Downloading artifact from $download_url"
            filepath = Downloads.download(download_url, headers = Dict("Authorization" => "token $(ENV["GITHUB_TOKEN"])"))
            @info "Download successful"
            tmpdir = mktempdir()
            unzip(filepath, tmpdir)

            folders = readdir(tmpdir)
            if length(folders) == 0
                error("No folder in zip")
            elseif length(folders) == 1
                folder = only(folders)
            else
                menu = REPL.TerminalMenus.RadioMenu(folders, pagesize=4)
                choice = REPL.TerminalMenus.request("Choose a reference image set:", menu)

                if choice == -1
                    error("Cancelled")
                end
                folder = folders[choice]
            end

            @info "Serving update page for folder $folder."
            serve_update_page_from_dir(joinpath(tmpdir, folder))
            return
        end
    end
    error("""
        Reference image artifacts found for commit $headsha and job id $(check["id"]).
        This could be because the job's workflow run ($(job["run_url"])) has not completed, yet.
        Artifacts are only available for complete runs.
    """)
end

function unzip(file, exdir = "")
    fileFullPath = isabspath(file) ?  file : joinpath(pwd(),file)
    basePath = dirname(fileFullPath)
    outPath = (exdir == "" ? basePath : (isabspath(exdir) ? exdir : joinpath(pwd(),exdir)))
    isdir(outPath) ? "" : mkdir(outPath)
    @info "Extracting zip file $file to $outPath"
    zarchive = ZipFile.Reader(fileFullPath)
    for f in zarchive.files
        fullFilePath = joinpath(outPath,f.name)
        if (endswith(f.name,"/") || endswith(f.name,"\\"))
            mkdir(fullFilePath)
        else
            write(fullFilePath, read(f))
        end
    end
    close(zarchive)
    @info "Extracted zip file"
end


text:
using BenchmarkConfigSweeps

BenchmarkConfigSweeps.run(
    joinpath(@__DIR__, "build"),
    joinpath(@__DIR__, "benchmarks.jl"),
    BenchmarkConfigSweeps.nthreads.([1]),
)


text:
using PauliStrings
using Test

@testset "PauliString - Structure & Basic Arithmetic" begin

	p_string = PauliString(4,[1,3],[2,3],3.3)

	#
	#	check construction
	#
	p_string.N == 4
	p_string.sites == [1,3]
	p_string.baseIdx == [2,3]
	@test coef(p_string) == 3.3

	#
	#	check set_coef!
	#
	set_coef!(p_string, 2.0)
	@test coef(p_string) == 2.0

	#
	#	check multiplication
	#
	p_start = deepcopy(p_string)
	p_string = 3.0*p_string
	@test coef(p_string) == 6.0

	p_string = p_string*2.0
	@test coef(p_string) == 12.0

	#
	#	check equal and parallel:
	#	parallel means same site same baseIdx different coef
	#
	p_string∥p_start == true
	@test iszero(p_string == p_start)

	# check addition
	p_start += p_start # 4
	p_start *= 3.0
	@test isone(p_string == p_start)

	#
	#	check conj
	#
	p_string *= 1.0im
	p_string_conj = conj(p_string)

	@test coef(p_string_conj) == -12.0im

	conj!(p_string)
	@test isone(p_string_conj == p_string)

end

@testset "PauliStrings - Convertion to SparseArray" begin

	@show nothing

end

@testset "PauliStrings - Commutation & Multiplication" begin

	# test multiplication
	p1 = PauliString(1,[1],[2],1.0)
	p2 = PauliString(1,[1],[3],1.0)
	p3 = PauliString(1,[1],[1],1.0im)
	@test p3 ==  p1*p2
	@test p3 == conj(p2*p1)

	p1 = PauliString(2,[1,2],[2,1],1.0)
	p2 = PauliString(2,[1],[3],1.0)
	p3 = PauliString(2,[1,2],[1,1],1.0im)

	@test p3 ==  p1*p2
	@test p3 == conj(p2*p1)

	#commutation
	@test commutePauli(p1,p2) == 2*p3

	p1 = PauliString(2,[1,2],[2,1],1.0)
	p2 = PauliString(2,[1,2],[3,2],1.0)
	p3 = PauliString(2,[1,2],[1,3],-1.0)

	@test p3 ==  p1*p2
	@test p3 == p2*p1

	# commutation
	@test commutePauli(p1,p2) == nothing

end

max_stars_repo_name: Gawatz/PauliStrings

max_issues_repo_name: Gawatz/PauliStrings

max_forks_repo_name: Gawatz/PauliStrings


text:
module simulation
using Distributions, DataFrames, Compat
using MathProgBase, Gurobi
using ..BSBP

""" Generates data from a normal distribution.
Args:
    n (int):        Number of observations.
    beta (float):   The upper limit of the range to generate, from 0 to `n` - 1.
    sigma (float):  The upper limit of the range to generate, from 0 to `n` - 1.
    typ (int):      1 for heteroskedasticity
Yields:
    y:              outcome
    datax:          covariates
"""
function simdata(n::Integer,theta::AbstractArray,sigma::AbstractArray,typ::Integer,seed::Integer)

    #Set seed
    if seed != nothing
        rng = MersenneTwister(seed)
        reg = Distributions.rand(rng,MvNormal(sigma),n)'
    else
        reg = Distributions.rand(MvNormal(sigma),n)'
    end

    k = length(theta)
    datax = [reg[:,1] ones(n,1) reg[:,2:k-1]]
    if typ == 1
        sigmaW = reg[:,1] .+ reg[:,2]
        e = 0.25*(ones(n) .+ 2*sigmaW.^2 .+ sigmaW.^4).*randn(n)
    else
        e = 0.25*randn(n)
    end
    y = Int.((datax*theta).>=e)
    return y,datax
end

""" Simulates best subset binary prediction on simulated data.
Args:
    warm_start (int,optional):  Set to 1 for warm start strategy. Default = 0.
    tau (float,optional):       Tuning paramater to construct the refined bound used in the warm start approach. Default = 1.5.
    mio (int,optional):         Set to 1 for using Method 1 for the MIO formulation. Set to 2 for Method 2. Default to 1.
    q (int, optional):          Value of the variable selection bound. Default = 1.
    N (int, optional):          Size of training sample. Default = 100.
    N_val (int, optional):      Size of validation sample. Default = 5000.
    R (int, optional):          Simulation repetitions. Default = 10.
    typ (int, optional):        Set to 1 for heteroskedastic error design. Set to 0 for homoskedastic error design. Default = 1.
"""
function Simulation(;warm_start=0,tau=1.5,mio=1,q=1,p=10,N=10,N_val=50,R=1,typ=1,beta0=1,rho=0.25,maxT=0,seed=1)

    # theta
    if typ==1
        beta_s = -1.5
    else
        beta_s = -0.35
    end

    beta = [beta0; [0]; beta_s; zeros(p-1,1)]

    K=length(beta)

    bhat=zeros(K-1,R)
    # sigma matrix
    row = zeros(p+1)
    for i in 1:p+1
        row[i]= rho^(i-1)
    end
    sigma = zeros(p+1, p+1)
        for i in 1:p+1
            sigma[i,:] = [row[i:-1:1];row[2:p-i+2]]
        end

    y, datax = simdata(N,beta,sigma,typ,seed)

    gap=zeros(R) # MIO gap
    rtime=zeros(R) # MIO running time
    ncount=zeros(R) # MIO node count
    score=zeros(R) # MIO score

    DGP_score=zeros(R) # in-sample score at the DGP parameter vector
    val_score=zeros(R) # in-sample score at the estimated parameter vector
    DGP_score_test=zeros(R) # out-of-sample score at the DGP parameter vector
    val_score_test=zeros(R) # out-of-sample score at the estimated parameter vector

    bnd=[-10*ones(K-1,1) 10*ones(K-1,1)]
    bnd_h = zeros(size(bhat)[1],2,R)

    if q>=1 && p>N
        tol=minimum([0.5*sqrt((1+q)*log(p)*N),0.05*N])  # early stopping rule
    else
        tol=0
    end

    println("Start")
    for i = 1:R
        println(i)

        y,datax = simdata(N,beta,sigma,typ,seed)

        try
            if warm_start == 1 # warm start MIO
                println("Start")
                bhat[:,i],score[i],gap[i],rtime[i],ncount[i]  = WarmStartMaxScore(y,datax[:,1:2],datax[:,3:size(datax)[2]],beta0,q,maxT,tol,bnd,mio,tau)
            else # cold start MIO
                bhat[:,i],score[i],gap[i],rtime[i],ncount[i]  = MaxScore(y,datax[:,1:2],datax[:,3:size(datax)[2]],beta0,q,maxT,tol,bnd,mio)
                println("Done")
            end
            catch y
                println(y.msg)
                throw(ErrorException("Error"))
        end

        DGP_score[i] = mean(Int.(y .== Int.(datax*beta.>=0)))
        val_score[i] = mean(Int.(y .== Int.(datax*[beta0;bhat[:,i]].>=0)))
        if N_val>0
            y_val,datax_val = simdata(N_val,beta,sigma,typ,seed)
            DGP_score_test[i] = mean(Int.(y_val .== Int.(datax_val*beta.>=0)))
            val_score_test[i] = mean(Int.(y_val .== Int.(datax_val*[beta0;bhat[:,i]].>=0)))
         end

    end

    println(mean([val_score[:,:] DGP_score[:,:] val_score_test[:,:] DGP_score_test[:,:]],2))
    println(mean([(val_score./DGP_score)[:,:] (val_score_test./DGP_score_test)[:,:]],2))
    return score,gap,rtime,ncount
end

export Simulation
end


text:
include("../functions/Functions.jl")
include("../evaluate/Evaluate.jl")
include("../opt/GradientDescent.jl")
include("../data/data.jl")

X, y = getOldData()

lambda = 0.00001
iter = 1000

x0 = zeros(size(X,2),1)

f = f_logreg(X,y,x0,lambda)
g = g_logreg(X,y,x0,lambda)

#f = f_svm(X,y,x0,lambda)
#g = sub_g_svm()

@time (theta, status,vals1,stops1 ) = qn(x0,f,g,1e-12,iter,20,20,1e-4,.9,"bt","gd")

grad = g(theta)
grad_norm = norm(grad, Inf)
@printf("Inf-Norm of gradient is %f.\n",grad_norm)
println(evaluate(X,y,theta,"logreg",false))

max_stars_repo_name: zagazao/num-opt

max_issues_repo_name: zagazao/num-opt

max_forks_repo_name: zagazao/num-opt


text:
"""

    DeltaKIMolBlockData{FT<:AbstractFloat}

Represent the internuclear couplings between a pair of molecules, assumed to be computed via
a Morse potential, with the possibility of Bloch periodicity.

See also: [`DeltaKISystemData{FT}`](@ref)

"""
struct DeltaKIMolBlockData{FT<:AbstractFloat}

    mol1::Integer
    mol2::Integer
    k0molblock::Array{FT, 2}
    a0molblock::Array{FT, 2}
    Einfmolblock::Array{FT, 2}
    DeltaKImolblock::SharedArray{FT, 2}

    numblocks1::Integer
    blocklist1::UnitRange{<:Integer}
    numblocks2::Integer
    blocklist2::UnitRange{<:Integer}

end

"""

    DeltaKISystemData{FT<:AbstractFloat}

Represent the internuclear couplings between every pair of molecules, assumed to be computed
via a Morse potential, with the possibility of Bloch periodicity, and represent the diagonal
elements of that internuclear spring constant matrix for the full system based on the
off-diagonal blocks.

See also: [`DeltaKIMolBlockData{FT}`](@ref)

"""
struct DeltaKISystemData{FT<:AbstractFloat}

    DeltaKImolblockdataarray::Array{DeltaKIMolBlockData{FT}, 1}
    DeltaKIdiagelements::Array{FT, 1}
    numcouplings::Integer

end

"""

    keffMorse(r::FT, k0::FT, a0::FT, Einf::FT) where FT<:AbstractFloat

Compute and return the effective Morse spring constant between two particles at distance `r`
(assumed to be nonnegative), with equilibrium spring constant `k0` and separation `a0` and
with dissociation energy `Einf`. If ``U(r)`` is the potential function, then this returns
the quantity ``(1/(r - a_{0}))\\partial U/\\partial r``.

"""
function keffMorse(r::FT, k0::FT, a0::FT, Einf::FT) where FT<:AbstractFloat

    if (k0 == zero(k0) || a0 == zero(a0) || Einf == zero(Einf))
        return zero(r);
    elseif (r == a0)
        return -1*k0;
    end

    x = sqrt(k0/(2*Einf))*(r - a0);
    return k0*exp(-1*x)*expm1(-1*x)/x;

end

"""

    readDeltaKISystemData(mollistfilename::AbstractString,
                          myAllMolData::MolSystem{OneMol_WithPhonons{FT}}) where FT<:AbstractFloat

Read a configuration file `mollistfilename` (a string) to extract data detailing which pairs
of molecules have nontrivial internuclear couplings between them, handling duplicate entries
appropriately.

"""
function readDeltaKISystemData(mollistfilename::AbstractString,
                               myAllMolData::MolSystem{OneMol_WithPhonons{FT}}) where FT<:AbstractFloat

    molstrarray = strip.(chomp.(readlines(mollistfilename)));
    couplingstartinds = findall(map(u -> startswith(u, "COUPLING"), molstrarray));
    couplingendinds = findall(molstrarray .== "ENDCOUPLING");
    numcouplings = length(couplingendinds);
    if (length(couplingstartinds) != numcouplings)
        error(string(mollistfilename, ": number of COUPLING and ENDCOUPLING statements",
                     " should be the same (may be zero)"));
    end

    # if there are no couplings, return a structure with an empty array
    if (numcouplings == 0)
        return DeltaKISystemData{FT}(Array{DeltaKIMolBlockData{FT}, 1}(undef, 0),
                                     zeros(FT, myAllMolData.cumnumatomslist[end]),
                                     numcouplings);
    end

    DeltaKImolblockdataarray = Array{DeltaKIMolBlockData{FT}, 1}(undef, numcouplings);
    molpairarray = Array{Integer, 2}(undef, numcouplings, 2);
    for nn=1:numcouplings

        currcouplingstrarray = strip.(chomp.(molstrarray[couplingstartinds[nn]+1:couplingendinds[nn]-1]));
        DeltaKImolblockdataarray[nn] = readDeltaKIBlockData(currcouplingstrarray,
                                                            myAllMolData, mollistfilename);
        molpairarray[nn, 1] = DeltaKImolblockdataarray[nn].mol1;
        molpairarray[nn, 2] = DeltaKImolblockdataarray[nn].mol2;

    end

    # need to figure out which blocks have been duplicated or have self-couplings,
    # in order to remove them
    rowisdupe = falses(numcouplings);

    for nn1=2:numcouplings

        for nn2=1:(nn1-1)

            if (molpairarray[nn1, 1] == molpairarray[nn2, 1] &&
                molpairarray[nn1, 2] == molpairarray[nn2, 2])

                rowisdupe[nn1] = true;
                break;

            end

        end

    end
    deleteidxarr = findall(rowisdupe .| (molpairarray[:, 1] .== molpairarray[:, 2]));
    if (length(deleteidxarr) > 0)

        deleteat!(DeltaKImolblockdataarray, deleteidxarr);
        println(string("If ", mollistfilename, " has multiple consecutive pairs of ",
                       "COUPLING and ENDCOUPLING statements, and there are multiple such ",
                       "pairs where the lines of the form: \"MOL1 \${mol1label}\" and ",
                       "\"MOL2 \${mol2label}\", have the same (or transposed) molecular ",
                       "label pairs mol1label and mol2label, only the first occurrence ",
                       "will be used, even if later occurrences specify different input ",
                       "data files; additionally, any self-coupling blocks specifying ",
                       "\${mol1label} to be the same as \${mol2label} will be ignored"));

    end

    return DeltaKISystemData{FT}(DeltaKImolblockdataarray,
                                 zeros(FT, myAllMolData.cumnumatomslist[end]),
                                 length(DeltaKImolblockdataarray));

end

function readDeltaKIBlockData(couplingstrarray::Array{<:AbstractString, 1},
                              myAllMolData::MolSystem{OneMol_WithPhonons{FT}},
                              mollistfilename::AbstractString) where FT<:AbstractFloat

    mol1, mol2, flipmol = extractMolPair(couplingstrarray, myAllMolData, mollistfilename);
    k0, a0, Einf = extractGeneralCouplingData(couplingstrarray, mollistfilename, FT)
    if (flipmol)
        k0 = transpose(k0) .+ zero(eltype(k0));
        a0 = transpose(a0) .+ zero(eltype(a0));
        Einf = transpose(Einf) .+ zero(eltype(Einf));
    end

    numblocks1 = div(size(k0, 1), myAllMolData.numatomslist[mol1]);
    numblocks2 = div(size(k0, 2), myAllMolData.numatomslist[mol2]);
    blocklist1 = div(1 - numblocks1, 2):div(numblocks1 - 1, 2);
    blocklist2 = div(1 - numblocks2, 2):div(numblocks2 - 1, 2);

    return DeltaKIMolBlockData{FT}(mol1, mol2, k0, a0, Einf, zero(k0),
                                   numblocks1, blocklist1, numblocks2, blocklist2);

end

function extractMolPair(couplingstrarray::Array{<:AbstractString, 1},
                        myAllMolData::MolSystem{OneMol_WithPhonons{FT}},
                        mollistfilename::AbstractString) where FT<:AbstractFloat

    curridx1 = findfirst(map(u -> startswith(u, "MOL1 "), couplingstrarray));
    if (curridx1 == nothing)
        error(string("If ", mollistfilename, " has at least one consecutive pair of ",
                     "COUPLING and ENDCOUPLING statements, then it must contain a line of ",
                     "the form: MOL1 \${mol1label}"));
    end
    curridx2 = findfirst(map(u -> startswith(u, "MOL2 "), couplingstrarray));
    if (curridx2 == nothing)
        error(string("If ", mollistfilename, " has at least one consecutive pair of ",
                     "COUPLING and ENDCOUPLING statements, then it must contain a line of ",
                     "the form: MOL2 \${mol2label}"));
    end

    mol1label = replace(couplingstrarray[curridx1], "MOL1 " => "");
    mol2label = replace(couplingstrarray[curridx2], "MOL2 " => "");

    mol1 = findfirst(myAllMolData.molLabels .== mol1label);
    if (mol1 == nothing)
        error(string("If ", mollistfilename, " has at least one consecutive pair of ",
                     "COUPLING and ENDCOUPLING statements, then in the line of ",
                     "the form: \"MOL1 \${mol1label}\", mol1label must correspond to a ",
                     "valid molecular label from elsewhere in the same file"));
    end
    mol2 = findfirst(myAllMolData.molLabels .== mol2label);
    if (mol2 == nothing)
        error(string("If ", mollistfilename, " has at least one consecutive pair of ",
                     "COUPLING and ENDCOUPLING statements, then in the line of ",
                     "the form: \"MOL2 \${mol2label}\", mol1label must correspond to a ",
                     "valid molecular label from elsewhere in the same file"));
    end

    flipmol = (mol2 < mol1);
    if (mol2 == mol1)
       error(string("If ", mollistfilename, " has at least one consecutive pair of ",
                     "COUPLING and ENDCOUPLING statements, then in the lines of ",
                     "the form: \"MOL1 \${mol1label}\" and \"MOL2 \${mol2label}\", ",
                     "mol1label and mol2label must be distinct molecular labels"));
    end
    if (flipmol)
        println(string("If ", mollistfilename, " has at least one consecutive pair of ",
                     "COUPLING and ENDCOUPLING statements, then in the lines of ",
                     "the form: \"MOL1 \${mol1label}\" and \"MOL2 \${mol2label}\", ",
                     "it is customary that mol1label refer to a molecular label that ",
                     "comes before mol2label in the same file; all data matrices will ",
                     "therefore be transposed"));
    end

    return mol1, mol2, flipmol;

end

function extractGeneralCouplingData(couplingstrarray::Array{<:AbstractString, 1},
                                    mollistfilename::AbstractString, T::DataType=Float64)

    k0 = extractCouplingDataMat("k0", couplingstrarray, mollistfilename, T);
    a0 = extractCouplingDataMat("a0", couplingstrarray, mollistfilename, T);
    Einf = extractCouplingDataMat("Einf", couplingstrarray, mollistfilename, T);
    return k0, a0, Einf;

end

function extractCouplingDataMat(datalabel::AbstractString,
                                couplingstrarray::Array{<:AbstractString, 1},
                                mollistfilename::AbstractString, T::DataType=Float64)

    curridx = findfirst(map(u -> startswith(u, string(datalabel, " ")), couplingstrarray));
    if (curridx == nothing)
        error(string("If ", mollistfilename, " has at least one consecutive pair of ",
                     "COUPLING and ENDCOUPLING statements, then it must contain a line of ",
                     "the form: ", datalabel, "\${", datalabel, "filename}"));
    end
    return readdlm(replace(couplingstrarray[curridx], string(datalabel, " ") => ""), ',', T);

end

function assembleDeltaKImolblock!(myDeltaKImolblockdata::DeltaKIMolBlockData{FT},
                                  myAllMolData::MolSystem{OneMol_WithPhonons{FT}},
                                  myTransData::TransData{FT},
                                  myPeriodicData::PeriodicData{FT}) where FT<:AbstractFloat

    mymol1 = myDeltaKImolblockdata.mol1;
    mymol2 = myDeltaKImolblockdata.mol2;

    mynumblocks1 = myDeltaKImolblockdata.numblocks1;
    mynumblocks2 = myDeltaKImolblockdata.numblocks2;
    myblocklist1 = myDeltaKImolblockdata.blocklist1;
    myblocklist2 = myDeltaKImolblockdata.blocklist2;

    mynumatoms1 = myAllMolData.numatomslist[mymol1];
    mynumatoms2 = myAllMolData.numatomslist[mymol2];

    fill!(myDeltaKImolblockdata.DeltaKImolblock, zero(FT));
    atomposqq = zeros(FT, 3);
    atompospp = zeros(FT, 3);
    mylatticevec = zeros(FT, 3);

    for nn2=1:mynumblocks2
        for nn1=1:mynumblocks1

            mylatticevec .= myblocklist1[nn1] .* myPeriodicData.latticevecs[:,1] .+
            myblocklist2[nn2] .* myPeriodicData.latticevecs[nn2];

            A = view(myDeltaKImolblockdata.DeltaKImolblock,
                     (nn1-1)*mynumatoms1+1:nn1*mynumatoms1,
                     (nn2-1)*mynumatoms2+1:nn2*mynumatoms2);

            k0 = view(myDeltaKImolblockdata.k0molblock,
                      (nn1-1)*mynumatoms1+1:nn1*mynumatoms1,
                      (nn2-1)*mynumatoms2+1:nn2*mynumatoms2);

            a0 = view(myDeltaKImolblockdata.a0molblock,
                      (nn1-1)*mynumatoms1+1:nn1*mynumatoms1,
                      (nn2-1)*mynumatoms2+1:nn2*mynumatoms2);

            Einf = view(myDeltaKImolblockdata.Einfmolblock,
                        (nn1-1)*mynumatoms1+1:nn1*mynumatoms1,
                        (nn2-1)*mynumatoms2+1:nn2*mynumatoms2);

            for qq=1:mynumatoms2

                atomposqq .= myTransData.allAtomPos[:, qq];
                if (mymol2 > 1)
                    atomposqq .= myTransData.allAtomPos[:, myAllMolData.cumnumatomslist[mymol2-1]+qq];
                end

                for pp=1:mynumatoms1

                    atompospp .= myTransData.allAtomPos[:, pp];
                    if (mymol1 > 1)
                        atompospp .= myTransData.allAtomPos[:, myAllMolData.cumnumatomslist[mymol1-1]+pp];
                    end
                    atompospp .+= mylatticevec;
                    A[pp, qq] = keffMorse(norm(atompospp .- atomposqq),
                                          k0[pp, qq], a0[pp, qq], Einf[pp, qq]);
                end
            end
        end
    end

end

function assembleDeltaKIdiagblocks!(myDeltaKIsystemdata::DeltaKISystemData{FT},
                                    myAllMolData::MolSystem{OneMol_WithPhonons{FT}}) where FT<:AbstractFloat

    for myDeltaKImolblockdata in myDeltaKIsystemdata.DeltaKImolblockdataarray

        mymol1 = myDeltaKImolblockdata.mol1;
        mymol2 = myDeltaKImolblockdata.mol2;

        mynumblocks1 = myDeltaKImolblockdata.numblocks1;
        mynumblocks2 = myDeltaKImolblockdata.numblocks2;
        myblocklist1 = myDeltaKImolblockdata.blocklist1;
        myblocklist2 = myDeltaKImolblockdata.blocklist2;

        mynumatoms1 = myAllMolData.numatomslist[mymol1];
        mynumatoms2 = myAllMolData.numatomslist[mymol2];

        currDeltaKIdiagelementsblock = view(myDeltaKIsystemdata.DeltaKIdiagelements,
                                            1:myAllMolData.cumnumatomslist[mymol2]);
        if (mymol2 > 1)
            currDeltaKIdiagelementsblock = view(myDeltaKIsystemdata.DeltaKIdiagelements,
                                                myAllMolData.cumnumatomslist[mymol2-1]+1:myAllMolData.cumnumatomslist[mymol2]);
        end

        fill!(currDeltaKIdiagelementsblock, zero(FT));

        for nn2=1:mynumblocks2

            for nn1=1:mynumblocks1

                A = view(myDeltaKImolblockdata.DeltaKImolblock,
                         (nn1-1)*mynumatoms1+1:nn1*mynumatoms1,
                         (nn2-1)*mynumatoms2+1:nn2*mynumatoms2);

                for qq=1:mynumatoms2

                    for pp=1:mynumatoms1

                        currDeltaKIdiagelementsblock[qq] -= A[pp, qq];

                    end

                end

            end

        end

        currDeltaKIdiagelementsblock = view(myDeltaKIsystemdata.DeltaKIdiagelements,
                                            1:myAllMolData.cumnumatomslist[mymol1]);
        if (mymol1 > 1)
            currDeltaKIdiagelementsblock = view(myDeltaKIsystemdata.DeltaKIdiagelements,
                                                myAllMolData.cumnumatomslist[mymol1-1]+1:myAllMolData.cumnumatomslist[mymol1]);
        end

        fill!(currDeltaKIdiagelementsblock, zero(FT));

        for nn2=1:mynumblocks2

            for nn1=1:mynumblocks1

                A = view(myDeltaKImolblockdata.DeltaKImolblock,
                         (nn1-1)*mynumatoms1+1:nn1*mynumatoms1,
                         (nn2-1)*mynumatoms2+1:nn2*mynumatoms2);

                for qq=1:mynumatoms2

                    for pp=1:mynumatoms1

                        currDeltaKIdiagelementsblock[pp] -= A[pp, qq];

                    end

                end

            end

        end

    end

end

"""

    assembleDeltaKI!(DeltaKI::AbstractArray{Complex{FT}, 2},
                     myDeltaKIsystemdata::DeltaKISystemData{FT},
                     myAllMolData::MolSystem{OneMol_WithPhonons{FT}},
                     myPeriodicData::PeriodicData{FT},
                     blochk::AbstractArray{FT, 1}=zeros(FT, 3)) where FT<:AbstractFloat

Assemble the internuclear spring constant matrix of the collection of molecules (which, at
any wavevector `blochk`, will be a Hermitian positive-definite matrix) for the whole system,
including the diagonal contributions, and stamp this into `DeltaKI`.

"""
function assembleDeltaKI!(DeltaKI::AbstractArray{Complex{FT}, 2},
                          myDeltaKIsystemdata::DeltaKISystemData{FT},
                          myAllMolData::MolSystem{OneMol_WithPhonons{FT}},
                          myPeriodicData::PeriodicData{FT},
                          blochk::AbstractArray{FT, 1}=zeros(FT, 3)) where FT<:AbstractFloat

    mylatticevec = zeros(FT, 3);
    fill!(DeltaKI, zero(eltype(DeltaKI)));

    # first, do the off-diagonal blocks
    for myDeltaKImolblockdata in myDeltaKIsystemdata.DeltaKImolblockdataarray

        mymol1 = myDeltaKImolblockdata.mol1;
        mymol2 = myDeltaKImolblockdata.mol2;

        if (mymol2 > mymol1) # only assemble the upper-right part, then use Hermiticity for the rest
            startidx1 = 1;
            if (mymol1 > 1)
                startidx1 = 3*myAllMolData.cumnumatomslist[mymol1-1] + 1;
            end
            startidx2 = 1;
            if (mymol2 > 1)
                startidx2 = 3*myAllMolData.cumnumatomslist[mymol2-1] + 1;
            end

            currKIview = view(DeltaKI, startidx1:3*myAllMolData.cumnumatomslist[mymol1],
                              startidx2:3*myAllMolData.cumnumatomslist[mymol2]);
            fill!(currKIview, zero(eltype(currKIview)));

            mynumblocks1 = myDeltaKImolblockdata.numblocks1;
            mynumblocks2 = myDeltaKImolblockdata.numblocks2;
            myblocklist1 = myDeltaKImolblockdata.blocklist1;
            myblocklist2 = myDeltaKImolblockdata.blocklist2;

            mynumatoms1 = myAllMolData.numatomslist[mymol1];
            mynumatoms2 = myAllMolData.numatomslist[mymol2];

            for nn2=1:mynumblocks2

                for nn1=1:mynumblocks1

                    mylatticevec .= myblocklist1[nn1] .* myPeriodicData.latticevecs[:,1] .+
                    myblocklist2[nn2] .* myPeriodicData.latticevecs[nn2];
                    phaseterm = exp(-1im*dot(blochk, mylatticevec));
                    A = view(myDeltaKImolblockdata.DeltaKImolblock,
                             (nn1-1)*mynumatoms1+1:nn1*mynumatoms1,
                             (nn2-1)*mynumatoms2+1:nn2*mynumatoms2);

                    for qq=1:mynumatoms2

                        for pp=1:mynumatoms1

                            currKIview[3*pp-2:3*pp, 3*qq-2:3*qq] .+=
                            Matrix{Complex{eltype(A)}}(phaseterm * A[pp, qq] * I, 3, 3);

                        end

                    end

                end

            end

        end

    end

    # now make it Hermitian, and multiply by 2 to cancel the division by 2
    hermpart!(DeltaKI);
    lmul!(2, DeltaKI);

    # now add the diagonal blocks, which should have no dependence on wavevector
    for qq=1:myAllMolData.cumnumatomslist[end]

        DeltaKI[3*qq-2:3*qq, 3*qq-2:3*qq] .+=
        Matrix{Complex{eltype(myDeltaKIsystemdata.DeltaKIdiagelements)}}(myDeltaKIsystemdata.DeltaKIdiagelements[qq] * I, 3, 3);

    end

    # make it Hermitian again, just to be sure
    hermpart!(DeltaKI);

end


text:
using Gurobi

include(joinpath(Pkg.dir("MathProgBase"),"test","linprog.jl"))
linprogtest(GurobiSolver(InfUnbdInfo=1))

include(joinpath(Pkg.dir("MathProgBase"),"test","linproginterface.jl"))
linprogsolvertest(GurobiSolver())

include(joinpath(Pkg.dir("MathProgBase"),"test","mixintprog.jl"))
mixintprogtest(GurobiSolver())

include(joinpath(Pkg.dir("MathProgBase"),"test","quadprog.jl"))
quadprogtest(GurobiSolver())
socptest(GurobiSolver())
qpdualtest(GurobiSolver(QCPDual=1))

include(joinpath(Pkg.dir("MathProgBase"),"test","conicinterface.jl"))
coniclineartest(GurobiSolver())
conicSOCtest(GurobiSolver())
conicSOCRotatedtest(GurobiSolver())
conicSOCINTtest(GurobiSolver())


text:
using SparseArrays
import SparseArrays: dropzeros!

"""
dropzeros!(sparseMatrix, drop_criteria)
# Description
- Drops machine zeros in sparse matrix
# Arguments
- `A`: a sparse matrix
- `drop_criteria`: criteria for dropping entries
# return
- nothing
"""
function dropzeros!(A, drop_criteria)
    i, j = findnz(A)
    for loop in 1:length(i)
        if abs(A[i[loop], j[loop]]) < drop_criteria
            A[i[loop], j[loop]] = 0.0
        end
    end
    dropzeros!(A)
end

"""
droprelativezeros!(A; scale=100)
# Description
- Drops machine zeros in sparse matrix
# Arguments
- `A`: a sparse matrix
- `drop_criteria`: criteria for dropping entries
# return
- nothing
# This function calls
dropzeros!(A, eps(scale * maximum(da)))
"""
function droprelativezeros!(A; scale = 100)
    dropzeros!(A, eps(scale * maximum(abs.(A))))
end

max_stars_repo_name: sandreza/StatisticalNonlocality

max_issues_repo_name: sandreza/StatisticalNonlocality

max_forks_repo_name: sandreza/StatisticalNonlocality


text:
#script that estimates the dir bin gas network model on emid data and evaluates GAS
# forecasting performances

using Utilities,AReg,StaticNets,JLD,MLBase,StatsBase
using PyCall; pygui(:qt); using PyPlot
using DynNets
#estimate and save for half the datase (after LTRO) or whole data?
halfPeriod = true

## Load data
fold_Path =  "/home/Domenico/Dropbox/Dynamic_Networks/data/emid_data/juliaFiles/"
loadFilePartialName = "Weekly_eMid_Data_from_"
halfPeriod ? periodEndStr =  "2012_03_12_to_2015_02_27.jld" : periodEndStr =  "2009_06_22_to_2015_02_27.jld"
@load(fold_Path*loadFilePartialName*periodEndStr, AeMidWeekly_T,banksIDs,inactiveBanks, YeMidWeekly_T,weekInd,datesONeMid ,degsIO_T,strIO_T)

#
# ## Binary  estimates
# divide the data in Train and test samples and estimate the gas model on the train
N2,T = size(degsIO_T);N = round(Int,N2/2)
Ttrain = round(Int,T/2)
Ttest = T - Ttrain
modGasDirBin1_eMid = DynNets.SdErgmDirBin1(degsIO_T)
modGasDirBin1_eMidTrain = DynNets.SdErgmDirBin1(degsIO_T[:,1:Ttrain])
SSestForTargTrain = DynNets.estimateSnapSeq(modGasDirBin1_eMidTrain)
# ## test targeted estimate

estTargDirBin1_eMidTrain,~ = DynNets.estimateTarg(modGasDirBin1_eMidTrain;SSest = SSestForTargTrain)
#estTargDirBin1_eMid,~ = DynNets.estimateTarg(modGasDirBin1_eMid)

# ## save estimates
#@save(fold_Path*"EstDirBin1Gas"*loadFilePartialName*periodEndStr,estTargDirBin1_eMidTrain,modGasDirBin1_eMidTrain,SSestForTargTrain)
##load the estimates instead of re estimating every time
#@load(fold_Path*"EstDirBin1Gas"*loadFilePartialName*periodEndStr,estTargDirBin1_eMidTrain,SSestForTargTrain)

## Forecast functions

N2,T = size(degsIO_T);N = round(Int,N2/2)

t=1

ftotIO_t_Start = identify!(modGasDirBin1_eMid,gasPar[1].*(1-gasPar[2]))
tvPartp1 = predict_score_driven_par( modGasDirBin1_eMid,N,degsIO_T[:,t],
                         copy(ftotIO_t_Start),zeros(N2,N2),trues(N2),
                         identify!(modGasDirBin1_eMid,gasPar[1]),gasPar[2],gasPar[3])[1]

## --------------Arrivato qui al test della funzione che valuta il foreasting
gasParEstOnTrain = estTargDirBin1_eMidTrain
testObsNet_T = AeMidWeekly_T[:,:,Ttrain+1:end]

foreEval = forecastEvalGasNetDirBin1(gasParEstOnTrain,testObsNet_T)
rocCurve(foreEval[1],foreEval[2])

#check that forecasted fitnesses and degrees are reasonable
foreFit,~ = score_driven_filter_or_dgp( DynNets.SdErgmDirBin1(degsIO_T),[gasParEstOnTrain[1];gasParEstOnTrain[2];gasParEstOnTrain[3]])
plot(foreFit')
expDegsIO_Ttest = zeros(N2,Ttest)
for t=1:Ttest
    expDegsIO_Ttest[:,t] = [sumSq(StaticNets.expMatrix(StaticNets.fooErgmDirBin1,foreEval[3][:,t]),2);sumSq(StaticNets.expMatrix(StaticNets.fooErgmDirBin1,foreEval[3][:,t]),1)]
end
plot(expDegsIO_Ttest' )

## Wrap everithing up in a function

#

#

max_stars_repo_name: domenicodigangi/SDERGM

max_issues_repo_name: domenicodigangi/SDERGM

max_forks_repo_name: domenicodigangi/SDERGM


text:
module Mimi

using DataFrames
using DataStructures
using Distributions
using Electron
using JSON
using NamedArrays
using StringBuilders

export
    @defcomp,
    MarginalModel,
    Model,
    add_comp!,
    components,
    connect_param!,
    create_marginal_model,
    disconnect_param!,
    explore,
    getdataframe,
    gettime,
    get_param_value,
    get_var_value,
    hasvalue,
    is_first,
    is_last,
    is_time,
    is_timestep,
    load_comps,
    modeldef,
    name,
    new_comp,
    parameters,
    plot_comp_graph,
    replace_comp!,
    set_dimension!,
    set_leftover_params!,
    set_param!,
    update_param!,
    update_params!,
    variables

include("core/types.jl")

# After loading types, the rest can just be alphabetical
include("core/build.jl")
include("core/connections.jl")
include("core/defs.jl")
include("core/defcomp.jl")
include("core/dimensions.jl")
include("core/instances.jl")
include("core/references.jl")
include("core/time.jl")
include("core/model.jl")
include("explorer/explore.jl")
include("mcs/mcs.jl")
include("utils/graph.jl")
include("utils/plotting.jl")
include("utils/getdataframe.jl")
include("utils/lint_helper.jl")
include("utils/misc.jl")

"""
    load_comps(dirname::String="./components")

Call include() on all the files in the indicated directory `dirname`.
This avoids having modelers create a long list of include()
statements. Just put all the components in a directory.
"""
function load_comps(dirname::String="./components")
    files = readdir(dirname)
    for file in files
        if endswith(file, ".jl")
            pathname = joinpath(dirname, file)
            include(pathname)
        end
    end
end

# Components are defined here to allow pre-compilation to work
function __init__()
    compdir = joinpath(dirname(@__FILE__), "components")
    load_comps(compdir)
end

end # module


text:
using Test, Random, Pkg, Registrator

tmpdir = joinpath(tempdir(), Random.randstring())
mkdir(tmpdir)
println(tmpdir)

try
    cd(tmpdir)

    # Create an empty registry
    Registrator.create_registry("TestRegistry", "not_a_real_repo")
    @test isfile(joinpath(tmpdir, "TestRegistry", "Registry.toml"))

    # Create a package to register
    Pkg.generate("TestPackage")
    pkgpath = joinpath(tmpdir, "TestPackage")
    run(`git -C $pkgpath init`)
    run(`git -C $pkgpath add --all`)
    run(`git -C $pkgpath commit -m 'initial commit'`)
    Pkg.activate("TestPackage")
    using TestPackage

    Registrator.register(TestPackage, "TestRegistry"; repo="not_a_real_repo")
finally
    rm(tmpdir, recursive=true)
end


text:
"""
This is borrowed from DataFrames.jl. Suggest up to 8 names in which differ little
from what the user typed in Levenshtein distance.
"""
function fuzzymatch(colnames, name::Symbol)
    ucname = uppercase(string(name))
    dist = [(levenshtein(uppercase(string(x)), ucname), x) for x in colnames]
    sort!(dist)
    c = [count(x -> x[1] <= i, dist) for i in 0:2]
    maxd = max(0, searchsortedlast(c, 8) - 1)
    return [s for (d, s) in dist if d <= maxd]
end

colnames(table) = string.(Tables.columnnames(table))

vecmap(f, iter) = [f(el) for el in iter]

function data_options(t::Observable; keywords=[""], suffix="")
    return lift(t) do table
        names = map(name -> name * suffix, colnames(table))
        return [keyword => names for keyword in keywords]
    end
end

iscontinuous(v::AbstractVector) = false
iscontinuous(v::AbstractVector{<:Number}) = true
iscontinuous(v::AbstractVector{<:Bool}) = false

function indexoftype(::Type{T}, list, el) where {T}
    idx = 0
    for x in list
        idx += x isa T
        isequal(x, el) && return idx
    end
    return nothing
end

macro maybereturn(x)
    expr = quote
        local tmp = $(esc(x))
        isnothing(tmp) && return nothing
        tmp
    end
    return expr
end

function buttonclass(positive)
    class = "text-xl font-semibold rounded text-left py-2 px-4 bg-opacity-75"
    class *= positive ? " bg-blue-100 hover:bg-blue-200 text-blue-800 hover:text-blue-900 mr-8" :
        " bg-red-100 hover:bg-red-200 text-red-800 hover:text-red-900"
    return class
end

struct Call
    fs::Vector{String}
    positional::Vector{String}
    named::Vector{Pair{String, String}}
end

Call() = Call(String[], String[], Pair{String, String}[])

function compute_calls(str::AbstractString)
    calls, call = Call[], Call()
    for chunk in split(str, ' ')
        isempty(chunk) && continue
        s = split(chunk, ':')
        if length(s) == 1
            if only(s) == "+"
                push!(calls, call)
                call = Call()
            else
                push!(call.fs, only(s))
            end
        else
            pre, post = s
            positional, named = call.positional, call.named
            isempty(pre) ? push!(positional, post) : push!(named, pre => post)
        end
    end
    push!(calls, call)
    return calls
end

for sym in (:on, :onany)
    trysim = Symbol(:try, sym)
    @eval function $trysim(f, session::Session, obs::Observable...)
        error_msg = Observable("")
        onjs(session, error_msg, js"""
            function (value) {
                alert(value);
            }
        """)
        return $sym(session, obs...) do args...
            try
                f(args...)
            catch err
                io = IOBuffer()
                print(io, "Could not complete command due to the following error.")
                print(io, "\n\n")
                print(io, err)
                error_msg[] = String(take!(io))
            end
        end
    end
end

function move_item(v, (old, new))
    return map(1:length(v)) do i
        i == new && return v[old]
        old ≤ i ≤ new && return v[i+1]
        old ≥ i ≥ new && return v[i-1]
        return v[i]
    end
end

function remove_item(v, idx)
    return map(1:length(v)-1) do i
        return i < idx ? v[i] : v[i+1]
    end
end

function insert_item(v, idx, value)
    return map(1:length(v)+1) do i
        i < idx && return v[i]
        i > idx && return v[i-1]
        return value
    end
end

function scrollable_component(args...; kwargs...)
    return DOM.div(
        DOM.div(args...; class="absolute left-0 right-8");
        class="overflow-y-scroll h-full w-full relative",
        kwargs...
    )
end


text:
# runtests.jl

using AppliMaster

using Test

using AppliAR, AppliSales, AppliGeneralLedger

using Query

using DataFrames

using Dates

@testset "Test AppliSales" begin
    orders = AppliSales.process()
    @test length(orders) == 3
    @test length(orders[1].students) == 1
    @test length(orders[2].students) == 2
    @test length(orders[3].students) == 1
    @test orders[1].training.price == 1000
end

@testset " Test AppliAR - unpaid invoices" begin
    orders = AppliSales.process()
    AppliAR.process(orders)
    unpaid_invoices = retrieve_unpaid_invoices()
    @test length(unpaid_invoices) == 3
    @test id(unpaid_invoices[1]) == "A1001"
    cmd = `rm test_invoicing.txt invoicenbr.txt`
    run(cmd)
end

@testset "Test AppliAR - entries unpaid invoices" begin
    orders = AppliSales.process()
    entries = AppliAR.process(orders)
    @test length(entries) == 3
    @test entries[1].from == 1300
    @test entries[1].to == 8000
    @test entries[1].debit == 1000
    @test entries[1].credit == 0
    @test entries[1].vat == 210.0
    cmd = `rm test_invoicing.txt invoicenbr.txt`
    run(cmd)
end

@testset "Test GeneralLedger - accounts receivable, bank, vat, sales" begin
    orders = AppliSales.process()

    journal_entries_unpaid_invoices = AppliAR.process(orders)
    AppliGeneralLedger.process(journal_entries_unpaid_invoices)

    unpaid_invoices = AppliAR.retrieve_unpaid_invoices()
    stm1 = BankStatement(Date(2020-01-15), "Duck City Chronicals Invoice A1002", "NL39INGB", 2420.0)
    stm2 = BankStatement(Date(2020-01-15), "Donalds Hardware Store Bill A1003", "NL39INGB", 1210.0)
    stms = [stm1, stm2]

    journal_entries_paid_invoices = AppliAR.process(unpaid_invoices, stms)
    AppliGeneralLedger.process(journal_entries_paid_invoices)

    df = DataFrame(AppliGeneralLedger.read_from_file("./test_ledger.txt"))

    df2 = df |> @filter(_.accountid == 1300) |> DataFrame
    @test sum(df2.debit - df2.credit) == 1210

    df2 = df |> @filter(_.accountid == 1150) |> DataFrame # bank
    @test sum(df2.debit - df2.credit) == 3630

    df2 = df |> @filter(_.accountid == 4000) |> DataFrame # vat
    @test sum(df2.credit - df2.debit) == 840

    df2 = df |> @filter(_.accountid == 8000) |> DataFrame # sales
    @test sum(df2.credit - df2.debit) == 4000

    @test sum(df.debit - df.credit) == 0.0

    cmd = `rm test_invoicing.txt test_invoicing_paid.txt test_journal.txt test_ledger.txt invoicenbr.txt`
    run(cmd)
end


text:
@isdefined(plotresults) || (plotresults = false)

if plotresults

    postproc || error("Cannot plot results without going through post-processing.")
    @eval begin
        using Plots
        using LinearAlgebra

        p1 = scatter(Θ, real.(norm.(ffd)))
        p2 = heatmap(clamp.(real.(norm.(nfd)), 0.0, 2.0))
        p3 = contour(clamp.(real.(norm.(nfd)), 0.0, 2.0))
        plot(p1,p2,p3,layout=(3,1))
    end
end


text:
# This structure contains necessary information for bundle methods.
mutable struct BundleModel
    n::Int              # dimension in the solution space
    N::Int              # number of separable functions in the objective
    ncuts_per_iter::Int # number of cuts per iteration
    cut_indices::Dict{Int,Vector{Int}}
    model::JuMP.Model   # Bundle model

    # user-defined function to evaluate f
    # and return the value and its subgradients
    evaluate_f

    time::Vector{Float64}

    user_data

    function BundleModel(n::Int=0, N::Int=0, ncuts_per_iter::Int=0, func=nothing)
        bundle = new()
        bundle.n = n
        bundle.N = N
        bundle.ncuts_per_iter = ncuts_per_iter
        bundle.cut_indices = Dict(i => [j for j = i:ncuts_per_iter:N] for i = 1:ncuts_per_iter)
        bundle.model = JuMP.Model()
        # bundle.splitvars = splitvars
        bundle.evaluate_f = func
        bundle.time = []
        bundle.user_data = nothing
        return bundle
    end
end

# This returns BundleModel object.
get_model(method::AbstractMethod)::BundleModel = BundleModel()

get_model(model::BundleModel)::JuMP.Model = model.model

set_optimizer(model::BundleModel, optimizer) = JuMP.set_optimizer(get_model(model), optimizer)

function set_user_data(model::BundleModel, user_data)
    model.user_data = user_data
end

function solve_model!(model::BundleModel)
    m = get_model(model)
    @assert !isnothing(m.moi_backend.optimizer)
    JuMP.optimize!(m)
    push!(model.time, JuMP.solve_time(m))
end


text:
# Functions

import MutableArithmetics

"""
    AbstractFunction

Abstract supertype for function objects.
"""
abstract type AbstractFunction <: MutableArithmetics.AbstractMutable end

"""
    output_dimension(f::AbstractFunction)

Return 1 if `f` has a scalar output and the number of output components if `f`
has a vector output.
"""
function output_dimension end

"""
    AbstractScalarFunction

Abstract supertype for scalar-valued function objects.
"""
abstract type AbstractScalarFunction <: AbstractFunction end
output_dimension(::AbstractScalarFunction) = 1

Base.broadcastable(f::AbstractScalarFunction) = Ref(f)
Base.ndims(::Type{<:AbstractScalarFunction}) = 0
Base.ndims(::AbstractScalarFunction) = 0

"""
    AbstractVectorFunction

Abstract supertype for vector-valued function objects.
"""
abstract type AbstractVectorFunction <: AbstractFunction end

"""
    SingleVariable(variable)

The function that extracts the scalar variable referenced by `variable`, a `VariableIndex`.
This function is naturally be used for single variable bounds or integrality constraints.
"""
struct SingleVariable <: AbstractScalarFunction
    variable::VariableIndex
end

"""
    VectorOfVariables(variables)

The function that extracts the vector of variables referenced by `variables`, a `Vector{VariableIndex}`.
This function is naturally be used for constraints that apply to groups of variables, such as an "all different" constraint, an indicator constraint, or a complementarity constraint.
"""
struct VectorOfVariables <: AbstractVectorFunction
    variables::Vector{VariableIndex}
end
output_dimension(f::VectorOfVariables) = length(f.variables)

"""
    struct ScalarAffineTerm{T}
        coefficient::T
        variable_index::VariableIndex
    end

Represents ``c x_i`` where ``c`` is `coefficient` and ``x_i`` is the variable
identified by `variable_index`.
"""
struct ScalarAffineTerm{T}
    coefficient::T
    variable_index::VariableIndex
end

# Note: ScalarAffineFunction is mutable because its `constant` field is likely of an immutable
# type, while its `terms` field is of a mutable type, meaning that creating a `ScalarAffineFunction`
# allocates, and it is desirable to provide a zero-allocation option for working with
# ScalarAffineFunctions. See https://github.com/jump-dev/MathOptInterface.jl/pull/343.
"""
    ScalarAffineFunction{T}(terms, constant)

The scalar-valued affine function ``a^T x + b``, where:
* ``a`` is a sparse vector specified by a list of
  [`ScalarAffineTerm`](@ref) structs.
* ``b`` is a scalar specified by `constant::T`

Duplicate variable indices in `terms` are accepted, and the corresponding
coefficients are summed together.
"""
mutable struct ScalarAffineFunction{T} <: AbstractScalarFunction
    terms::Vector{ScalarAffineTerm{T}}
    constant::T
end

"""
    struct VectorAffineTerm{T}
        output_index::Int64
        scalar_term::ScalarAffineTerm{T}
    end

A `ScalarAffineTerm` plus its index of the output component of a
`VectorAffineFunction` or `VectorQuadraticFunction`.
`output_index` can also be interpreted as a row index into a sparse matrix,
where the `scalar_term` contains the column index and coefficient.
"""
struct VectorAffineTerm{T}
    output_index::Int64
    scalar_term::ScalarAffineTerm{T}
end

function VectorAffineTerm(
    output_index::Base.Integer,
    scalar_term::ScalarAffineTerm,
)
    return VectorAffineTerm(convert(Int64, output_index), scalar_term)
end

"""
    VectorAffineFunction{T}(terms, constants)

The vector-valued affine function ``A x + b``, where:
* ``A`` is a sparse matrix specified by a list of `VectorAffineTerm` objects.
* ``b`` is a vector specified by `constants`

Duplicate indices in the ``A`` are accepted, and the corresponding coefficients
are summed together.
"""
struct VectorAffineFunction{T} <: AbstractVectorFunction
    terms::Vector{VectorAffineTerm{T}}
    constants::Vector{T}
end
output_dimension(f::VectorAffineFunction) = length(f.constants)

"""
    struct ScalarQuadraticTerm{T}
        coefficient::T
        variable_index_1::VariableIndex
        variable_index_2::VariableIndex
    end

Represents ``c x_i x_j`` where ``c`` is `coefficient`, ``x_i`` is the variable
identified by `variable_index_1` and ``x_j`` is the variable identified by
`variable_index_2`.
"""
struct ScalarQuadraticTerm{T}
    coefficient::T
    variable_index_1::VariableIndex
    variable_index_2::VariableIndex
end

# Note: ScalarQuadraticFunction is mutable because its `constant` field is likely of an immutable
# type, while its other fields are of mutable types, meaning that creating a `ScalarQuadraticFunction`
# allocates, and it is desirable to provide a zero-allocation option for working with
# ScalarQuadraticFunctions. See https://github.com/jump-dev/MathOptInterface.jl/pull/343.
"""
    ScalarQuadraticFunction{T}(affine_terms, quadratic_terms, constant)

The scalar-valued quadratic function ``\\frac{1}{2}x^TQx + a^T x + b``, where:
* ``a`` is a sparse vector specified by a list of `ScalarAffineTerm` structs.
* ``b`` is a scalar specified by `constant`.
* ``Q`` is a symmetric matrix specified by a list of `ScalarQuadraticTerm`
  structs.

Duplicate indices in ``a`` or ``Q`` are accepted, and the corresponding
coefficients are summed together. "Mirrored" indices `(q,r)` and `(r,q)` (where
`r` and `q` are `VariableIndex`es) are considered duplicates; only one need be
specified.

For example, for two scalar variables ``y, z``, the quadratic expression
``yz + y^2`` is represented by the terms
`ScalarQuadraticTerm.([1.0, 2.0], [y, y], [z, y])`.
"""
mutable struct ScalarQuadraticFunction{T} <: AbstractScalarFunction
    affine_terms::Vector{ScalarAffineTerm{T}}
    quadratic_terms::Vector{ScalarQuadraticTerm{T}}
    constant::T
end

"""
    struct VectorQuadraticTerm{T}
        output_index::Int64
        scalar_term::ScalarQuadraticTerm{T}
    end

A [`ScalarQuadraticTerm`](@ref) plus its index of the output component of a
`VectorQuadraticFunction`. Each output component corresponds to a
distinct sparse matrix ``Q_i``.
"""
struct VectorQuadraticTerm{T}
    output_index::Int64
    scalar_term::ScalarQuadraticTerm{T}
end

function VectorQuadraticTerm(
    output_index::Base.Integer,
    scalar_term::ScalarQuadraticTerm,
)
    return VectorQuadraticTerm(convert(Int64, output_index), scalar_term)
end

"""
    VectorQuadraticFunction{T}(affine_terms, quadratic_terms, constants)

The vector-valued quadratic function with i`th` component ("output index")
defined as ``\\frac{1}{2}x^TQ_ix + a_i^T x + b_i``, where:
* ``a_i`` is a sparse vector specified by the `VectorAffineTerm`s with
  `output_index == i`.
* ``b_i`` is a scalar specified by `constants[i]`
* ``Q_i`` is a symmetric matrix specified by the `VectorQuadraticTerm` with
  `output_index == i`.

Duplicate indices in ``a_i`` or ``Q_i`` are accepted, and the corresponding
coefficients are summed together. "Mirrored" indices `(q,r)` and `(r,q)` (where
`r` and `q` are `VariableIndex`es) are considered duplicates; only one need be
specified.
"""
struct VectorQuadraticFunction{T} <: AbstractVectorFunction
    affine_terms::Vector{VectorAffineTerm{T}}
    quadratic_terms::Vector{VectorQuadraticTerm{T}}
    constants::Vector{T}
end
output_dimension(f::VectorQuadraticFunction) = length(f.constants)

# Function modifications

"""
    AbstractFunctionModification

An abstract supertype for structs which specify partial modifications to functions, to be used for making small modifications instead of replacing the functions entirely.
"""
abstract type AbstractFunctionModification end

"""
    ScalarConstantChange{T}(new_constant::T)

A struct used to request a change in the constant term of a scalar-valued function.
Applicable to `ScalarAffineFunction` and `ScalarQuadraticFunction`.
"""
struct ScalarConstantChange{T} <: AbstractFunctionModification
    new_constant::T
end

"""
    VectorConstantChange{T}(new_constant::Vector{T})

A struct used to request a change in the constant vector of a vector-valued function.
Applicable to `VectorAffineFunction` and `VectorQuadraticFunction`.
"""
struct VectorConstantChange{T} <: AbstractFunctionModification
    new_constant::Vector{T}
end

"""
    ScalarCoefficientChange{T}(variable::VariableIndex, new_coefficient::T)

A struct used to request a change in the linear coefficient of a single variable
in a scalar-valued function.
Applicable to `ScalarAffineFunction` and `ScalarQuadraticFunction`.
"""
struct ScalarCoefficientChange{T} <: AbstractFunctionModification
    variable::VariableIndex
    new_coefficient::T
end

# Note: MultiRowChange is mutable because its `variable` field of an immutable
# type, while `new_coefficients` is of a mutable type, meaning that creating a `MultiRowChange`
# allocates, and it is desirable to provide a zero-allocation option for working with
# MultiRowChanges. See https://github.com/jump-dev/MathOptInterface.jl/pull/343.
"""
    MultirowChange{T}(variable::VariableIndex, new_coefficients::Vector{Tuple{Int64, T}})

A struct used to request a change in the linear coefficients of a single
variable in a vector-valued function. New coefficients are specified by
`(output_index, coefficient)` tuples. Applicable to `VectorAffineFunction` and
`VectorQuadraticFunction`.
"""
mutable struct MultirowChange{T} <: AbstractFunctionModification
    variable::VariableIndex
    new_coefficients::Vector{Tuple{Int64,T}}
end

function MultirowChange(
    variable::VariableIndex,
    new_coefficients::Vector{Tuple{Ti,T}},
) where {Ti<:Base.Integer,T}
    return MultirowChange(
        variable,
        [(convert(Int64, i), j) for (i, j) in new_coefficients],
    )
end

# Implementation of comparison for functions
function Base.:(==)(f::VectorOfVariables, g::VectorOfVariables)
    return f.variables == g.variables
end

function Base.isapprox(
    f::Union{SingleVariable,VectorOfVariables},
    g::Union{SingleVariable,VectorOfVariables};
    kwargs...,
)
    return f == g
end

# For affine and quadratic functions, terms are compressed in a dictionary using `_dicts` and then the dictionaries are compared with `dict_compare`
function dict_compare(d1::Dict, d2::Dict{<:Any,T}, compare::Function) where {T}
    return all(kv -> compare(kv.second, Base.get(d2, kv.first, zero(T))), d1)
end

# Build a dictionary where the duplicate keys are summed
function sum_dict(kvs::Vector{Pair{K,V}}) where {K,V}
    d = Dict{K,V}()
    for kv in kvs
        key = kv.first
        d[key] = kv.second + Base.get(d, key, zero(V))
    end
    return d
end

"""
    coefficient(t::Union{ScalarAffineTerm, ScalarQuadraticTerm
                         VectorAffineTerm, VectorQuadraticTerm})

Finds the coefficient stored in the term `t`.
"""
function coefficient end

function coefficient(t::Union{ScalarAffineTerm,ScalarQuadraticTerm})
    return t.coefficient
end
function coefficient(t::Union{VectorAffineTerm,VectorQuadraticTerm})
    return t.scalar_term.coefficient
end

"""
    term_indices(t::Union{ScalarAffineTerm, ScalarQuadraticTerm,
                         VectorAffineTerm, VectorQuadraticTerm})

Returns the indices of the input term `t` as a tuple of `Int`s.

* For `t::ScalarAffineTerm`, this is a 1-tuple of the variable index.
* For `t::ScalarQuadraticTerm`, this is a 2-tuple of the variable indices
  in non-decreasing order.
* For `t::VectorAffineTerm`, this is a 2-tuple of the row/output and
  variable indices.
* For `t::VectorQuadraticTerm`, this is a 3-tuple of the row/output and
  variable indices in non-decreasing order.
"""
term_indices(t::ScalarAffineTerm) = (t.variable_index.value,)
function term_indices(t::ScalarQuadraticTerm)
    return minmax(t.variable_index_1.value, t.variable_index_2.value)
end
function term_indices(t::Union{VectorAffineTerm,VectorQuadraticTerm})
    return (t.output_index, term_indices(t.scalar_term)...)
end

"""
    term_pair(t::Union{ScalarAffineTerm, ScalarQuadraticTerm,
                       VectorAffineTerm, VectorQuadraticTerm})

Returns the pair [`term_indices`](@ref) `=>` [`coefficient`](@ref) of the term.
"""
function term_pair(
    t::Union{
        ScalarAffineTerm,
        ScalarQuadraticTerm,
        VectorAffineTerm,
        VectorQuadraticTerm,
    },
)
    return term_indices(t) => coefficient(t)
end

function _dicts(f::Union{ScalarAffineFunction,VectorAffineFunction})
    return (sum_dict(term_pair.(f.terms)),)
end

function _dicts(f::Union{ScalarQuadraticFunction,VectorQuadraticFunction})
    return (
        sum_dict(term_pair.(f.affine_terms)),
        sum_dict(term_pair.(f.quadratic_terms)),
    )
end

"""
    constant(f::Union{ScalarAffineFunction, ScalarQuadraticFunction})

Returns the constant term of the scalar function
"""
constant(f::Union{ScalarAffineFunction,ScalarQuadraticFunction}) = f.constant

"""
    constant(f::Union{VectorAffineFunction, VectorQuadraticFunction})

Returns the vector of constant terms of the vector function
"""
constant(f::Union{VectorAffineFunction,VectorQuadraticFunction}) = f.constants

function Base.isapprox(
    f::F,
    g::G;
    kwargs...,
) where {
    F<:Union{
        ScalarAffineFunction,
        ScalarQuadraticFunction,
        VectorAffineFunction,
        VectorQuadraticFunction,
    },
    G<:Union{
        ScalarAffineFunction,
        ScalarQuadraticFunction,
        VectorAffineFunction,
        VectorQuadraticFunction,
    },
}
    return isapprox(constant(f), constant(g); kwargs...) && all(
        dict_compare.(
            _dicts(f),
            _dicts(g),
            (α, β) -> isapprox(α, β; kwargs...),
        ),
    )
end

function constant(
    f::Union{ScalarAffineFunction,ScalarQuadraticFunction},
    T::DataType,
)
    return constant(f)
end
function constant(
    f::Union{VectorAffineFunction,VectorQuadraticFunction},
    T::DataType,
)
    return constant(f)
end

"""
    constant(f::SingleVariable, T::DataType)

The constant term of a `SingleVariable` function is
the zero value of the specified type `T`.
"""
constant(f::SingleVariable, T::DataType) = zero(T)

"""
    constant(f::VectorOfVariables, T::DataType)

The constant term of a `VectorOfVariables` function is a
vector of zero values of the specified type `T`.
"""
constant(f::VectorOfVariables, T::DataType) = zeros(T, length(f.variables))

# isbits type, nothing to copy
Base.copy(func::SingleVariable) = func

Base.copy(func::VectorOfVariables) = VectorOfVariables(copy(func.variables))

"""
    copy(func::Union{ScalarAffineFunction, VectorAffineFunction})

Return a new affine function with a shallow copy of the terms and constant(s)
from `func`.
"""
function Base.copy(
    func::F,
) where {F<:Union{ScalarAffineFunction,VectorAffineFunction}}
    return F(copy(func.terms), copy(constant(func)))
end

"""
    copy(func::Union{ScalarQuadraticFunction, VectorQuadraticFunction})

Return a new quadratic function with a shallow copy of the terms and constant(s)
from `func`.
"""
function Base.copy(
    func::F,
) where {F<:Union{ScalarQuadraticFunction,VectorQuadraticFunction}}
    return F(
        copy(func.affine_terms),
        copy(func.quadratic_terms),
        copy(constant(func)),
    )
end

# Define shortcuts for
# SingleVariable -> ScalarAffineFunction
function ScalarAffineFunction{T}(f::SingleVariable) where {T}
    return ScalarAffineFunction([ScalarAffineTerm(one(T), f.variable)], zero(T))
end
# VectorOfVariables -> VectorAffineFunction
function VectorAffineFunction{T}(f::VectorOfVariables) where {T}
    n = length(f.variables)
    return VectorAffineFunction(
        map(
            i -> VectorAffineTerm(i, ScalarAffineTerm(one(T), f.variables[i])),
            1:n,
        ),
        zeros(T, n),
    )
end

# Conversion between scalar functions
# Conversion to SingleVariable
function Base.convert(::Type{SingleVariable}, f::ScalarAffineFunction)
    if (
        !iszero(f.constant) ||
        !isone(length(f.terms)) ||
        !isone(f.terms[1].coefficient)
    )
        throw(InexactError(:convert, SingleVariable, f))
    end
    return SingleVariable(f.terms[1].variable_index)
end

function Base.convert(
    ::Type{SingleVariable},
    f::ScalarQuadraticFunction{T},
) where {T}
    return convert(SingleVariable, convert(ScalarAffineFunction{T}, f))
end

# Conversion to ScalarAffineFunction
function Base.convert(::Type{ScalarAffineFunction{T}}, α::T) where {T}
    return ScalarAffineFunction{T}(ScalarAffineTerm{T}[], α)
end

function Base.convert(
    ::Type{ScalarAffineFunction{T}},
    f::SingleVariable,
) where {T}
    return ScalarAffineFunction{T}(f)
end

function Base.convert(
    ::Type{ScalarAffineTerm{T}},
    t::ScalarAffineTerm{T},
) where {T}
    return t
end

function Base.convert(
    ::Type{ScalarAffineTerm{T}},
    t::ScalarAffineTerm,
) where {T}
    return ScalarAffineTerm{T}(t.coefficient, t.variable_index)
end

function Base.convert(
    ::Type{ScalarAffineFunction{T}},
    f::ScalarAffineFunction{T},
) where {T}
    return f
end

function Base.convert(
    ::Type{ScalarAffineFunction{T}},
    f::ScalarAffineFunction,
) where {T}
    return ScalarAffineFunction{T}(f.terms, f.constant)
end

function Base.convert(
    ::Type{ScalarAffineFunction{T}},
    f::ScalarQuadraticFunction{T},
) where {T}
    if !Base.isempty(f.quadratic_terms)
        throw(InexactError(:convert, ScalarAffineFunction{T}, f))
    end
    return ScalarAffineFunction{T}(f.affine_terms, f.constant)
end

# Conversion to ScalarQuadraticFunction
function Base.convert(::Type{ScalarQuadraticFunction{T}}, α::T) where {T}
    return ScalarQuadraticFunction{T}(
        ScalarAffineTerm{T}[],
        ScalarQuadraticTerm{T}[],
        α,
    )
end

function Base.convert(
    ::Type{ScalarQuadraticFunction{T}},
    f::SingleVariable,
) where {T}
    return convert(
        ScalarQuadraticFunction{T},
        convert(ScalarAffineFunction{T}, f),
    )
end

function Base.convert(
    ::Type{ScalarQuadraticFunction{T}},
    f::ScalarAffineFunction{T},
) where {T}
    return ScalarQuadraticFunction{T}(
        f.terms,
        ScalarQuadraticTerm{T}[],
        f.constant,
    )
end

function Base.convert(::Type{VectorOfVariables}, g::SingleVariable)
    return VectorOfVariables([g.variable])
end

function Base.convert(
    ::Type{VectorAffineFunction{T}},
    g::SingleVariable,
) where {T}
    return VectorAffineFunction{T}(
        [VectorAffineTerm(1, ScalarAffineTerm(one(T), g.variable))],
        [zero(T)],
    )
end

function Base.convert(
    ::Type{VectorQuadraticFunction{T}},
    g::SingleVariable,
) where {T}
    return VectorQuadraticFunction{T}(
        [VectorAffineTerm(1, ScalarAffineTerm(one(T), g.variable))],
        VectorQuadraticTerm{T}[],
        [zero(T)],
    )
end

function Base.convert(
    ::Type{VectorAffineFunction{T}},
    g::ScalarAffineFunction,
) where {T}
    return VectorAffineFunction{T}(
        VectorAffineTerm{T}[VectorAffineTerm(1, term) for term in g.terms],
        [g.constant],
    )
end

function Base.convert(
    ::Type{VectorQuadraticFunction{T}},
    g::ScalarAffineFunction,
) where {T}
    return VectorQuadraticFunction{T}(
        VectorAffineTerm{T}[VectorAffineTerm(1, term) for term in g.terms],
        VectorQuadraticTerm{T}[],
        [g.constant],
    )
end

function Base.convert(
    ::Type{VectorQuadraticFunction{T}},
    g::ScalarQuadraticFunction,
) where {T}
    return VectorQuadraticFunction{T}(
        VectorAffineTerm{T}[
            VectorAffineTerm(1, term) for term in g.affine_terms
        ],
        VectorQuadraticTerm{T}[
            VectorQuadraticTerm(1, term) for term in g.quadratic_terms
        ],
        [g.constant],
    )
end