text:
using Revise
using PyPlot
using Infiltrator
using LinearAlgebra
using Bem2d

"""
    gravityparticularfunctions()

From Pape and Bannerjee 1987
"""
function gravityparticularfunctions(x, y, g, rho, lambda, mu)
    U = zeros(length(x), 2)
    S = zeros(length(x), 3)
    U[:, 1] = @. -lambda * rho * g / (4 * mu * (lambda + mu)) * x * y  # Pape and Banerjee (1987) equation (6a)
    U[:, 2] = @. (rho * g) / (8 * mu * (lambda + mu)) * (lambda * x^2 + (lambda + 2 * mu) * y^2) # Pape and Banerjee (1987) equation (6b)
    S[:, 1] .= 0 # Pape and Banerjee (1987) equation (6c)
    S[:, 2] = @. rho * g * y  # Pape and Banerjee (1987) equation (6d)
    S[:, 3] .= 0 # Pape and Banerjee (1987) equation (6e)
    return U, S
end

"""
    gravitysquareparticular()

Experiments with gravity body force.
"""
function gravitysquareparticular()
    # TODO: Move particular solution to Bem2d.jl
    # TODO: It's strange that the top of the model has to be at zero.  Can we generalize this?
    # TODO: Rule of thumb for choosing precondtioner value (alpha)?
    # TODO: Make a version of gravityparticularfunctions() for BC generation

    close("all")
    alpha = 7e-8 # scalar preconditioner for traction terms
    fontsize = 20
    mu = 3e10
    lambda = 3e10
    nu = 0.25
    rho = 2700
    g = 9.81
    nels = 20
    npts = 100
    L = 1e4
    offset = 10
    x, y = obsgrid(-L+offset, -2*L+offset, L-offset, 0-offset, npts)

    # Define BEM geometry
    els = Elements(Int(1e5))
    x1, y1, x2, y2 = discretizedline(-L, -2*L, L, -2*L, nels) # Bottom
    addelsez!(els, x1, y1, x2, y2, "B")
    x1, y1, x2, y2 = discretizedline(L, -2*L, L, 0, nels) # Right hand side
    addelsez!(els, x1, y1, x2, y2, "R")
    x1, y1, x2, y2 = discretizedline(L, 0, -L, 0, nels) # Top
    addelsez!(els, x1, y1, x2, y2, "T")
    x1, y1, x2, y2 = discretizedline(-L, 0, -L, -2*L, nels) # Left hand side
    addelsez!(els, x1, y1, x2, y2, "L")

    # Common indexing
    idx = getidxdict(els) # Should this return "all" - YES TODO
    bcidxU = idx["B"] # Boundaries with *displacement* BCs
    bcidxT = [idx["R"] ; idx["T"]; idx["L"]] # Boundaries with *traction* BCs
    bcidxall = collect(1:1:els.endidx) # All boundaries

    # Gravity square problem with quadratic elements
    T_pU_qall, _ = PUTQ(slip2dispstress, els, bcidxU, bcidxall, mu, nu)
    _, H_pT_qall = PUTQ(slip2dispstress, els, bcidxT, bcidxall, mu, nu)

    TH = [T_pU_qall ; alpha .* H_pT_qall] # Assemble combined linear operator

    # Particular solution and effective boundary conditions
    xnodes = transpose(els.xnodes[idx["B"], :])[:]
    ynodes = transpose(els.ynodes[idx["B"], :])[:]
    UB, _ = gravityparticularfunctions(xnodes, ynodes, g, rho, lambda, mu)
    bcs = zeros(6 * els.endidx)
    bcs[1:2:6*nels] = UB[:, 1] # Bottom boundary (x-component)
    bcs[2:2:6*nels] = UB[:, 2] # Bottom boundary (y-component)
    bcs *= -1 # This is neccesary for the right answer and is consistent with derivation

    # BEM solve to get particular solution
    Ueffparticular = inv(TH) * bcs

    # Evaluate and plot interior solution
    Uinteriorcomplementary, Sinteriorcomplementary = quaddispstress(slip2dispstress, x, y, els, bcidxall, quadstack(Ueffparticular[1:2:end]), quadstack(Ueffparticular[2:2:end]), mu, nu)
    Uinteriorparticular, Sinteriorparticular = gravityparticularfunctions(x, y, g, rho, lambda, mu)
    U = @. Uinteriorcomplementary + Uinteriorparticular
    S = @. Sinteriorcomplementary + Sinteriorparticular
    plotfields(els, reshape(x, npts, npts), reshape(y, npts, npts), Uinteriorcomplementary, Sinteriorcomplementary, "Complementary solution")
    plotfields(els, reshape(x, npts, npts), reshape(y, npts, npts), Uinteriorparticular, Sinteriorparticular, "Particular solution")
    plotfields(els, reshape(x, npts, npts), reshape(y, npts, npts), U, S, "Complementary + Particular solutions")

    # Interior point solution for nice vector plot
    npts = 10
    offset = 500
    x, y = obsgrid(-L+offset, -2*L+offset, L-offset, 0-offset, npts)
    Uinteriorcomplementary, Sinteriorcomplementary = quaddispstress(slip2dispstress, x, y, els, bcidxall, quadstack(Ueffparticular[1:2:end]), quadstack(Ueffparticular[2:2:end]), mu, nu)
    Uinteriorparticular, Sinteriorparticular = gravityparticularfunctions(x, y, g, rho, lambda, mu)
    U = @. Uinteriorcomplementary + Uinteriorparticular
    S = @. Sinteriorcomplementary + Sinteriorparticular

    # Plot of interior points
    figure()
    fill([-10000, 10000, 10000, -10000], [-20000, -20000, 0, 0] .+ 20000, color="lightgray", zorder=1)
    quiver(x[:], y[:].+20000, U[:, 1], U[:, 2], units="width", color="cyan", edgecolor="k", linewidths=0.25, zorder=10)
    xlabel("x (m)")
    ylabel("y (m)")
    gca().set_aspect("equal")
end
gravitysquareparticular()


text:
module ConvexBodyProximityQueries

export closest_points
export minimum_distance
export tolerance_verification
export collision_detection

export ConvexPolygon
export point, line, rect, square, randpoly
export @point, @line, @rect, @square

using LinearAlgebra
using StaticArrays
using Random: shuffle!

include("support.jl")
include("obstacles.jl")
include("simplex.jl")
include("proximity.jl")
include("plot.jl")

using Requires

function __init__()
    @require LazySets = "b4f0291d-fe17-52bc-9479-3d1a343d9043" include("lazysets.jl")
end

end # module

ext: jl


text:
function findpkgpath()
    for depot in DEPOT_PATH
        path = joinpath(depot, "dev", "Pkg", "src", "Pkg.jl")
        isfile(path) && return path
    end
    error("Pkg.jl is not checked out at `dev/Pkg` in a depot (e.g., `~/.julia/dev/Pkg`).")
end

code = """
const pkgpath = $(repr(findpkgpath()))
"""

write(joinpath(@__DIR__, "deps.jl"), code)


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.

"""
     ADAM

The solver described in Diederik Kingma, Jimmy Ba: *Adam: A Method for
Stochastic Optimization*. arXiv:1412.6980 [cs.LG].

    ADAM(; kwargs...)

### Arguments
* `η`: default `0.001`, learning rate.
* `β1`: default `0.9`.
* `β2`: default `0.999`.
* `ϵ`: default `1e-8`.
* `clip`: default `0`, gradient clipping.
  If positive, will clip the gradient into the range `[-clip, clip]`.
* `scale`: default `0`, gradient rescaling.
  If != 0, multiply the gradient with `scale` before updating.
  Often choose to be `1.0 / batch_size`.
  If leave it default, high-level API like `fit!` will set it to
  `1.0 / batch_size`, since `fit!` knows the `batch_size`.
* `λ`: default `0.00001`, weight decay is equivalent
  to adding a global l2 regularizer for all the parameters.
* `η_sched::AbstractLearningRateScheduler`: default `LearningRate.Fixed(η)`, a
  dynamic learning rate scheduler. If set, will overwrite the `η` parameter.
"""
ADAM

@defstruct ADAM <: AbstractOptimizer (
  (η      :: Real = 0.001, η > 0),
  (β1     :: Real = 0.9,   0 <= β1 < 1),
  (β2     :: Real = 0.999, 0 <= β2 < 1),
  (ϵ      :: Real = 1e-8,  ϵ > 0),
  (clip   :: Real = 0,     clip >= 0),
   scale  :: Real = 0,
  (λ      :: Real = 1e-5,  λ >= 0),
  η_sched :: Any  = initlrsched(η)
)

mutable struct ADAMState
  η   :: Float64  # current learning rate
  mₜ  :: NDArray
  vₜ  :: NDArray
  β1ᵗ :: Float64
  β2ᵗ :: Float64
end

create_state(adam::ADAM, ::Int, W::NDArray) =
  ADAMState(get(adam.η_sched),
            zeros(size(W), context(W)),
            zeros(size(W), context(W)),
            adam.β1, adam.β2)

function update!(adam::ADAM, ::Int, W::NDArray, ∇:: NDArray, s::ADAMState)
  η = s.η
  β1 = adam.β1
  β2 = adam.β2
  ϵ = adam.ϵ

  normgrad!(adam, W, ∇)

  s.mₜ = β1 * s.mₜ + (1 - β1) .* ∇
  s.vₜ = β2 * s.vₜ + (1 - β2) .* ∇.^2

  aₜ= sqrt(1.0 - s.β2ᵗ)/(1.0 - s.β1ᵗ)

  # update βᵗ to βᵗ⁺¹
  s.β1ᵗ *= β1
  s.β2ᵗ *= β2

  @inplace W .+= -η * aₜ * s.mₜ ./ (sqrt(s.vₜ) .+ ϵ)
end

max_stars_repo_name: Vikas-kum/incubator-mxnet

max_issues_repo_name: Vikas-kum/incubator-mxnet

max_forks_repo_name: Vikas-kum/incubator-mxnet


text:
#=
test_linked_lists:
- Julia version: 1.5.0
- Author: shisa
- Date: 2020-09-03
=#

module TestLinkedLists
    # packages
    using Test

    # external modules
    include("create_linked_list.jl")
    include("circular_linked_list.jl")
    include("delete_node.jl")

    # methods
    function test()
        @testset "LinkedLists" begin
            @testset "LinkedListCreation" begin
                @test_nowarn LinkedListCreation.main()
            end
            @testset "NodeDelete" begin
                @test_nowarn NodeDelete.main()
            end
            @testset "CircularLinkedList" begin
                @test_nowarn CircularLinkedList.main()
            end
        end
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    using .TestLinkedLists
    TestLinkedLists.test()
end

max_stars_repo_name: ShisatoYano/JuliaPractice

max_issues_repo_name: ShisatoYano/JuliaPractice

max_forks_repo_name: ShisatoYano/JuliaPractice


text:
"""
	$(SIGNATURES)

Abstract Uniform type.
"""
abstract type AbstractUniform{T1} <: AbstractDistributionLH{T1} end

"""
	$(SIGNATURES)

Abstract type for switches governing how Uniform distributions are parameterized.
"""
abstract type AbstractUniformSwitches{T1} <: AbstractDistributionSwitches{T1} end

## -------------------  Generic

draw(u :: AbstractUniform{T1}, nDims, rng :: AbstractRNG) where
    T1 <: AbstractFloat =
    min_value(u) .+ value_range(u) .* rand(rng, T1, nDims...);

draw(u :: AbstractUniform{T1}, rng :: AbstractRNG)  where T1 <: AbstractFloat =
    min_value(u) + value_range(u) * rand(rng, T1);

quantiles(u :: AbstractUniform{T1}, pctV)  where T1 <: AbstractFloat =
    min_value(u) .+ value_range(u) .* T1.(pctV);

Base.eltype(::AbstractUniform{T1}) where T1 <: AbstractFloat = T1;

isbounded(u :: AbstractUniform{T1}) where T1 = true;

Base.show(io :: IO, um :: AbstractUniform{T1}) where T1 =
    print(io, "Uniform{$T1} over range $(min_value(um)) to $(max_value(um))");

# Validate `Uniform` draws
function check_draws(u :: AbstractUniform{T1}, drawM) where T1
    isValid = true;
    if any(x -> x < min_value(u), drawM)
        @warn "Below $(min_value(u)):  $(minimum(drawM))"
        isValid = false;
    end

    if any(x -> x > max_value(u), drawM)
        @warn "Above $(max_value(u)):  $(maximum(drawM))"
        isValid = false;
    end
    return isValid
end

## ----------------  Uniform with minimum and range

# The defaults really don't make all that much sense, but we want a keyword constructor
Base.@kwdef mutable struct UniformSwitches{T1 <: AbstractFloat} <: AbstractUniformSwitches{T1}
    xMin :: T1 = 1.0
    xMinLb :: T1 = 0.1
    xMinUb :: T1 = 100.0
    xMinDescription :: String = "xMin"
    xMinLatex :: String = "xMin"
    xRange :: T1 = 3.0
    xRangeLb :: T1 = 0.1
    xRangeUb :: T1 = 100.0
    xRangeDescription :: String = "xRange"
    xRangeLatex :: String = "xRange"
    calXMin :: Bool = true
    calXRange :: Bool = true
end

"""
	$(SIGNATURES)

Unform distribution on interval `xMin + [0, xRange]`.
"""
mutable struct Uniform{T1} <: AbstractUniform{T1}
    objId :: ObjectId
    xMin :: T1
    # Difference between h0Max and h0Min
    xRange :: T1
    pvec :: ParamVector
end

min_value(u :: Uniform{T1}) where T1 = u.xMin;
max_value(u :: Uniform{T1}) where T1 = u.xMin + u.xRange;
value_range(u :: Uniform{T1}) where T1 = u.xRange;

## ----------  Construction

init_distribution(objId :: ObjectId, switches :: AbstractUniformSwitches{T1}) where T1 =
    init_uniform(objId, switches);

"""
	$(SIGNATURES)

Construct a Uniform distribution from its switches.
"""
function init_uniform(objId :: ObjectId, switches :: UniformSwitches{T1}) where T1
    pMin = init_xmin(switches);
    pRange = init_xrange(switches);
    pvec = ParamVector(objId, [pMin, pRange]);
    xMin = ModelParams.value(pMin);
    return Uniform(objId, xMin, ModelParams.value(pRange), pvec)
end

init_xmin(s :: UniformSwitches{T1}) where T1 =
    Param(:xMin, s.xMinDescription, s.xMinLatex, s.xMin, s.xMin,
        s.xMinLb, s.xMinUb, s.calXMin);

init_xrange(s :: AbstractUniformSwitches{T1}) where T1 =
    Param(:xRange, s.xRangeDescription, s.xRangeLatex, s.xRange, s.xRange,
        s.xRangeLb, s.xRangeUb, s.calXRange);

make_test_uniform_switches() = UniformSwitches();
make_test_uniform() = init_uniform(ObjectId(:none),  make_test_uniform_switches());

"""
	$(SIGNATURES)

Uniform distribution with fixed bounds (not calibrated).
"""
function UniformFixedBounds(objId :: ObjectId, lb :: T1, ub :: T1) where
    T1 <: AbstractFloat

    dx = ub - lb;
    @assert dx > 0
    s = UniformSwitches(
        xMin = lb, xMinLb = lb - 0.1, xMinUb = lb + 0.1, calXMin = false,
        xRange = dx, xRangeLb = dx - 0.1, xRangeUb = dx + 0.1, calXRange = false);
    return init_uniform(objId, s);
end

make_test_uniform_fixed() = UniformFixedBounds(ObjectId(:test), -1.0, 2.0);

## ----------------  Uniform with mean and range

"""
	$(SIGNATURES)

Switches for Uniform distribution that is centered around a mean.
"""
Base.@kwdef mutable struct UniformCenteredSwitches{T1 <: AbstractFloat} <: AbstractUniformSwitches{T1}
    xMean :: T1 = 1.0
    xMeanLb :: T1 = 0.1
    xMeanUb :: T1 = 100.0
    xMeanDescription :: String = "xMean"
    xMeanLatex :: String = "xMean"
    xRange :: T1 = 3.0
    xRangeLb :: T1 = 0.1
    xRangeUb :: T1 = 100.0
    xRangeDescription :: String = "xRange"
    xRangeLatex :: String = "xRange"
    calXMean :: Bool = true
    calXRange :: Bool = true
end

"""
	$(SIGNATURES)

Unform distribution on interval `xMean +/- 0.5 * xRange`.
"""
mutable struct UniformCentered{T1} <: AbstractUniform{T1}
    objId :: ObjectId
    xMean :: T1
    # Difference between h0Max and h0Min
    xRange :: T1
    pvec :: ParamVector
end

min_value(u :: UniformCentered{T1}) where T1 = u.xMean - 0.5 * u.xRange;
max_value(u :: UniformCentered{T1}) where T1 = (u.xMean + 0.5 * u.xRange);
value_range(u :: UniformCentered{T1}) where T1 = u.xRange;

## ----------  Construction

function init_uniform(objId :: ObjectId, switches :: UniformCenteredSwitches)
    pMean = init_xmean(switches);
    pRange = init_xrange(switches);
    pvec = ParamVector(objId, [pMean, pRange]);
    xMean = ModelParams.value(pMean);
    return UniformCentered(objId, xMean, ModelParams.value(pRange), pvec)
end

init_xmean(s :: UniformCenteredSwitches) =
    Param(:xMean, s.xMeanDescription, s.xMeanLatex, s.xMean, s.xMean,
        s.xMeanLb, s.xMeanUb, s.calXMean);

make_test_uniform_centered_switches() = UniformCenteredSwitches();
make_test_uniform_centered() = init_uniform(ObjectId(:none),  make_test_uniform_centered_switches());

# Count no of observations in quantile groups
# Include 0 as lower bound of first percentile group.
# function count_quantiles(u :: Uniform, drawV, pctV)
#     @assert pctV[1] == 0.0  "Include 0 lower bound"
#     h = fit(Histogram, vec(drawV), quantiles(u, pctV));
#     return h.weights
# end

# function scale_normal_draws(hmS :: AbstractUniform{F},  inM :: T)  where {F, T}
#     cdfM = cdf.(Normal(), inM);
#     drawM = quantile.(Distributions.Uniform(min_value(hmS), max_value(hmS)),  cdfM);
#     return drawM
# end

# -------------


text:
function application()
  WebRTCApplication(randstring())
end


text:
function heatmap(mat::Array{Float64,2})
  Nrow, Ncol = size(mat)
  xvals = Array(Float64,length(mat))
  yvals = Array(Float64,length(mat))
  colorvals = Array(Float64,length(mat))
  index = 1
  for x in 1:Ncol
    for y in 1:Nrow
      xvals[index] = x
      yvals[index] = y
      colorvals[index] = mat[y,x]
      index+=1
    end
  end
  plot(x=xvals,y=yvals,color=colorvals, Geom.rectbin,
  	Scale.x_discrete, Scale.y_discrete,
	Coord.cartesian(yflip=true,fixed=true),
	Guide.xlabel(""),Guide.ylabel(""),
	Scale.ContinuousColorScale(Scale.lab_gradient(colorant"white",colorant"red"),minvalue=0.))
end

function writeMatrix(filename::ASCIIString,matrix::Matrix)
  stream = open(filename,"w+")
  write(stream,mat)
  close(stream)
end

function readMatrix{T<:Real}(filename::ASCIIString,format::Type{T},dims::Tuple{Vararg{Int64}})
  stream = open(filename,"r")
  matrix = read(stream,format,dims)
  close(stream)
  return matrix
end


text:
abstract type BibliographyBlocks <: Expanders.ExpanderPipeline end

Selectors.order(::Type{BibliographyBlocks}) = 6.0
Selectors.matcher(::Type{BibliographyBlocks}, node, page, doc) = Expanders.iscode(node, "@bibliography")
function Selectors.runner(::Type{BibliographyBlocks}, x, page, doc)
    node = Documents.buildnode(BibliographyNode, x, doc, page)
    settings = Documents.getplugin(doc, CitationPlugin)
    push!(settings.bibnodes, node)
    node.style === nothing && (node.style = settings.style)
    page.mapping[x] = node
end


text:
# an artificial example illustrating the necessity of an extra projection in the multi-output case
using Logging

using StructuralIdentifiability

logger = Logging.SimpleLogger(stdout, Logging.Info)
global_logger(logger)

ode = @ODEmodel(
    x1'(t) = (1 + x1(t)^2) // 2,
    x2'(t) = (1 - x1(t)^2) // (1 + x1(t)^2),
    y1(t) = 2 * x1(t) // (b * (1 + x1(t)^2)),
    y2(t) = x2(t)
)

@time println(assess_global_identifiability(ode))


text:
#using Plots
#println("Creating animation of the obstruction")
#anim = @animate for i=1:2*p.N2-1
#               plot(real(Uint[:,abs(p.N2-i)+1,2,2]),imag(Uint[:,abs(p.N2-i)+1,2,2]),xlims=[-1,1],ylims=[-1,1])
#end

#gif(anim,"obs.gif",fps=25)
#close("all")
#Ucol1 = zeros(size(Uint)[1],size(Uint)[2],3)
#
#for i = 1:p.N1, j=1:p.N2
#	if Uint[i,j,2,1] != 0.0
#		phase = imag(log(im*Uint[i,j,2,1]/(abs(Uint[i,j,2,1]))))
#	else
#		phase = 0.0
#	end
#	array = Uint[i,j,:,1].*exp(-im*phase)
#
#	Ucol1[i,j,:] = [real(array[1]),imag(array[1]),imag(array[2])]
#end

max_stars_repo_name: antoine-levitt/wannier

max_issues_repo_name: antoine-levitt/wannier

max_forks_repo_name: antoine-levitt/wannier


text:
#Benchmarking of RBC and FVGQ variants
using DifferenceEquations, BenchmarkTools, LinearAlgebra
using DelimitedFiles, Distributions, Zygote
# General likelihood calculation
function make_problem_2(A_0, A_1, A_2, B, C_0, C_1, C_2, u0, noise, observables, D; kwargs...)
    prob = QuadraticStateSpaceProblem(A_0, A_1, A_2, B, u0, (0, size(observables, 2)); C_0, C_1,
                                      C_2, observables_noise = MvNormal(Diagonal(abs2.(D))), noise,
                                      observables, kwargs...)
    return prob.A_1[1, 1] + prob.B[1, 1]
end

function joint_likelihood_2(A_0, A_1, A_2, B, C_0, C_1, C_2, u0, noise, observables, D; kwargs...)
    prob = QuadraticStateSpaceProblem(A_0, A_1, A_2, B, u0, (0, size(observables, 2)); C_0, C_1,
                                      C_2, observables_noise = MvNormal(Diagonal(abs2.(D))), noise,
                                      observables, kwargs...)
    return solve(prob).logpdf
end

function simulate_model_no_noise_2(A_0, A_1, A_2, B, C_0, C_1, C_2, u0, observables, D; kwargs...)
    prob = QuadraticStateSpaceProblem(A_0, A_1, A_2, B, u0, (0, size(observables, 2)); C_0, C_1,
                                      C_2, observables_noise = MvNormal(Diagonal(abs2.(D))),
                                      observables, kwargs...)
    sol = solve(prob)
    return sol.retcode
end

function simulate_model_no_observations_2(A_0, A_1, A_2, B, C_0, C_1, C_2, u0, T; kwargs...)
    prob = QuadraticStateSpaceProblem(A_0, A_1, A_2, B, u0, (0, T); C_0, C_1, C_2, kwargs...)
    sol = solve(prob)
    return sol.retcode
end

const QUADRATIC = BenchmarkGroup()

# Matrices from RBC
const A_0_rbc = [-7.824904812740593e-5, 0.0]
const A_1_rbc = [0.9568351489231076 6.209371005755285; 3.0153731819288737e-18 0.20000000000000007]
const A_2_rbc = cat([-0.00019761505863889124 0.03375055315837927; 0.0 0.0],
                    [0.03375055315837913 3.128758481817603; 0.0 0.0]; dims = 3)
const B_2_rbc = reshape([0.0; -0.01], 2, 1) # make sure B is a matrix
const C_0_rbc = [7.824904812740593e-5, 0.0]
const C_1_rbc = [0.09579643002426148 0.6746869652592109; 1.0 0.0]
const C_2_rbc = cat([-0.00018554166974717046 0.0025652363153049716; 0.0 0.0],
                    [0.002565236315304951 0.3132705036896446; 0.0 0.0]; dims = 3)
const D_2_rbc = [0.1, 0.1]
const u0_2_rbc = zeros(2)

const observables_2_rbc = readdlm(joinpath(pkgdir(DifferenceEquations),
                                           "test/data/RBC_observables.csv"), ',')' |> collect
const noise_2_rbc = readdlm(joinpath(pkgdir(DifferenceEquations), "test/data/RBC_noise.csv"),
                            ',')' |> collect
const T_2_rbc = size(observables_2_rbc, 2)
# Matrices from FVGQ
# Load FVGQ data for checks
A_0_raw = readdlm(joinpath(pkgdir(DifferenceEquations), "test/data/FVGQ20_A_0.csv"), ',')
const A_0_FVGQ = vec(A_0_raw)
const A_1_FVGQ = readdlm(joinpath(pkgdir(DifferenceEquations), "test/data/FVGQ20_A_1.csv"), ',')
A_2_raw = readdlm(joinpath(pkgdir(DifferenceEquations), "test/data/FVGQ20_A_2.csv"), ',')
const A_2_FVGQ = reshape(A_2_raw, length(A_0_FVGQ), length(A_0_FVGQ), length(A_0_FVGQ))
const B_2_FVGQ = readdlm(joinpath(pkgdir(DifferenceEquations), "test/data/FVGQ20_B.csv"), ',')
C_0_raw = readdlm(joinpath(pkgdir(DifferenceEquations), "test/data/FVGQ20_C_0.csv"), ',')
const C_0_FVGQ = vec(C_0_raw)
const C_1_FVGQ = readdlm(joinpath(pkgdir(DifferenceEquations), "test/data/FVGQ20_C_1.csv"), ',')
C_2_raw = readdlm(joinpath(pkgdir(DifferenceEquations), "test/data/FVGQ20_C_2.csv"), ',')
const C_2_FVGQ = reshape(C_2_raw, length(C_0_FVGQ), length(A_0_FVGQ), length(A_0_FVGQ))
# D_raw = readdlm(joinpath(pkgdir(DifferenceEquations), "FVGQ_D.csv"); header = false)))
D_2_FVGQ = ones(6) * 1e-3
u0_2_FVGQ = zeros(size(A_1_FVGQ, 1))
const observables_2_FVGQ = readdlm(joinpath(pkgdir(DifferenceEquations),
                                            "test/data/FVGQ20_observables.csv"), ',')' |> collect

const noise_2_FVGQ = readdlm(joinpath(pkgdir(DifferenceEquations), "test/data/FVGQ20_noise.csv"),
                             ',')' |> collect
const T_2_FVGQ = size(observables_2_FVGQ, 2)

# RBC sized specific tests
# Verifying code prior to benchmark
# executing gradients once to avoid compilation time in benchmarking
make_problem_2(A_0_rbc, A_1_rbc, A_2_rbc, B_2_rbc, C_0_rbc, C_1_rbc, C_2_rbc, u0_2_rbc, noise_2_rbc,
               observables_2_rbc, D_2_rbc)
gradient((args...) -> make_problem_2(args..., observables_2_rbc, D_2_rbc), A_0_rbc, A_1_rbc,
         A_2_rbc, B_2_rbc, C_0_rbc, C_1_rbc, C_2_rbc, u0_2_rbc, noise_2_rbc)
make_problem_2(A_0_FVGQ, A_1_FVGQ, A_2_FVGQ, B_2_FVGQ, C_0_FVGQ, C_1_FVGQ, C_2_FVGQ, u0_2_FVGQ,
               noise_2_FVGQ, observables_2_FVGQ, D_2_FVGQ)
gradient((args...) -> make_problem_2(args..., observables_2_FVGQ, D_2_FVGQ), A_0_FVGQ, A_1_FVGQ,
         A_2_FVGQ, B_2_FVGQ, C_0_FVGQ, C_1_FVGQ, C_2_FVGQ, u0_2_FVGQ, noise_2_FVGQ)

joint_likelihood_2(A_0_rbc, A_1_rbc, A_2_rbc, B_2_rbc, C_0_rbc, C_1_rbc, C_2_rbc, u0_2_rbc,
                   noise_2_rbc, observables_2_rbc, D_2_rbc)
gradient((args...) -> joint_likelihood_2(args..., observables_2_rbc, D_2_rbc), A_0_rbc, A_1_rbc,
         A_2_rbc, B_2_rbc, C_0_rbc, C_1_rbc, C_2_rbc, u0_2_rbc, noise_2_rbc)
joint_likelihood_2(A_0_FVGQ, A_1_FVGQ, A_2_FVGQ, B_2_FVGQ, C_0_FVGQ, C_1_FVGQ, C_2_FVGQ, u0_2_FVGQ,
                   noise_2_FVGQ, observables_2_FVGQ, D_2_FVGQ)
gradient((args...) -> joint_likelihood_2(args..., observables_2_FVGQ, D_2_FVGQ), A_0_FVGQ, A_1_FVGQ,
         A_2_FVGQ, B_2_FVGQ, C_0_FVGQ, C_1_FVGQ, C_2_FVGQ, u0_2_FVGQ, noise_2_FVGQ)

const QUADRATIC["rbc"] = BenchmarkGroup()
const QUADRATIC["rbc"]["make_problem_2"] = @benchmarkable make_problem_2($A_0_rbc, $A_1_rbc,
                                                                         $A_2_rbc, $B_2_rbc,
                                                                         $C_0_rbc, $C_1_rbc,
                                                                         $C_2_rbc, $u0_2_rbc,
                                                                         $noise_2_rbc,
                                                                         $observables_2_rbc,
                                                                         $D_2_rbc)
const QUADRATIC["rbc"]["make_problem_2_gradient"] = @benchmarkable gradient((args...) -> make_problem_2(args...,
                                                                                $observables_2_rbc,
                                                                                $D_2_rbc),
                                                                            $A_0_rbc, $A_1_rbc,
                                                                            $A_2_rbc, $B_2_rbc,
                                                                            $C_0_rbc, $C_1_rbc,
                                                                            $C_2_rbc, $u0_2_rbc,
                                                                            $noise_2_rbc)

const QUADRATIC["rbc"]["simulate_model_no_noise_2"] = @benchmarkable simulate_model_no_noise_2($A_0_rbc,
                                                                                $A_1_rbc,
                                                                                $A_2_rbc,
                                                                                $B_2_rbc,
                                                                                $C_0_rbc,
                                                                                $C_1_rbc,
                                                                                $C_2_rbc,
                                                                                $u0_2_rbc,
                                                                                $observables_2_rbc,
                                                                                $D_2_rbc)
const QUADRATIC["rbc"]["simulate_model_no_observations_2"] = @benchmarkable simulate_model_no_observations_2($A_0_rbc,
                                                                                $A_1_rbc,
                                                                                $A_2_rbc,
                                                                                $B_2_rbc,
                                                                                $C_0_rbc,
                                                                                $C_1_rbc,
                                                                                $C_2_rbc,
                                                                                $u0_2_rbc,
                                                                                $T_2_rbc)
const QUADRATIC["rbc"]["joint_2"] = @benchmarkable joint_likelihood_2($A_0_rbc, $A_1_rbc, $A_2_rbc,
                                                                      $B_2_rbc, $C_0_rbc, $C_1_rbc,
                                                                      $C_2_rbc, $u0_2_rbc,
                                                                      $noise_2_rbc,
                                                                      $observables_2_rbc, $D_2_rbc)
const QUADRATIC["rbc"]["joint_2_gradient"] = @benchmarkable gradient((args...) -> joint_likelihood_2(args...,
                                                                                $observables_2_rbc,
                                                                                $D_2_rbc),
                                                                     $A_0_rbc, $A_1_rbc, $A_2_rbc,
                                                                     $B_2_rbc, $C_0_rbc, $C_1_rbc,
                                                                     $C_2_rbc, $u0_2_rbc,
                                                                     $noise_2_rbc)

# FVGQ sized specific test
const QUADRATIC["FVGQ"] = BenchmarkGroup()
const QUADRATIC["FVGQ"]["make_problem_2"] = @benchmarkable make_problem_2($A_0_FVGQ, $A_1_FVGQ,
                                                                          $A_2_FVGQ, $B_2_FVGQ,
                                                                          $C_0_FVGQ, $C_1_FVGQ,
                                                                          $C_2_FVGQ, $u0_2_FVGQ,
                                                                          $noise_2_FVGQ,
                                                                          $observables_2_FVGQ,
                                                                          $D_2_FVGQ)
const QUADRATIC["FVGQ"]["make_problem_2_gradient"] = @benchmarkable gradient((args...) -> make_problem_2(args...,
                                                                                $observables_2_FVGQ,
                                                                                $D_2_FVGQ),
                                                                             $A_0_FVGQ, $A_1_FVGQ,
                                                                             $A_2_FVGQ, $B_2_FVGQ,
                                                                             $C_0_FVGQ, $C_1_FVGQ,
                                                                             $C_2_FVGQ, $u0_2_FVGQ,
                                                                             $noise_2_FVGQ)

const QUADRATIC["FVGQ"]["simulate_model_no_noise_2"] = @benchmarkable simulate_model_no_noise_2($A_0_FVGQ,
                                                                                $A_1_FVGQ,
                                                                                $A_2_FVGQ,
                                                                                $B_2_FVGQ,
                                                                                $C_0_FVGQ,
                                                                                $C_1_FVGQ,
                                                                                $C_2_FVGQ,
                                                                                $u0_2_FVGQ,
                                                                                $observables_2_FVGQ,
                                                                                $D_2_FVGQ)
const QUADRATIC["FVGQ"]["simulate_model_no_observations_2"] = @benchmarkable simulate_model_no_observations_2($A_0_FVGQ,
                                                                                $A_1_FVGQ,
                                                                                $A_2_FVGQ,
                                                                                $B_2_FVGQ,
                                                                                $C_0_FVGQ,
                                                                                $C_1_FVGQ,
                                                                                $C_2_FVGQ,
                                                                                $u0_2_FVGQ,
                                                                                $T_2_FVGQ)

const QUADRATIC["FVGQ"]["joint_2"] = @benchmarkable joint_likelihood_2($A_0_FVGQ, $A_1_FVGQ,
                                                                       $A_2_FVGQ, $B_2_FVGQ,
                                                                       $C_0_FVGQ, $C_1_FVGQ,
                                                                       $C_2_FVGQ, $u0_2_FVGQ,
                                                                       $noise_2_FVGQ,
                                                                       $observables_2_FVGQ,
                                                                       $D_2_FVGQ)
const QUADRATIC["FVGQ"]["joint_2_gradient"] = @benchmarkable gradient((args...) -> joint_likelihood_2(args...,
                                                                                $observables_2_FVGQ,
                                                                                $D_2_FVGQ),
                                                                      $A_0_FVGQ, $A_1_FVGQ,
                                                                      $A_2_FVGQ, $B_2_FVGQ,
                                                                      $C_0_FVGQ, $C_1_FVGQ,
                                                                      $C_2_FVGQ, $u0_2_FVGQ,
                                                                      $noise_2_FVGQ)
# return for the test suite
QUADRATIC


text:
"""
    train(df, ycol, scaled_prefix, scaled_features,
        sample_size, input_size, batch_size, output_size, epoch_size,
        total_wd_idxs, test_wd_idxs,
        train_chnk, valid_idxs, test_idxs,
        statvals, filename, test_dates)

"""
function train_season_DNN(train_wd::Array{DataFrame, 1}, valid_wd::Array{DataFrame, 1}, test_wd::Array{DataFrame, 1},
    ycol::Symbol, scaled_ycol::Symbol, scaled_features::Array{Symbol}, scaling_method::Symbol,
    train_size::Integer, valid_size::Integer, test_size::Integer,
    sample_size::Integer, input_size::Integer, batch_size::Integer, output_size::Integer,
    epoch_size::Integer, _eltype::DataType,
    test_dates::Array{ZonedDateTime,1},
    statvals::AbstractNDSparse, season_table::AbstractNDSparse,
    output_prefix::String, filename::String) where I <: Integer

    @info "DNN training starts.."

    # extract from ndsparse
    total_μ = statvals[String(ycol) * "_res", "μ"].value
    total_σ = statvals[String(ycol) * "_res", "σ"].value
    total_max = float(statvals[String(ycol) * "_res", "maximum"].value)
    total_min = float(statvals[String(ycol) * "_res", "minimum"].value)

    # compute mean and std by each train/valid/test set
    statval = ndsparse((
        dataset = ["total", "total", "total", "total"],
        type = ["μ", "σ", "maximum", "minimum"]),
        (value = [total_μ, total_σ, total_max, total_min],))

    # construct compile function symbol
    compile = eval(Symbol(:compile, "_", ycol, "_season_DNN"))
    model, loss, accuracy, opt = compile(input_size, batch_size, output_size, statval)

    # modify scaled_features to train residuals

    # |> gpu doesn't work to *_set directly
    # construct minibatch for train_set
    # https://github.com/FluxML/Flux.jl/issues/704
    @info "    Construct Training Set batch..."
    p = Progress(div(length(train_wd), batch_size) + 1, dt=1.0, barglyphs=BarGlyphs("[=> ]"), barlen=40, color=:yellow)

    train_set = [(ProgressMeter.next!(p);
        make_batch_DNN(collect(dfs), scaled_ycol, scaled_features,
            sample_size, output_size, batch_size, 0.5, _eltype))
        for dfs in Base.Iterators.partition(train_wd, batch_size)]

    # don't construct minibatch for valid & test sets
    @info "    Construct Valid Set..."
    p = Progress(length(valid_wd), dt=1.0, barglyphs=BarGlyphs("[=> ]"), barlen=40, color=:yellow)
    valid_set = [(ProgressMeter.next!(p);
        make_batch_DNN(collect(dfs), scaled_ycol, scaled_features,
            sample_size, output_size, batch_size, 0.5, _eltype))
            for dfs in Base.Iterators.partition(valid_wd, batch_size)]

    @info "    Construct Test Set..."
    p = Progress(length(test_wd), dt=1.0, barglyphs=BarGlyphs("[=> ]"), barlen=40, color=:yellow)
    test_set = [(ProgressMeter.next!(p);
        make_pair_date_DNN(df, scaled_ycol, scaled_features, sample_size, output_size, _eltype)) for df in test_wd]

    train_set = gpu.(train_set)
    valid_set = gpu.(valid_set)
    test_set = gpu.(test_set)

    # *_set : normalized values
    df_evals = train_season_DNN!(model, train_set, valid_set,
        loss, accuracy, opt, epoch_size, statval, scaled_features, filename)

    # TODO : (current) validation with zscore, (future) validation with original value?
    @info "    Test ACC  : ", accuracy(test_set)
    @info "    Valid ACC : ", accuracy(valid_set)
    flush(stdout); flush(stderr)

    eval_metrics = [:RMSE, :MAE, :MSPE, :MAPE]
    # test set
    for metric in eval_metrics
        _eval = evaluation1(test_set, model, statval, metric)
        @info " $(string(ycol)) $(string(metric)) for test   : ", _eval
    end

    # valid set
    for metric in eval_metrics
        _eval = evaluation1(valid_set, model, statval, metric)
        @info " $(string(ycol)) $(string(metric)) for valid  : ", _eval
    end

    # create directory per each time
    for i = 1:output_size
        i_pad = lpad(i, 2, '0')
        Base.Filesystem.mkpath("/$(output_prefix)/$(i_pad)/")
    end

    # back to unnormalized for comparison
    if scaling_method == :zscore
        dnn_df = predict_DNN_model_zscore_season(test_set, model, ycol,
            total_μ, total_σ, season_table, _eltype, output_size, "/$(output_prefix)/")
    elseif scaling_method == :minmax
        # TODO : compose seasonality
        dnn_df = predict_DNN_model_minmax(test_set, model, ycol,
            total_min, total_max, 0.0, 10.0, _eltype, output_size, "/$(output_prefix)/")
    elseif scaling_method == :logzscore
        # TODO : compose seasonality
        log_μ = statvals[string(:log_, ycol), "μ"].value
        log_σ = statvals[string(:log_, ycol), "σ"].value
        dnn_df = predict_DNN_model_logzscore(test_set, model, ycol,
            log_μ, log_σ, _eltype, output_size, "/$(output_prefix)/")
    elseif scaling_method == :invzscore
        # TODO : compose seasonality
        inv_μ = statvals[string(:inv_, ycol), "μ"].value
        inv_σ = statvals[string(:inv_, ycol), "σ"].value
        dnn_df = predict_DNN_model_invzscore(test_set, model, ycol,
            inv_μ, inv_σ, _eltype, output_size, "/$(output_prefix)/")
    elseif scaling_method == :logminmax
        log_max = statvals[string(:log_, ycol), "maximum"].value
        log_min = statvals[string(:log_, ycol), "minimum"].value
        dnn_df = predict_DNN_model_logminmax(test_set, model, ycol,
            log_min, log_max, 0.0, 10.0, _eltype, output_size, "/$(output_prefix)/")
    end

    # 3 months plot
    # TODO : how to generalize date range? how to split based on test_dates?
    # 1/4 : because train size is 3 days, result should be start from 1/4
    # 12/29 : same reason 1/4, but this results ends with 12/31 00:00 ~ 12/31 23:00
    push!(eval_metrics, :learn_rate)
    plot_evaluation(df_evals, ycol, eval_metrics, "/$(output_prefix)/")

    model, dnn_df, statval
end

function train_season_DNN!(model::C,
    train_set::Array{T2, 1}, valid_set::Array{T2, 1},
    loss, accuracy, opt,
    epoch_size::Integer, statval::AbstractNDSparse, scaled_features::Array{Symbol},
    filename::String) where {C <: Flux.Chain, F <: AbstractFloat, T2 <: Tuple{AbstractArray{F, 2}, AbstractArray{F, 2}},
    T1 <: Tuple{AbstractArray{F, 1}, AbstractArray{F, 1}}}

    @info("    Beginning training loop...")
    flush(stdout); flush(stderr)

    _acc = Inf
    # for adjR2 -∞ is the worst
    best_acc = Inf
    last_improvement = 0

    df_eval = DataFrame(epoch = Int64[], learn_rate = Float64[], ACC = Float64[],
        RMSE = Float64[], MAE = Float64[], MSPE = Float64[], MAPE = Float64[])

    for epoch_idx in 1:epoch_size
        best_acc, last_improvement
        # train model with normalized data set
        Flux.train!(loss, Flux.params(model), train_set, opt)

        # record evaluation
        rmse, mae, mspe, mape =
            evaluations1(valid_set, model, statval,
            [:RMSE, :MAE, :MSPE, :MAPE])
        push!(df_eval, [epoch_idx opt.eta _acc rmse mae mspe mape])

        # Calculate accuracy:
        _acc = accuracy(valid_set) |> cpu
        _loss = loss(train_set[1][1], train_set[1][2]) |> cpu
        @info(@sprintf("epoch [%d]: loss[1]: %.8E Valid accuracy: %.8f Time: %s", epoch_idx, _loss, _acc, now()))
        flush(stdout); flush(stderr)

        # If our accuracy is good enough, quit out.
        if _acc < 0.001
            @info("    -> Early-exiting: We reached our target accuracy")
            break
        end

        # If this is the best accuracy we've seen so far, save the model out
        if _acc < best_acc
            @info "    -> New best accuracy! Saving model out to " * filename
            flush(stdout)

            cpu_model = model |> cpu
            weights = Flux.params(cpu_model)
            filepath = "/mnt/" * filename * ".bson"
            μ, σ = statval["total", "μ"].value, statval["total", "σ"].value
            total_max, total_min =
                float(statval["total", "maximum"].value), float(statval["total", "minimum"].value)
            # BSON can't save weights now (2019/12), disable temporarily
            BSON.@save filepath cpu_model epoch_idx _acc μ σ total_max total_min

            best_acc = _acc
            last_improvement = epoch_idx
        end

        # If we haven't seen improvement in 5 epochs, drop our learning rate:
        if epoch_idx - last_improvement >= 10 && opt.eta > 1e-6
            opt.eta /= 10.0
            @warn("    -> Haven't improved in a while, dropping learning rate to $(opt.eta)!")
            flush(stdout); flush(stderr)

            # After dropping learning rate, give it a few epochs to improve
            last_improvement = epoch_idx
        end

        if epoch_idx - last_improvement >= 20
            @warn("     -> We're calling this converged.")
            flush(stdout); flush(stderr)

            break
        end
    end

    df_eval
end

function compile_PM10_season_DNN(input_size::Integer, batch_size::Integer, output_size::Integer, statval::AbstractNDSparse)
    @info("    Compiling model...")

    unit_size = 16

    # https://machinelearningmastery.com/dropout-regularization-deep-learning-models-keras/
    # not sigmoid, elu used to reduce vanishing gradient problem
    # predict low concentration is not important than high concentration
    model = Chain(
        Dense(input_size, unit_size, leakyrelu),

        Dense(unit_size, unit_size, leakyrelu),

        Dense(unit_size, output_size)
    ) |> gpu

    GC.@preserve model

    @info "Unit size in PM10: ", unit_size
    @info "Model     in PM10: ", model

    # no regularization
    loss(x, y) = Flux.mse(model(x), y)

    # L1
    #loss(x, y) = Flux.mse(model(x), y) +
    #    lambda(x, y) *
    #    sum(x1 -> LinearAlgebra.norm(x1, 1), Flux.params(model))

    # L2
    # disable due to Flux#930 issue
    #lambda(x, y) = 10^(log10(Flux.mse(model(x), y)) - 1)
    #loss(x, y) = Flux.mse(model(x), y) + lambda(x, y) * sum(LinearAlgebra.norm, Flux.params(model))

    # Changhoon lee's method to capture outliers
    #loss(x, y) = sum(exp.(y .* 0.000001)) * Flux.mse(model(x), y)
    #loss(x, y) = maximum(10.0.^(y .* 0.0001)) * Flux.mse(model(x), y)

    accuracy(data) = evaluation1(data, model, statval, :RMSE)
    opt = Flux.ADAM()

    model, loss, accuracy, opt
end

function compile_PM25_season_DNN(input_size::Integer, batch_size::Integer, output_size::Integer, statval::AbstractNDSparse)
    @info("    Compiling model...")

    unit_size = 32

    # https://machinelearningmastery.com/dropout-regularization-deep-learning-models-keras/
    # not sigmoid, elu used to reduce vanishing gradient problem
    # predict low concentration is not important than high concentration
    model = Chain(
        Dense(input_size, unit_size, leakyrelu),

        Dense(unit_size, unit_size, leakyrelu),

        Dense(unit_size, output_size)
    ) |> gpu

    GC.@preserve model

    @info "Unit size in PM25: ", unit_size
    @info "Model     in PM25: ", model

    # no regularization
    loss(x, y) = Flux.mse(model(x), y)

    # L1 regularization
    #loss(x, y) = Flux.mse(model(x), y) +
    #    lambda(x, y) *
    #    sum(x1 -> LinearAlgebra.norm(x1, 1), Flux.params(model))

    # L2 regularization
    # disable regularization due to Flux#930 issue
    #lambda(x, y) = 10^(log10(Flux.mse(model(x), y)) - 1)
    #loss(x, y) = Flux.mse(model(x), y) + lambda(x, y) * sum(LinearAlgebra.norm, Flux.params(model))

    # Changhoon lee's method to capture outliers
    #loss(x, y) = sum(exp.(y .* 0.000001)) * Flux.mse(model(x), y)
    #loss(x, y) = maximum(10.0.^(y .* 0.0001)) * Flux.mse(model(x), y)

    accuracy(data) = evaluation1(data, model, statval, :RMSE)
    opt = Flux.ADAM()

    model, loss, accuracy, opt
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "rocm_cmake"
version = v"4.0.0"

# Collection of sources required to build
sources = [
    ArchiveSource("https://github.com/RadeonOpenCompute/rocm-cmake/archive/rocm-$(version).tar.gz",
                  "4577487acaa6e041a1316145867584f31caaf0d4aa2dd8fd7f82f81c269cada6"), # 4.0.0
]

# Bash recipe for building across all platforms
script = raw"""
cd ${WORKSPACE}/srcdir/rocm-cmake*/

mkdir build && cd build
cmake -DCMAKE_INSTALL_PREFIX=${prefix} ..
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Platform("x86_64", "linux"; libc="glibc"),
    Platform("x86_64", "linux"; libc="musl"),
]
platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
    FileProduct("share/rocm/cmake", :cmake_dir),
]

# Dependencies that must be installed before this package can be built
dependencies = [
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
# ---
# title: 147. Insertion Sort List
# id: problem147
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Linked List, Sort
# link: <https://leetcode.com/problems/insertion-sort-list/description/>
# hidden: true
# ---
#
# Sort a linked list using insertion sort.
#
# ![](https://upload.wikimedia.org/wikipedia/commons/0/0f/Insertion-sort-
# example-300px.gif)
# A graphical example of insertion sort. The partial sorted list (black)
# initially contains only the first element in the list.
# With each iteration one element (red) is removed from the input data and
# inserted in-place into the sorted list
#
#
# **Algorithm of Insertion Sort:**
#
#   1. Insertion sort iterates, consuming one input element each repetition, and growing a sorted output list.
#   2. At each iteration, insertion sort removes one element from the input data, finds the location it belongs within the sorted list, and inserts it there.
#   3. It repeats until no input elements remain.
#
#
# **Example 1:**
#
#
#
#     Input: 4->2->1->3
#     Output: 1->2->3->4
#
#
# **Example 2:**
#
#
#
#     Input: -1->5->3->4->0
#     Output: -1->0->3->4->5
#
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
function clean(phone_number)
    new = replace(phone_number, r"[^0-9]"  => "")
    noone = new[1] == '1' ? new[2:end] : new
    if length(noone) != 10
        return nothing
    elseif '1' ∈ [noone[1], noone[4]]
        return nothing
    elseif '0' ∈ [noone[1], noone[4]]
        return nothing
    else
        noone
    end
end

max_stars_repo_name: Rigante-pl/Exercism-Julia

max_issues_repo_name: Rigante-pl/Exercism-Julia

max_forks_repo_name: Rigante-pl/Exercism-Julia


text:
module MPI

include("compile-time.jl")
include("mpi-base.jl")

end


text:
#################################################################################
# Type Alias for long type signatures
const MinMax = NamedTuple{(:min, :max), NTuple{2, Float64}}
const NamedMinMax = Tuple{String, MinMax}
const UpDown = NamedTuple{(:up, :down), NTuple{2, Float64}}
const InOut = NamedTuple{(:in, :out), NTuple{2, Float64}}
const StartUpStages = NamedTuple{(:hot, :warm, :cold), NTuple{3, Float64}}

const BUILD_PROBLEMS_TIMER = TimerOutputs.TimerOutput()
const RUN_SIMULATION_TIMER = TimerOutputs.TimerOutput()

# Type Alias for JuMP and PJ containers
const JuMPExpressionMatrix = Matrix{<:JuMP.AbstractJuMPScalar}
const PGAE = PJ.ParametrizedGenericAffExpr{Float64, JuMP.VariableRef}
const GAE = JuMP.GenericAffExpr{Float64, JuMP.VariableRef}
const JuMPAffineExpressionArray = Matrix{GAE}
const JuMPAffineExpressionVector = Vector{GAE}
const JuMPConstraintArray = JuMP.Containers.DenseAxisArray{JuMP.ConstraintRef}
const JuMPVariableArray = JuMP.Containers.DenseAxisArray{JuMP.VariableRef}
const JuMPParamArray = JuMP.Containers.DenseAxisArray{PJ.ParameterRef}
const DenseAxisArrayContainer = Dict{Symbol, JuMP.Containers.DenseAxisArray}

# Settings constants
const UNSET_HORIZON = 0
const UNSET_INI_TIME = Dates.DateTime(0)

# Tolerance of comparisons
# MIP gap tolerances in most solvers are set to 1e-4
const ABSOLUTE_TOLERANCE = 1.0e-3
const BALANCE_SLACK_COST = 1e6
const SERVICES_SLACK_COST = 1e5
const COST_EPSILON = 1e-3
const MISSING_INITIAL_CONDITIONS_TIME_COUNT = 999.0
const SECONDS_IN_MINUTE = 60.0
const MINUTES_IN_HOUR = 60.0
const SECONDS_IN_HOUR = 3600.0
const MAX_START_STAGES = 3
const OBJECTIVE_FUNCTION_POSITIVE = 1.0
const OBJECTIVE_FUNCTION_NEGATIVE = -1.0
# The DEFAULT_RESERVE_COST value is used to avoid degeneracy of the solutions, reseve cost isn't provided.
const DEFAULT_RESERVE_COST = 1.0e-4
const KiB = 1024
const MiB = KiB * KiB
const GiB = MiB * KiB

const UNSUPPORTED_POWERMODELS =
    [PM.SOCBFPowerModel, PM.SOCBFConicPowerModel, PM.IVRPowerModel]

const PSI_NAME_DELIMITER = "__"

const M_VALUE = 1e6

const NO_SERVICE_NAME_PROVIDED = ""

# File Names definitions
const PROBLEM_SERIALIZATION_FILENAME = "operations_problem.bin"
const PROBLEM_BUILD_LOG_FILENAME = "operations_problem_build.log"
const HASH_FILENAME = "check.sha256"
const SIMULATION_SERIALIZATION_FILENAME = "simulation.bin"
const SIMULATION_LOG_FILENAME = "simulation.log"
const REQUIRED_RECORDERS = (:simulation_status, :simulation)
const KNOWN_SIMULATION_PATHS = [
    "data_store",
    "logs",
    "models_json",
    "problems",
    "recorder",
    "results",
    "simulation_files",
]

# Enums
IS.@scoped_enum(BuildStatus, IN_PROGRESS = -1, BUILT = 0, FAILED = 1, EMPTY = 2,)
IS.@scoped_enum(RunStatus, READY = -1, SUCCESSFUL = 0, RUNNING = 1, FAILED = 2,)
IS.@scoped_enum(SOSStatusVariable, NO_VARIABLE = 1, PARAMETER = 2, VARIABLE = 3,)

const ENUMS = (BuildStatus, RunStatus, SOSStatusVariable)

const ENUM_MAPPINGS = Dict()

for enum in ENUMS
    ENUM_MAPPINGS[enum] = Dict()
    for value in instances(enum)
        ENUM_MAPPINGS[enum][lowercase(string(value))] = value
    end
end

"""Get the enum value for the string. Case insensitive."""
function get_enum_value(enum, value::String)
    if !haskey(ENUM_MAPPINGS, enum)
        throw(ArgumentError("enum=$enum is not valid"))
    end

    val = lowercase(value)
    if !haskey(ENUM_MAPPINGS[enum], val)
        throw(ArgumentError("enum=$enum does not have value=$val"))
    end

    return ENUM_MAPPINGS[enum][val]
end

Base.convert(::Type{BuildStatus}, val::String) = get_enum_value(BuildStatus, val)
Base.convert(::Type{RunStatus}, val::String) = get_enum_value(RunStatus, val)
Base.convert(::Type{SOSStatusVariable}, x::String) = get_enum_value(SOSStatusVariable, x)


text:
import Base: *

"""
    *(A::AbstractOperator, ex::AbstractExpression)

Multiply an 'AbstractExpression` by an `AbstractOperator`.

# Example

```julia
julia> A = FiniteDiff(Float64, (10,))
δx  ℝ^10 -> ℝ^9

julia> x = Variable(10);

julia> ex = A*x;

julia> B = DCT(9)
ℱc  ℝ^9 -> ℝ^9

julia> ex2 = B*ex;

julia> affine(ex2)
ℱc*δx  ℝ^10 -> ℝ^9
```

"""
function (*)(L::AbstractOperator, a::AbstractExpression)
  A = convert(Expression,a)
  Expression{length(A.x)}(A.x,L*affine(A))
end

"""
    *(A::AbstractMatrix, ex::AbstractExpression)

Multiply an `AbstractExpression` by an `AbstractMatrix`.

```julia
julia> A = randn(10,5);

julia> x = Variable(5)
Variable(Float64, (5,))

julia> A*x

```

Other types of multiplications are also possible:

* left array multiplication
```julia
julia> X = Variable(10,5)
Variable(Float64, (10, 5))

julia> X*randn(5,6)

```

* scalar multiplication:
```julia
julia> π*X

```

* elementwise multiplication:
```julia
julia> randn(10,5).*X

```

"""
function (*)(m::T, a::Union{AbstractVector,AbstractMatrix}) where {T<:AbstractExpression}
  M = convert(Expression,m)
  op = LMatrixOp(codomainType(affine(M)),size(affine(M),1),a)
  return op*M
end
#LMatrixOp

function (*)(M::AbstractMatrix, a::T) where {T<:AbstractExpression}
  A = convert(Expression,a)
  op = MatrixOp(codomainType(affine(A)),size(affine(A),1),M)
  return op*A
end
#MatrixOp

function Broadcast.broadcasted(::typeof(*), d::D, a::T) where {D <: Union{Number,AbstractArray}, T<:AbstractExpression}
  A = convert(Expression,a)
  op = DiagOp(codomainType(affine(A)),size(affine(A),1),d)
  return op*A
end
Broadcast.broadcasted(::typeof(*), a::T, d::D) where {D <: Union{Number,AbstractArray}, T<:AbstractExpression} =
d.*a
#DiagOp

function (*)(coeff::T1, a::T) where {T1<:Number, T<:AbstractExpression}
  A = convert(Expression,a)
  return Expression{length(A.x)}(A.x,coeff*affine(A))
end
(*)(a::T, coeff::T1) where {T1<:Number, T<:AbstractExpression} = coeff*a
##Scale

"""
    *(A::AbstractExpression, ex::AbstractExpression)

Multiply an `AbstractExpression` by another `AbstractExpression`.

# Examples

```julia
julia> W1 = Variable(10,5)
Variable(Float64, (10, 5))

julia> W2 = Variable(5,15)
Variable(Float64, (5, 15))

julia> ex = W1*σ(W2);

julia> affine(ex)
I*σ  ℝ^(10, 5)  ℝ^(5, 15) -> ℝ^(10, 15)

```

`.*(A::AbstractExpression, ex::AbstractExpression)`

Elementwise multiplication between `AbstractExpression` (i.e. Hadamard product).

"""
function (*)(ex1::AbstractExpression, ex2::AbstractExpression)
  ex1 = convert(Expression,ex1)
  ex2 = convert(Expression,ex2)
  x = extract_variables((ex1,ex2))
  A = extract_affines(x, ex1)
  B = extract_affines(x, ex2)
  op = Ax_mul_Bx(A,B)
  exp3 = Expression{length(x)}(x,op)
  return exp3
end
# Ax_mul_Bx

function (*)(ex1::AdjointExpression, ex2::AbstractExpression)
  ex1 = ex1.ex
  ex2 = convert(Expression,ex2)
  x = extract_variables((ex1,ex2))
  A = extract_affines(x, ex1)
  B = extract_affines(x, ex2)
  op = Axt_mul_Bx(A,B)
  exp3 = Expression{length(x)}(x,op)
  return exp3
end
# Axt_mul_Bx

function (*)(ex1::AbstractExpression, ex2::AdjointExpression)
  ex1 = convert(Expression,ex1)
  ex2 = ex2.ex
  x = extract_variables((ex1,ex2))
  A = extract_affines(x, ex1)
  B = extract_affines(x, ex2)
  op = Ax_mul_Bxt(A,B)
  exp3 = Expression{length(x)}(x,op)
  return exp3
end
# Ax_mul_Bxt

function Broadcast.broadcasted(::typeof(*), ex1::AbstractExpression, ex2::AbstractExpression)
  ex1 = convert(Expression,ex1)
  ex2 = convert(Expression,ex2)
  x = extract_variables((ex1,ex2))
  A = extract_affines(x, ex1)
  B = extract_affines(x, ex2)
  op = HadamardProd(A,B)
  exp3 = Expression{length(x)}(x,op)
  return exp3
end
# Hadamard


text:
using SpatialFields
using Base.Test

function hrbf_2d()
	points = SVector{2, Float64}[[1; 0], [0; 1], [-1; 0], [0; -1]]
	normals = SVector{2, Float64}[[1; 1], [0; 1], [-1; 1], [0; -1]]

	field = HermiteRadialField(points, normals)

	X = linspace(-2, 2)
	Y = linspace(-2, 2)
	Z = zeros(length(X), length(Y))

	for i = 1:length(X)
	    for j = 1:length(Y)
	        Z[j,i] = field([X[i], Y[j]])
	    end
	end

	for i in 1:size(points, 2)
		@test isapprox(field(points[i]), 0, atol=1e-6)
		eps = 1e-4
		nudged_point = points[i] + eps * SVector(normals[i])
		@test isapprox(dot(SpatialFields.gradient(field)(points[i]), SVector(normals[i]) * eps), field(nudged_point), atol=1e-6)
	end

	g = SpatialFields.gradient(field)
	@test isapprox(g(SVector(1.,0)), [1;1])
end

hrbf_2d()


text:
## Prototype for generic control deflection on panel mesh
#==========================================================================================#

abstract type AbstractControlSurface end

struct BoundingBox{T <: Real} <: AbstractControlSurface
    location :: SVector{3,T}
    length   :: T
    width    :: T
    function BoundingBox{T}(x_o, s_r, c_r) where T <: Real
        @assert all(0. < x < 1., x_o) "Origin of bounding box must lie within [0,1] ⊂ ℝ!"
        @assert 0. <= s_r <= 1. "The normalized span length must lie within [0,1] ⊂ ℝ!"
        @assert 0. <= c_r <= 1. "The normalized chord length must lie within [0,1] ⊂ ℝ!"

        new{T}(x_o, s_r, c_r)
    end
end


text:
using NirsGadgets
using Test

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


text:
"""
`skipwrite(fpath, df; kwargs...)` use `CSV.write` to write the dataframe `df` to csv, but skip and warn when file exists (`isfile(fpath)` is true). It takes the same keyword arguments as `CSV.write`.

"""
function skipwrite(fpath, df; args...)
    fname = basename(fpath)
    if isfile(fpath)
        @warn "$fname already exists. Skip the file."
    else
        CSV.write(fpath, df; args...)
        println("""File "$fname" created. """)
    end
end


text:
function part1(input)
    x1, x2, y1, y2 = parse.(Int, match(r"target area: x=([-\d]+)\.\.([-\d]+), y=([-\d]+)\.\.([-\d]+)", readchomp(input)).captures)
    @assert y1 < 0 && x1 > 0
    @assert any(x1 <= sum(1:t) <= x2 for t in 1:-y1)
    return sum(1:(-1 - y1))
end

function part2(input)
    x1, x2, y1, y2 = parse.(Int, match(r"target area: x=([-\d]+)\.\.([-\d]+), y=([-\d]+)\.\.([-\d]+)", readchomp(input)).captures)
    X = Dict{Int, Set{Int}}()
    Y = Dict{Int, Set{Int}}()
    for x = 1:x2
        X[x] = find_x(x, x1, x2, 5 - 2 * y1)
    end
    for y = y1:(-y1)
        Y[y] = find_y(y, y1, y2)
    end
    t = 0
    Z = Set{Tuple{Int, Int}}()
    XX = Int[]
    YY = Int[]
    while true
        t += 1
        empty!(XX)
        empty!(YY)
        for y in collect(keys(Y))
            if t in Y[y]
                push!(YY, y)
                pop!(Y[y], t)
            end
            if isempty(Y[y])
                delete!(Y, y)
            end
        end
        isempty(YY) && continue
        for x in collect(keys(X))
            if t in X[x]
                push!(XX, x)
                pop!(X[x], t)
            end
            if isempty(X[x])
                delete!(X, x)
            end
        end

        for y in YY, x in XX
            push!(Z, (x, y))
        end
        if isempty(X) || isempty(Y)
            return length(Z)
        end
    end
end

function find_x(v, x1, x2, max_t)
    T = Set{Int}()
    t = 0
    x = 0
    while true
        t += 1
        x += v
        v = max(0, v - 1)
        if x1 <= x <= x2
            push!(T, t)
        end
        if x > x2 || t > max_t
            return T
        end
    end
end

function find_y(v, y1, y2)
    T = Set{Int}()
    t = 0
    y = 0
    while true
        t += 1
        y += v
        v -= 1
        if y1 <= y <= y2
            push!(T, t)
        end
        if y < y1
            return T
        end
    end
end


text:
import LsqFit
import Distributions
import PyPlot.plt

#Fit result type
immutable JLabFitresult
    #Fit inputs
    xdata::Array{Float64,1}     #independent variable
    ydata::Array{Float64,1}     #dependent variable
    yerrors::Array{Float64,1}   #uncertainties in dependent variable
    model::Function             #model to be fitted against

    #Fit results
    residue::Array{Float64,1}   #weighted residue
    params::Array{Float64,1}    #best fit parameters for model
    perrors::Array{Float64,1}   #uncertainties in best fit parameters
    dof::Int64                  #number of degrees of freedom in fit
    chisq::Float64              #chi^2 value of fit
    probability::Float64        #probability of chi^2 value
end

#Fitting procedure
function getfit{T1<:Real,T2<:Real,T3<:Real}(model::Function,xdata::Array{T1,1},ydata::Array{T2,1},initparams::Array{T3,1};kwargs...)
    return getfit(model,xdata,ydata,ones(Float64,length(ydata)),initparams;kwargs...)
end

function getfit{T1<:Real,T2<:Real,T3<:Real,T4<:Real}(model::Function,xdata::Array{T1,1},ydata::Array{T2,1},yerrors::Array{T3,1},initparams::Array{T4,1};alpha=0.6827)
    #Takes a model function with experimental data and returns a fit result.
    #model is a function f(x,p) with x an array of values for the independent
    #variable, and p an array of parameters for the model, and returns values for the
    #dependent variable.
    #
    #Usage:
    #   fit = getfit(model,xdata,ydata,yerrors,initparams)
    #       model           Function on (x,params)
    #       xdata           independent variable values
    #       ydata           dependent variable values
    #       yerrors         uncertainty in dependent variable
    #       initparams      initial guess for model parameters

    #Conduct fit
    weights = 1./yerrors
    fit = LsqFit.curve_fit(model,float64(xdata),float64(ydata),float64(weights),float64(initparams))
    perrors = LsqFit.estimate_errors(fit,alpha)

    chisq = sum(fit.resid.^2)
    chisqdist = Distributions.Chisq(fit.dof)
    probability = Distributions.ccdf(chisqdist,chisq)

    return JLabFitresult(float64(xdata),float64(ydata),float64(yerrors),model,fit.resid,fit.param,perrors,fit.dof,chisq,probability)
end

#Plot fit result
function plotfit(fit::JLabFitresult;datakwargs=(),modelkwargs=())
    #Makes a plot of data points and fit curve with annotations on fit parameters

    #Plot experimental data
    datapts = plt.errorbar(fit.xdata,fit.ydata,fit.yerrors,fmt=".",label="Expt data";datakwargs...)

    #Plot Model
    x = linspace(minimum(fit.xdata),maximum(fit.xdata),200)
    y = fit.model(x,fit.params)
    modelline = plt.plot(x,y,"r",lw=2,label="Model";modelkwargs...)

    #Annotations
    println("@sprintf(\"\\\$\\\\chi^2\\\$ / ndf = %.2f / %d = %.2f\\np(\\\$\\\\chi^2\\\$) = %.3f\\n\",fit.chisq,fit.dof,fit.chisq/fit.dof,fit.probability)")
    @printf("\$\\chi^2\$ / ndf = %.2f / %d = %.2f\np(\$\\chi^2\$) = %.3f\n\n",fit.chisq,fit.dof,fit.chisq/fit.dof,fit.probability)
    annot = plt.text(0.5*(minimum(fit.xdata)+maximum(fit.xdata)),0.5*(minimum(fit.ydata)+maximum(fit.ydata)),@sprintf("\$\\chi^2\$ / ndf = %.2f / %d = %.2f\np(\$\\chi^2\$) = %.3f\n\n",fit.chisq,fit.dof,fit.chisq/fit.dof,fit.probability))

    #Plot
    plt.xlabel("\$x\$")
    plt.ylabel("\$y\$")
    plt.title("\$y\$ vs \$x\$")
    plt.legend()

    return datapts,modelline,annot
end

# #############################################################################
#Common fitting functions

#Linear
linear{T<:Real,N}(x::Array{T,N},p::Array{Float64,1}) = p[1]+p[2]*x

#Polynomial
polynomial{T<:Real,N}(x::Array{T,N},p::Array{Float64,1}) = broadcast(^,x,transpose([0:length(p)-1]))*p

#Gaussian
gaussian{T<:Real,N}(x::Array{T,N},p::Array{Float64,1}=[1/sqrt(2*pi),1.,0.]) = p[1]*exp(-0.5*(x-p[3]).^2/p[2]^2)

#Lorentzian
lorentzian{T<:Real,N}(x::Array{T,N},p::Array{Float64,1}=[2/pi,1.,0.]) = p[1]/(1+(x-p[3]).^2/p[2]^2)

#Exponential decay
decay{T<:Real,N}(x::Array{T,N},p::Array{Float64,1}=[1.,1.]) = p[1]*exp(-p[2]*x)


text:
function cep_coeffs(n::Order,q::Number,N::Integer)
    iszero(q) && return q0_coeffs1(n,q,N)
    C1 = mat_C1(q,N) # Generate matrix
    Au = coeffs(C1,n)
    Au[1,:] .*= √0.5 # Renormalize eigenvectors
    return Au
end

function cea_coeffs(n::Order,q::Number,N::Integer)
    iszero(q) && return q0_coeffs2(n,q,N)
    C2 = mat_C2(q,N) # Generate matrix
    return coeffs(C2,n)
end

function sea_coeffs(n::Order,q::Number,N::Integer)
    iszero(q) && return q0_coeffs2(n,q,N)
    C3 = mat_C3(q,N) # Generate matrix
    return coeffs(C3,n)
end

function sep_coeffs(n::Order,q::Number,N::Integer)
    iszero(q) && return q0_coeffs2(n,q,N)
    C4 = mat_C4(q,N) # Generate matrix
    return coeffs(C4,n)
end

function coeffs(C::AbstractMatrix,n::Order)
    u = sort!(unique(n))
    au = eigvals(C)[u.+1] # Calculate characteristic values
    #au = LAPACK.stebz!('A','E',0.0,0.0,0,0,2*eps(C.ev[1]),C.dv,C.ev)[1][u+1]
    Au = eigvecs(C,au) # Calculate eigenvectors
    Au .*= signp.(sum(Au,dims=1)) # Make sure we have the correct sign
    return Au[:,indexin(checkvec(n),u)]
end

function q0_coeffs1(n::Order,q::Number,N::Integer)
    A = zeros(typeof(blasfloat(q)),N,length(n))
    for (i,n) in enumerate(n); A[n+1,i] = n==0 ? 1/√(eltype(A)(2)) : 1; end
    return A
end

function q0_coeffs2(n::Order,q::Number,N::Integer)
    A = zeros(typeof(blasfloat(q)),N,length(n))
    for (i,n) in enumerate(n); A[n+1,i] = one(eltype(A)); end
    return A
end

function per_coeffs(m::Order,q::Number,N::Integer)
    ie = map(iseven,m)
    all(ie)  && return cep_coeffs(div.(m,2),q,N)
    !any(ie) && return sep_coeffs(div.(m,2),q,N)
    P1 = cep_coeffs(div.(filter(iseven,m),2),q,N)
    P2 = sep_coeffs(div.(filter(isodd, m),2),q,N)
    P = similar(P1,N,length(m))
    P[:,ie] = P1
    P[:,.!ie] = P2
    return P
end


text:
#unique!(push!(LOAD_PATH, "~/Documents/UGent/PhD/Code_cMPS/julia1.0/CMPSKit.jl/src/"))
unique!(push!(LOAD_PATH, joinpath(pwd(), "src")))

using Revise
using CMPSKit
using KrylovKit
using OptimKit
using LinearAlgebra
using JLD2
using TensorOperations
using Plots

D = 4
k = 1.
g = 1.
Λ = 1000.

gradienttest = false
tensorproduct = true
lagrangianmultiplier = false
optimisation = true

if tensorproduct
    #Tensors product ansatz
    Id = 1*Matrix(I,D,D)
    Q = Constant(randn(D^2,D^2))
    R1 = Constant(kron(Id,randn(D,D)))
    R2 = Constant(kron(randn(D,D),Id))
    #Put them in cMPS form
    Ψ = InfiniteCMPS(Q, (R1,R2))
else
    #Random Tensors
    Q = Constant(randn(D^2,D^2))
    R1 = Constant(randn(D^2,D^2))
    R2 = Constant(randn(D^2,D^2))
    #Put them in cMPS form
    Ψ = InfiniteCMPS(Q, (R1,R2))
end

function firstordercorrection(dQ,gradQ,dR1,gradR1,dR2,gradR2)
    #left gauge: dQ = -R1'*dR1 - R2'*dR2
    return 2*real(dot(dQ,gradQ) + dot(dR1,gradR1) + dot(dR2,gradR2));
end

alg1 = LBFGS(; verbosity = 2, maxiter = 1000000, gradtol = 1e-4);
alg2 = ConjugateGradient(; verbosity = 2, maxiter = 1000000, gradtol = 1e-4);
if tensorproduct
    linalg = GMRES(krylovdim = 80; tol = 1e-5)
else
    linalg = GMRES(krylovdim = 50; tol = 1e-5)
end

if lagrangianmultiplier
    h = k * (∂ψ[1]'*∂ψ[1] + ∂ψ[2]'*∂ψ[2]) - 2 * g * (ψ[1]'*ψ[2]'*ψ[2]*ψ[1] + ψ[2]'*ψ[1]'*ψ[1]*ψ[2]) + Λ * ((ψ[1]')^2*ψ[1]^2 + (ψ[2]')^2*ψ[2]^2 + (ψ[1]*ψ[2] - ψ[2]*ψ[1])' * (ψ[1]*ψ[2] - ψ[2]*ψ[1]))
    H = ∫(h, (-Inf,+Inf))
else
    h = k * (∂ψ[1]'*∂ψ[1] + ∂ψ[2]'*∂ψ[2]) - 2 * g * (ψ[1]'*ψ[2]'*ψ[2]*ψ[1] + ψ[2]'*ψ[1]'*ψ[1]*ψ[2])
    H = ∫(h, (-Inf,+Inf))
end

#h = k * (∂ψ[1]'*∂ψ[1] + ∂ψ[2]'*∂ψ[2]) - μ * (ψ[1]'*ψ[1] + ψ[2]'*ψ[2]) + g * ((ψ[1]')^2*ψ[1]^2 + (ψ[2]')^2*ψ[2]^2) + Λ * ((ψ[1]*ψ[2] - ψ[2]*ψ[1])' * (ψ[1]*ψ[2] - ψ[2]*ψ[1]))
#H = ∫(k * (∂ψ[1]'*∂ψ[1] + ∂ψ[2]'*∂ψ[2]) - μ * (ψ[1]'*ψ[1] + ψ[2]'*ψ[2]) + g * ((ψ[1]')^2*ψ[1]^2 + (ψ[2]')^2*ψ[2]^2) + Λ * ((ψ[1]*ψ[2] - ψ[2]*ψ[1])' * (ψ[1]*ψ[2] - ψ[2]*ψ[1])), (-Inf,+Inf))
#h = k * (∂ψ[1]'*∂ψ[1] + ∂ψ[2]'*∂ψ[2]) - μ * (ψ[1]'*ψ[1] + ψ[2]'*ψ[2]) + g * ((ψ[1]')^2*ψ[1]^2 + (ψ[2]')^2*ψ[2]^2)
#H = ∫(k * (∂ψ[1]'*∂ψ[1] + ∂ψ[2]'*∂ψ[2]) - μ * (ψ[1]'*ψ[1] + ψ[2]'*ψ[2]) + g * ((ψ[1]')^2*ψ[1]^2 + (ψ[2]')^2*ψ[2]^2), (-Inf,+Inf))

if gradienttest
    ρL, ρR, infoL, infoR = environments!(Ψ;linalg = GMRES(; tol = 1e-6))

    HL, E, e, hL, infoL = leftenv(H, (Ψ,ρL,ρR);linalg = GMRES(; tol = 1e-6))
    HR, E, e, hR, infoR = rightenv(H, (Ψ,ρL,ρR);linalg = GMRES(; tol = 1e-6))

    gradQ, gradRs = gradient(H, (Ψ, ρL, ρR))
    dQ = Constant(1e-4*randn(D^2,D^2))
    dR1 = Constant(kron(Id,1e-4*(randn(D,D))))
    dR2 = Constant(kron(1e-4*(randn(D,D)),Id))
    #left gauge dK = -1im*(dR1'*R1 - R1'*dR1 + dR2'*R2 - R2'*dR2)/2
    Ψ2 = InfiniteCMPS(Q+dQ,(R1+dR1,R2+dR2))
    foc1 = firstordercorrection(dQ,gradQ,dR1,gradRs[1],dR2,gradRs[2])
    dQ/=10
    dR1/=10
    dR2/=10
    foc2 = firstordercorrection(dQ,gradQ,dR1,gradRs[1],dR2,gradRs[2])
    Ψ3 = InfiniteCMPS(Q+dQ,(R1+dR1,R2+dR2))
    dQ/=10
    dR1/=10
    dR2/=10
    foc3 = firstordercorrection(dQ,gradQ,dR1,gradRs[1],dR2,gradRs[2])
    Ψ4 = InfiniteCMPS(Q+dQ,(R1+dR1,R2+dR2))

    @show E1 = expval(h,Ψ)
    E2 = expval(h,Ψ2)
    E3 = expval(h,Ψ3)
    E4 = expval(h,Ψ4)

    @show E2[]-E1[]-foc1
    @show E3[]-E1[]-foc2
    @show E4[]-E1[]-foc3
end

if optimisation
    Ψ, ρR, E, e, normgrad, numfg, history = groundstate5(H, Ψ; optalg = alg1, linalg = linalg)
    #αs,fs, dfs1, dfs2 = groundstate(H, Ψ; optalg = alg1, linalg = GMRES(; tol = 1e-5))
    #αs = (αs[1:end-1] + αs[2:end])/2
    #push!(αs,0.1)

    @show Q = Ψ.Q
    @show R1 = Ψ.Rs[1]
    @show R2 = Ψ.Rs[2]
    @show R1[]*R2[] + R2[]*R1[]
    @show expval(ψ[1]*ψ[2] + ψ[2]*ψ[1],Ψ)[]
    #DR1 = differentiate(R1) + Q * R1 - R1 * Q
    #DR2 = differentiate(R2) + Q * R2 - R2 * Q
    #R1² = R1 * R1
    #R2² = R2 * R2
end


text:
"""
    function get_reaction_idx(rn::Catalyst.ReactionSystem)

Get the rearranged order for a reaction system, returns a matrix `m` with first column as the rearranged reactions, second column is the order in the orign `ReactionSystem`
If you know the old idx for a certain reaction in `ReactionSystem`, you can use
`m[:,2][idx]` to get the new idx.
"""
function get_reaction_idx(rn::Catalyst.ReactionSystem)
    massactionjump_order = Int64[]
    constantjump_order = Int64[]
    rxvars = []
    for (i,rx) in enumerate(ModelingToolkit.get_eqs(rn))
        empty!(rxvars)
        (rx.rate isa ModelingToolkit.Symbolic) && ModelingToolkit.get_variables!(rxvars, rx.rate)
        @inbounds for j = 1:length(rxvars)
            if isequal(rxvars[j], ModelingToolkit.get_iv(rn))
                error("Does not support VariableRateJump")
            end
        end
        if Catalyst.ismassaction(rx, rn)
            push!(massactionjump_order, i)
        else
            push!(constantjump_order, i)
        end
    end
    order = vcat(massactionjump_order, constantjump_order)
    # println("The rearranged order :")
    vcat([[rn.eqs[order[i]] order[i]] for i in eachindex(order)]...)
end

export get_reaction_idx


text:
# OpenCL.Platform

struct Platform <: CLObject
    id::CL_platform_id
end

Base.pointer(p::Platform) = p.id

function info(p::Platform, pinfo::Symbol)
    info_map = Dict{Symbol, CL_platform_info}(
        :profile => CL_PLATFORM_PROFILE,
        :version => CL_PLATFORM_VERSION,
        :name    => CL_PLATFORM_NAME,
        :vendor  => CL_PLATFORM_VENDOR,
        :extensions => CL_PLATFORM_EXTENSIONS
    )
    try
        cl_info = info_map[pinfo]
        inf = info(p, cl_info)
        pinfo == :extensions && return split(inf)
        return inf
    catch err
        if isa(err, KeyError)
            throw(ArgumentError("OpenCL.Platform has no info for: $pinfo"))
        else
            throw(err)
        end
    end
end

Base.getindex(p::Platform, pinfo::Symbol) = info(p, pinfo)

function Base.show(io::IO, p::Platform)
    strip_extra_whitespace = r"\s+"
    platform_name = replace(p[:name], strip_extra_whitespace => " ")
    ptr_val = convert(UInt, Base.pointer(p))
    ptr_address = "0x$(string(ptr_val, base = 16, pad = Sys.WORD_SIZE>>2))"
    print(io, "OpenCL.Platform('$platform_name' @$ptr_address)")
end

function platforms()
    nplatforms = Ref{CL_uint}()
    @check api.clGetPlatformIDs(0, C_NULL, nplatforms)
    cl_platform_ids = Vector{CL_platform_id}(undef, nplatforms[])
    @check api.clGetPlatformIDs(nplatforms[], cl_platform_ids, C_NULL)
    return [Platform(id) for id in cl_platform_ids]
end

function num_platforms()
    nplatforms = Ref{CL_uint}()
    @check api.clGetPlatformIDs(0, C_NULL, nplatforms)
    return Int(nplatforms[])
end

function info(p::Platform, pinfo::CL_platform_info)
    size = Ref{Csize_t}()
    @check api.clGetPlatformInfo(p.id, pinfo, 0, C_NULL, size)
    result = Vector{CL_char}(undef, size[])
    @check api.clGetPlatformInfo(p.id, pinfo, size[], result, C_NULL)
    return CLString(result)
end

function devices(p::Platform, dtype::CL_device_type)
    try
        ndevices = Ref{CL_uint}()
        @check api.clGetDeviceIDs(p.id, dtype, 0, C_NULL, ndevices)
        if ndevices[] == 0
            return Device[]
        end
        result = Vector{CL_device_id}(undef, ndevices[])
        @check api.clGetDeviceIDs(p.id, dtype, ndevices[], result, C_NULL)
        return Device[Device(id) for id in result]
    catch err
        if err.desc == :CL_DEVICE_NOT_FOUND || err.code == -1
            return Device[]
        else
            throw(err)
        end
    end
end

devices(p::Platform) = devices(p, CL_DEVICE_TYPE_ALL)

function devices(p::Platform, dtype::Symbol)
    devices(p, cl_device_type(dtype))
end

function devices(dtype::CL_device_type)
    devs = Device[]
    for platform in platforms()
        append!(devs, devices(platform, dtype))
    end
    return devs
end
devices(dtype::Symbol) = devices(cl_device_type(dtype))

function devices()
    devs = Device[]
    for platform in platforms()
        append!(devs, devices(platform))
    end
    return devs
end

has_device_type(p::Platform, dtype) = length(devices(p, dtype)) > 0

available_devices(p::Platform, dtype::Symbol) = filter(d -> info(d, :available),  devices(p, dtype))
available_devices(p::Platform) = available_devices(p, :all)


text:
import Statistics
using LinearAlgebra
using Cassette
using Cassette: Tagged, tag, untag, istagged, metadata, hasmetadata,
    enabletagging, @overdub, overdub, canrecurse, similarcontext, fallback
using JuliaInterpreter
using Espresso
using CUDAapi

include("utils.jl")
include("helpers.jl")
include("devices.jl")
include("tape.jl")
include("tapeutils.jl")
include("trace/trace.jl")
include("diffrules/diffrules.jl")
include("grad.jl")
include("compile.jl")
include("update.jl")
include("transform.jl")

# function __init__()
#     @require CuArrays="3a865a2d-5b23-5a0f-bc46-62713ec82fae" include("cuda.jl")
# end

const BEST_AVAILABLE_DEVICE = Ref{AbstractDevice}(CPU())

if has_cuda()
    try
        using CuArrays
        using CUDAnative

        BEST_AVAILABLE_DEVICE[] = GPU(0)

        include("cuda.jl")
    catch ex
        # something is wrong with the user's set-up (or there's a bug in CuArrays)
        @warn "CUDA is installed, but CuArrays.jl fails to load" exception=(ex,catch_backtrace())

    end
end

best_available_device() = BEST_AVAILABLE_DEVICE[]


text:
using CAPI

prob = CAPI.Problem([1,2.0], [1.0,2])
para = CAPI.Parameter(1.0,1.0)

pprob = convert(Ptr{CAPI.Problem}, pointer_from_objref(prob))
ppara = convert(Ptr{CAPI.Parameter}, pointer_from_objref(para))
CAPI.compute_score(pprob, ppara)


text:
"""
    create_vector(nrows, nrows_loc = PETSC_DECIDE)

Create a `PetscVec` vector of global size `(nrows)`.
"""
function create_vector(nrows, nrows_loc = PETSC_DECIDE)
    vec = VecCreate()
    VecSetSizes(vec::PetscVec, nrows_loc, nrows)
    return vec
end

Base.ndims(::Type{PetscVec}) = 1

"""
    Base.setindex!(vec::PetscVec, value::Number, row::Integer)

`row` must be in [1,size(vec)], i.e indexing starts at 1 (Julia).

# Implementation
For some unkwnown reason, calling `VecSetValue` fails.
"""
function Base.setindex!(vec::PetscVec, value::Number, row::Integer)
    VecSetValues(vec, PetscInt[row .- 1], PetscScalar[value], INSERT_VALUES)
end

# This is stupid but I don't know how to do better yet
Base.setindex!(vec::PetscVec, values, rows) = VecSetValues(vec, collect(rows .- 1), values, INSERT_VALUES)

set_global_size!(vec::PetscVec, nrows) = VecSetSizes(vec, PETSC_DECIDE, nrows)
set_local_size!(vec::PetscVec, nrows) = VecSetSizes(vec, nrows, PETSC_DECIDE)

set_from_options!(vec::PetscVec) = VecSetFromOptions(vec)

set_up!(vec::PetscVec) = VecSetUp(vec)

"""
    get_range(vec::PetscVec)

Wrapper to `VecGetOwnershipRange`

However, the result `(rstart, rend)` is such that `mat[rstart:rend]` are the rows handled by the local processor.
This is different from the default `PETSc.VecGetOwnershipRange` result where the indexing starts at zero and where
`rend-1` is last row handled by the local processor.
"""
function get_range(vec::PetscVec)
    rstart, rend = VecGetOwnershipRange(vec)
    return (rstart + 1, rend)
end

# Discutable choice
Base.length(vec::PetscVec) = VecGetLocalSize(vec)
Base.size(vec::PetscVec) = (length(vec),)

function assemble!(vec::PetscVec)
    VecAssemblyBegin(vec)
    VecAssemblyEnd(vec)
end

duplicate(vec::PetscVec) = VecDuplicate(vec)

"""
    vec2array(vec::PetscVec)

Convert a `PetscVec` into a Julia `Array`. Allocation is involved in the process since the `PetscVec`
allocated by PETSC is copied into a freshly allocated array. If you prefer not to allocate memory,
use `VectGetArray` and `VecRestoreArray`
"""
function vec2array(vec::PetscVec)
    arrayFromC, array_ref = VecGetArray(vec)
    array = copy(arrayFromC)
    VecRestoreArray(vec, array_ref)
    return array
end

Base.show(::IO, vec::PetscVec) = VecView(vec)

destroy!(vec::PetscVec) = VecDestroy(vec)


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

# # Logistic Regression

# **Originally Contributed by**: François Pacaud

# This tutorial shows how to solve a logistic regression problem
# with JuMP. Logistic regression is a well known method in machine learning,
# useful when we want to classify binary variables with the help of
# a given set of features. To this goal,
# we find the optimal combination of features maximizing
# the (log)-likelihood onto a training set. From a modern optimization glance,
# the resulting problem is convex and differentiable. On a modern optimization
# glance, it is even conic representable.
#
# ## Formulating the logistic regression problem
#
# Suppose we have a set of training data-point $i = 1, \cdots, n$, where
# for each $i$ we have a vector of features $x_i \in \mathbb{R}^p$ and a
# categorical observation $y_i \in \{-1, 1\}$.
#
# The log-likelihood is given by
#
# ```math
# l(\theta) = \sum_{i=1}^n \log(\dfrac{1}{1 + \exp(-y_i \theta^\top x_i)})
# ```
#
# and the optimal $\theta$ minimizes the logistic loss function:
#
# ```math
# \min_{\theta}\; \sum_{i=1}^n \log(1 + \exp(-y_i \theta^\top x_i)).
# ```
#
# Most of the time, instead of solving directly the previous optimization problem, we
# prefer to add a regularization term:
#
# ```math
# \min_{\theta}\; \sum_{i=1}^n \log(1 + \exp(-y_i \theta^\top x_i)) + \lambda \| \theta \|
# ```
#
# with $\lambda \in \mathbb{R}_+$ a penalty and $\|.\|$ a norm function. By adding
# such a regularization term, we avoid overfitting on the training set and usually
# achieve a greater score in cross-validation.

# ## Reformulation as a conic optimization problem
#
# By introducing auxiliary variables $t_1, \cdots, t_n$ and $r$,
# the optimization problem is equivalent to
#
# ```math
# \begin{aligned}
# \min_{t, r, \theta} \;& \sum_{i=1}^n t_i + \lambda r \\
# \text{subject to } & \quad t_i \geq \log(1 + \exp(- y_i \theta^\top x_i)) \\
#                    & \quad r \geq \|\theta\|
# \end{aligned}
# ```
#
# Now, the trick is to reformulate the constraints $t_i \geq \log(1 + \exp(- y_i \theta^\top x_i))$
# with the help of the *exponential cone*
#
# ```math
# K_{exp} = \{ (x, y, z) \in \mathbb{R}^3 : \; y \exp(x / y) \leq z \} .
# ```
#
# Indeed, by passing to the exponential, we
# see that for all $i=1, \cdots, n$, the constraint $t_i \geq \log(1 + \exp(- y_i \theta^\top x_i))$
# is equivalent to
#
# ```math
# \exp(-t_i) + \exp(u_i - t_i) \leq 1
# ```
#
# with $u_i = -y_i \theta^\top x_i$. Then, by adding two auxiliary variables
# $z_{i1}$ and $z_{i2}$ such that $z_{i1} \geq \exp(u_i-t_i)$ and $z_{i2} \geq \exp(-t_i)$, we get
# the equivalent formulation
#
# ```math
# \left\{
# \begin{aligned}
# (u_i -t_i , 1, z_{i1}) & \in  K_{exp}  \\
# (-t_i , 1, z_{i2}) & \in  K_{exp}  \\
# z_{i1} + z_{i2} & \leq  1
# \end{aligned}
# \right.
# ```
#
# In this setting, the conic version of the logistic regression problems writes out
#
# ```math
# \begin{aligned}
# \min_{t, z, r, \theta}&  \; \sum_{i=1}^n t_i + \lambda r \\
# \text{subject to } & \quad  (u_i -t_i , 1, z_{i1})  \in  K_{exp}  \\
#                    & \quad  (-t_i , 1, z_{i2})  \in  K_{exp}  \\
#                    & \quad  z_{i1} + z_{i2}  \leq  1 \\
#                    & \quad u_i = -y_i x_i^\top \theta \\
#                    & \quad r \geq \|\theta\|
# \end{aligned}
# ```
#
# and thus encompasses $3n + p + 1$ variables and $3n + 1$ constraints ($u_i = -y_i \theta^\top x_i$
# is only a virtual constraint used to clarify the notation).
# Thus, if $n \gg 1$, we get a large number of variables and constraints.

# ## Fitting logistic regression with a conic solver
#
# It is now time to pass to the implementation. We choose SCS as a conic solver.
using JuMP
import Random
import SCS

Random.seed!(2713);

# We start by implementing a function to generate a fake dataset, and where
# we could tune the correlation between the feature variables. The function
# is a direct transcription of the one used in [this blog post](http://fa.bianp.net/blog/2013/numerical-optimizers-for-logistic-regression/).
function generate_dataset(n_samples=100, n_features=10; shift=0.0)
    X = randn(n_samples, n_features)
    w = randn(n_features)
    y = sign.(X * w)
    X .+= 0.8 * randn(n_samples, n_features) # add noise
    X .+= shift # shift the points in the feature space
    X = hcat(X, ones(n_samples, 1))
    return X, y
end

# We write a `softplus` function to formulate each constraint
# $t \geq \log(1 + \exp(u))$ with two exponential cones.
function softplus(model, t, u)
    z = @variable(model, [1:2], lower_bound=0.0)
    @constraint(model, sum(z) <= 1.0)
    @constraint(model, [u - t, 1, z[1]] in MOI.ExponentialCone())
    @constraint(model, [-t, 1, z[2]] in MOI.ExponentialCone())
end

# ### $\ell_2$ regularized logistic regression
# Then, with the help of the `softplus` function, we could write our
# optimization model. In the $\ell_2$ regularization case, the constraint
# $r \geq \|\theta\|_2$ rewrites as a second order cone constraint.
function build_logit_model(X, y, λ)
    n, p = size(X)
    model = Model()
    @variable(model, θ[1:p])
    @variable(model, t[1:n])
    for i in 1:n
        u = - (X[i, :]' * θ) * y[i]
        softplus(model, t[i], u)
    end
    ## Add ℓ2 regularization
    @variable(model, 0.0 <= reg)
    @constraint(model, [reg; θ] in MOI.SecondOrderCone(p+1))
    ## Define objective
    @objective(model, Min, sum(t) + λ * reg)
    return model
end

# We generate the dataset.
#
# !!! warning
#     Be careful here, for large n and p SCS could fail to converge!
#
n, p = 200, 10
X, y = generate_dataset(n, p, shift=10.0);

## We could now solve the logistic regression problem
λ = 10.0
model = build_logit_model(X, y, λ)
set_optimizer(model, SCS.Optimizer)
JuMP.optimize!(model)

#-

θ♯ = JuMP.value.(model[:θ])

# It appears that the speed of convergence is not that impacted by the correlation
# of the dataset, nor by the penalty $\lambda$.

# ### $\ell_1$ regularized logistic regression
#
# We now formulate the logistic problem with a $\ell_1$ regularization term.
# The $\ell_1$ regularization ensures sparsity in the optimal
# solution of the resulting optimization problem. Luckily, the $\ell_1$ norm
# is implemented as a set in `MathOptInterface`. Thus, we could easily formulate
# the sparse logistic regression problem with the help of a `MOI.NormOneCone`
# set.
function build_sparse_logit_model(X, y, λ)
    n, p = size(X)
    model = Model()
    @variable(model, θ[1:p])
    @variable(model, t[1:n])
    for i in 1:n
        u = - (X[i, :]' * θ) * y[i]
        softplus(model, t[i], u)
    end
    ## Add ℓ1 regularization
    @variable(model, 0.0 <= reg)
    @constraint(model, [reg; θ] in MOI.NormOneCone(p+1))
    ## Define objective
    @objective(model, Min, sum(t) + λ * reg)
    return model
end

## Auxiliary function to count non-null components:
count_nonzero(v::Vector; tol=1e-6) = sum(abs.(v) .>= tol)

## We solve the sparse logistic regression problem on the same dataset as before.
λ = 10.0
sparse_model = build_sparse_logit_model(X, y, λ)
set_optimizer(sparse_model, SCS.Optimizer)
JuMP.optimize!(sparse_model)

#-

θ♯ = JuMP.value.(sparse_model[:θ])
println("Number of non-zero components: ", count_nonzero(θ♯),
        " (out of ", p, " features)")

# ### Extensions
# A direct extension would be to consider the sparse logistic regression with
# *hard* thresholding, which, on contrary to the *soft* version using a $\ell_1$ regularization,
# adds an explicit cardinality constraint in its formulation:
#
# ```math
# \begin{aligned}
# \min_{\theta} & \; \sum_{i=1}^n \log(1 + \exp(-y_i \theta^\top x_i)) + \lambda \| \theta \|_2^2 \\
# \text{subject to } & \quad \| \theta \|_0 <= k
# \end{aligned}
# ```
#
# where $k$ is the maximum number of non-zero components in the vector $\theta$,
# and $\|.\|_0$ is the $\ell_0$ pseudo-norm:
#
# ```math
# \| x\|_0 = \#\{i : \; x_i \neq 0\}
# ```
#
# The cardinality constraint $\|\theta\|_0 \leq k$ could be reformulated with
# binary variables. Thus the hard sparse regression problem could be solved
# by any solver supporting mixed integer conic problems.
#

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
"""
    IBooleanValued(nexid)

An interface for boolean valued widgets.
"""
struct IBooleanValued <: AbstractINextion
    _nid::NexID

    function IBooleanValued(nid::NexID)
        new(nid)
    end
end

"""
    obj.value = val

Set boolean value from `val` to Nextion object `obj`.
"""
function setproperty!(obj::IBooleanValued, property::Symbol, new_val)

    if property == :value
        setnexproperty!(NexID(obj), :val, Bool(new_val))
    else
        setfield!(obj, property, new_val)
    end

end

"""
    obj.value

Get boolean value from Nextion object `obj`.
"""
function getproperty(obj::IBooleanValued, property::Symbol)

    if property == :value
        Bool(getnexproperty(NexID(obj), :val, Int))
    else
        getfield(obj, property)
    end

end


text:
# ps3000_overview_buffer_ps_status.jl

export ps3000_overview_buffer_status

function ps3000_overview_buffer_status(handle :: Int16)
  # returns true if buffer has overrun
  previous_buffer_overrun = Ref{Int16}()
  ps_status = ccall((:ps3000_overview_buffer_status, ps3000driver),Int16,
    (Int16, Ref{Int16}), handle, previous_buffer_overrun)
  if ps_status == 0
    error("ps3000 Error: invalid handle")
  end
  return previous_buffer_overrun[]==0?false:true
end


text:
module GdMn2O5
    using InlineExports
    using Reexport
    @reexport using Overseer
    @reexport using Optim
    @reexport using StaticArrays
    using LinearAlgebra
    @reexport using DrWatson
    using BlackBoxOptim
    using Unitful
    using Unitful: eV, μ0, μB, Å
    using DataFrames
    using MuladdMacro
    using FileIO
    using Parameters
    using ProgressMeter
    using Requires

    const Vec2 = SVector{2, Float64}
    Vec2(x::Number) = Vec2(x, x)

    const Vec3 = SVector{3, Float64}
    Vec3(x::Number) = Vec3(x, x, x)

    polar2vec(ϕ) = Vec2(cos(ϕ), sin(ϕ))
    polar2vec(ϕ, z) = Vec3(cos(ϕ), sin(ϕ), z)

    include("components.jl")
    export Spin, L, Pb, EasyAxis, Optimize, Position, H, E_field, Model, Etot, Bond, JBond

    include("systems.jl")
    include("energies.jl")
    include("ledger.jl")
    include("optimize.jl")
    include("plotting.jl")

    @export H_sweep_range(start, stop, len=100) =
        vcat(collect(repeat([range(start, stop, length=len);
                             range(stop, start, length=len)[2:end-1]],
                            2)),
             range(start, stop, length=len))

     @export function findclosest(val, options)
         min = typemax(eltype(options))
         minval = zero(eltype(options))
         imin = 0
         for (i, o) in enumerate(options)
             t = abs(o - val)
             if t < min
                 min = t
                 imin = i
                 minval = o
             end
         end
         return minval, imin
    end

    @export function winding_number(hsweep)
        state2_id = Int(floor(length(hsweep)/4))
        if abs(hsweep[1][L][1].ϕ - hsweep[end][L][1].ϕ) > 1.9π || abs(hsweep[1][L][2].ϕ - hsweep[end][L][2].ϕ) > 1.9π
            return 1
        elseif sign(Pb(hsweep[1])) ==  sign(Pb(hsweep[state2_id]))
            return 1.5
        else
            return 0.5
        end
    end
    function __init__()
        @require Glimpse="f6e19d58-12a4-5927-8606-ac30a9ce9b69" include("visualization.jl")
    end

end


text:
###
### A short immutable sequence type for representing Kmers and Skipmers.
###

abstract type AbstractMer{A<:NucleicAcidAlphabet{2},K} <: BioSequence{A} end

primitive type Mer{A<:NucleicAcidAlphabet{2},K} <: AbstractMer{A,K} 64 end
primitive type BigMer{A<:NucleicAcidAlphabet{2},K} <: AbstractMer{A,K} 128 end

# Aliases
const DNAMer{K} = Mer{DNAAlphabet{2},K}
const RNAMer{K} = Mer{RNAAlphabet{2},K}
const DNAKmer   = DNAMer{31}
const RNAKmer   = RNAMer{31}

const BigDNAMer{K} = BigMer{DNAAlphabet{2},K}
const BigRNAMer{K} = BigMer{RNAAlphabet{2},K}
const BigDNAKmer   = BigMer{DNAAlphabet{2},63}
const BigRNAKmer   = BigMer{RNAAlphabet{2},63}

const DNACodon = DNAMer{3}
const RNACodon = RNAMer{3}

###
### Base Functions
###

@inline encoded_data_type(::Type{Mer{A,K}}) where {A,K} = UInt64
@inline encoded_data_type(::Type{BigMer{A,K}}) where {A,K} = UInt128
@inline encoded_data_type(x::AbstractMer) = encoded_data_type(typeof(x))
@inline encoded_data(x::AbstractMer) = reinterpret(encoded_data_type(typeof(x)), x)
@inline capacity(x::Type{U}) where {U<:Unsigned} = div(8 * sizeof(U), 2)
@inline capacity(mer::Type{T}) where {T<:AbstractMer} = capacity(encoded_data_type(mer))
@inline capacity(mer::AbstractMer) = capacity(typeof(mer))
@inline ksize(::Type{T}) where {A,K,T<:AbstractMer{A,K}} = K
@inline n_unused(::Type{T}) where {T<:AbstractMer} = capacity(T) - ksize(T)
@inline n_unused(x::AbstractMer) = n_unused(typeof(x))

@inline function checkmer(::Type{T}) where {T<:AbstractMer}
    c = capacity(T)
    if !(1 ≤ ksize(T) ≤ c)
        throw(ArgumentError("K must be within 1..$c"))
    end
end

@inline Base.length(x::AbstractMer{A,K}) where {A,K} = ksize(typeof(x))
@inline Base.unsigned(x::AbstractMer) = encoded_data(x)
@inline Base.summary(x::AbstractMer{DNAAlphabet{2},K}) where {K} = string("DNA ", K, "-mer")
@inline Base.summary(x::AbstractMer{RNAAlphabet{2},K}) where {K} = string("RNA ", K, "-mer")

function Base.typemin(::Type{T}) where {T<:AbstractMer}
    checkmer(T)
    return T(typemin(encoded_data_type(T)))
end

function Base.typemax(::Type{T}) where {T<:AbstractMer}
    checkmer(T)
    U = encoded_data_type(T)
    return T(typemax(U) >> (8 * sizeof(U) - 2ksize(T)))
end

function Base.rand(::Type{T}) where {T<:AbstractMer}
    return T(rand(encoded_data_type(T)))
end

function Base.rand(::Type{T}, size::Integer) where {T<:AbstractMer}
    return [rand(T) for _ in 1:size]
end

Base.:-(x::AbstractMer, y::Integer) = typeof(x)(encoded_data(x) - y % encoded_data_type(x))
Base.:+(x::AbstractMer, y::Integer) = typeof(x)(encoded_data(x) + y % encoded_data_type(x))
Base.:+(x::AbstractMer, y::AbstractMer) = y + x

Alphabet(::Type{Mer{A,K} where A<:NucleicAcidAlphabet{2}}) where {K} = Any

###
### Conversion
###

function Mer{A,K}(x::UInt64) where {A,K}
    checkmer(Mer{A,K})
    mask = (one(UInt64) << (2 * K)) - 1
    return reinterpret(Mer{A,K}, x & mask)
end

function BigMer{A,K}(x::UInt128) where {A,K}
    checkmer(BigMer{A,K})
    mask = (one(UInt128) << (2 * K)) - 1
    return reinterpret(BigMer{A,K}, x & mask)
end

include("conversion.jl")
include("indexing.jl")
include("predicates.jl")
include("counting.jl")
include("transformations.jl")

###
### Mer de-bruijn neighbors
###

# Neighbors on a de Bruijn graph
struct MerNeighborIterator{S<:AbstractMer}
    x::S
end

"""
    neighbors(skipmer::S) where {S<:AbstractMer}

Return an iterator through skip-mers neighboring `skipmer` on a de Bruijn graph.
"""
neighbors(skipmer::AbstractMer) = MerNeighborIterator{typeof(skipmer)}(skipmer)

Base.length(::MerNeighborIterator) = 4
Base.eltype(::Type{MerNeighborIterator{S}}) where {S<:AbstractMer} = S

function Base.iterate(it::MerNeighborIterator{S}, i = zero(encoded_data_type(S))) where {S<:AbstractMer}
    if i == 4
        return nothing
    else
        return S((encoded_data(it.x) << 2) | i), i + 1
    end
end

###
### String literals
###

macro mer_str(seq, flag)
    seq′ = remove_newlines(seq)
    if flag == "dna" || flag == "d"
        return DNAMer(seq′)
    elseif flag == "rna" || flag == "r"
        return RNAMer(seq′)
    else
        error("Invalid type flag: '$(flag)'")
    end
end

macro mer_str(seq)
    return DNAMer(remove_newlines(seq))
end

macro bigmer_str(seq, flag)
    seq′ = remove_newlines(seq)
    if flag == "dna" || flag == "d"
        return BigDNAMer(seq′)
    elseif flag == "rna" || flag == "r"
        return BigRNAMer(seq′)
    else
        error("Invalid type flag: '$(flag)'")
    end
end

macro bigmer_str(seq)
    return BigDNAMer(remove_newlines(seq))
end