text:
using Revise
using ADCME
using NNFEM
using JLD2
using PyPlot
using MAT
using DelimitedFiles
using Interpolations
reset_default_graph()

include("nnutil.jl")

loss = constant(0.0)
for tid = [1,2,4,5]
    strain, stress = read_strain_stress("Data/$(tid).dat")
    X, Y = prepare_strain_stress_data1D(strain, stress )
    y = squeeze(nn(constant(X[:,1:1]), constant(X[:,2:2]), constant(X[:,3:3])))
    global loss += sum((y - Y)^2)
end
sess = Session(); init(sess)
@show run(sess,loss)

loss_all = []
for i = 1:10
    println("************************** Outer Iteration = $i ************************** ")
    loss_ = BFGS!(sess, loss, 3000)
    push!(loss_all, loss_)
end

using Random;
file = randstring(8)
ADCME.save(sess, "$file.mat")
nnname = "$file.mat"
s = ae_to_code(nnname, "default", activation = activation)
eval(Meta.parse(s))
rm("$file.mat")

testtype = "NeuralNetwork1D"
include("NNTrussPull_Domain.jl")
tid = 3

nodes, EBC, g, gt, FBC, fext, ft, npoints, node_to_point = BoundaryCondition(tid)
prop = Dict("name"=> testtype, "rho"=> 8000.0, "E"=> 200e3, "nu"=> 0.45,
           "sigmaY"=>0.3e3, "K"=>1/9*200e3, "B"=> 0.0, "A0"=> 0.005, "nn"=>post_nn)

elements = []
for i = 1:nx
    elnodes = [i, i+1]; coords = nodes[elnodes,:];
    push!(elements, FiniteStrainTruss(coords,elnodes, prop, ngp))
end
domain = Domain(nodes, elements, ndofs, EBC, g, FBC, fext)
state = zeros(domain.neqs)
∂u = zeros(domain.neqs)
globdat = GlobalData(state,zeros(domain.neqs), zeros(domain.neqs),∂u, domain.neqs, gt, ft)
assembleMassMatrix!(globdat, domain)
# need to update state in domain from globdat
updateStates!(domain, globdat)

adaptive_solver_args = Dict("Newmark_rho"=> 0.0,
                          "Newton_maxiter"=>10,
                          "Newton_Abs_Err"=>1e-4,
                          "Newton_Rel_Err"=>1e-6,
                          "damped_Newton_eta" => 1.0)
globdat, domain, ts = AdaptiveSolver("NewmarkSolver", globdat, domain, T, NT, adaptive_solver_args)

state_ = hcat(domain.history["state"]...)
state2 = zeros(size(state_,1), NT+1)
for i = 1:size(state_,1)
    state2[i,:] = LinearInterpolation(ts, state_[i,:])(LinRange(0,T,NT+1))
end

@load "Data/domain3.jld2" domain
state = hcat(domain.history["state"]...)
loss = sqrt(mean((state2 - state).^2))
open("Data/$FILEID2.txt", "a") do io
    writedlm(io, [loss])
end


text:
# -----------------------------------------------------
# This script gives the methods for maniupulating the
# Cartesian indices corresponding to multilevels that
# are either kept or cut in various schemes of
# "sparsification"
# -----------------------------------------------------

# Efficiency Criticality: LOW
# This is not evaluated very often

# Accuracy Criticality: N/A
# No float manipulation

# Gives the appropriate boolean cutoff corresponding
# to a given scheme, e.g. sparse basis, full basis,
# and in the future possibly the energy basis of
# Bungartz and Griebel

function get_cutoff(scheme::String, D::Int, n::Int)
	if scheme == "sparse"
		return x -> (sum(x.I) > n+D)
	elseif scheme == "full"
		return x -> false
	else
		throw(ArgumentError)
	end
end

ext: jl


text:
include("linear_elasticity.jl")

E = Variable(1.0)
ν = Variable(0.0)
H11 = E*(1. -ν)/((1+ν)*(1. -2. *ν));
H12 = H11*ν/(1-ν);
H21 = H12;
H22 = H11;
H33 = H11*0.5*(1. -2. *ν)/(1. -ν);
H = tensor(
  [H11 H12 0.0
  H21 H22 0.0
  0.0 0.0 H33]
)
d, v, a= ExplicitSolver(globaldata, domain, d0, v0, a0, Δt, NT, H, Fext, ubd, abd)
idx = [1;domain.nnodes+1]
loss = sum((d[:, idx] - d_[:, idx])^2)
sess = Session(); init(sess)
BFGS!(sess, loss, var_to_bounds=Dict(E=>(0.0,100.0), ν=>(-0.5,0.49)))
run(sess, [E, ν])


text:
@inline function check_index(dict::CardinalPairDict{V,K}, idx1::K, idx2::K) where {K<:Signed, V}
    return (one(K) <= idx1 <= dict.idx1max) && (one(K) <= idx2 <= dict.idx2max)
end
@inline function check_index(dict::CardinalPairDict{V,K}, idx1::J, idx2::J) where {J<:Signed,K,V}
    return check_index(dict, idx1%K, idx2%K)
end

@inline function guarded_index(dict::CardinalPairDict{V,K}, idx1::K, idx2::K) where {K<:Signed, V}
    check_index(dict, idx1, idx2) || throw(ErrorException("invalid index ($(idx1), $(idx2))"))
    return nothing
end
@inline function guarded_index(dict::CardinalPairDict{V,K}, idx1::J, idx2::J) where {J<:Signed,K,V}
    return guarded_index(dict, idx1%K, idx2%K)
end

@inline function check_present(dict::CardinalPairDict{V,K}, idx1::K, idx2::K) where {K<:Signed, V}
    return getindex(dict.present, pair(idx1, idx2))
end
@inline function check_present(dict::CardinalPairDict{V,K}, idx1::J, idx2::J) where {J<:Signed,K,V}
    return check_present(dict, idx1%K, idx2%K)
end
@inline function check_present(dict::CardinalPairDict{V,K}, paired::K) where {K<:Signed, V}
    return getindex(dict.present, paired)
end

@inline function guarded_present(dict::CardinalPairDict{V,K}, idx1::K, idx2::K) where {K<:Signed, V}
    check_present(dict, idx1, idx2) || throw(ErrorException("there is no content at index ($(idx1), $(idx2))"))
    return nothing
end
@inline function guarded_present(dict::CardinalPairDict{V,K}, idx1::J, idx2::J) where {J<:Signed,K,V}
    return guarded_present(dict, idx1%K, idx2%K)
end

@inline function guarded_present(dict::CardinalPairDict{V,K}, paired::K) where {K<:Signed, V}
    if !check_present(dict, paired)
        idx1, idx2 = unpair(paired)
        throw(ErrorException("there is no content at index ($(idx1), $(idx2))"))
    end
    return nothing
end

@inline function keyval(dict::CardinalPairDict{V,K}, idx1::K, idx2::K) where {K<:Signed, V}
    guarded_index(dict, idx1, idx2)
    return pair(idx1, idx2)
end
@inline function keyval(dict::CardinalPairDict{V,K}, idx1::J, idx2::J) where {J<:Signed,K,V}
    return keyval(dict, idx1%K, idx2%K)
end

@inline function unsafe_keyval(dict::CardinalPairDict{V,K}, idx1::K, idx2::K) where {K<:Signed, V}
    return pair(idx1, idx2)
end
@inline function unsafe_keyval(dict::CardinalPairDict{V,K}, idx1::J, idx2::J) where {J<:Signed,K,V}
    return idx1%K, idx2%K
end

# some essential Julia associative, dict API functions

function Base.haskey(dict::CardinalPairDict{V,K}, idx1::K, idx2::K) where {K<:Signed, V}
    guarded_index(dict, idx1, idx2)
    return check_present(dict, idx1, idx2)
end
@inline function Base.haskey(dict::CardinalPairDict{V,K}, idx1::J, idx2::J) where {J<:Signed, K<:Signed, V}
    return haskey(dict, idx1%K, idx2%K)
end
@inline function Base.haskey(dict::CardinalPairDict{V,K}, idx12::Tuple{K,K}) where {K<:Signed, V}
    return haskey(dict, idx12[1], idx12[2])
end
@inline function Base.haskey(dict::CardinalPairDict{V,K}, idx12::Tuple{J,J}) where {J<:Signed, K<:Signed, V}
    return haskey(dict, idx12[1]%K, idx12[2]%K)
end

# this is quite unsafe
function unsafe_haskey(dict::CardinalPairDict{V,K}, idx1::K, idx2::K) where {K<:Signed, V}
    return check_present(dict, idx1, idx2)
end
@inline function unsafe_haskey(dict::CardinalPairDict{V,K}, idx1::J, idx2::J) where {J<:Signed, K<:Signed, V}
    return unsafe_haskey(dict, idx1%K, idx2%K)
end
function unsafe_haskey(dict::CardinalPairDict{V,K}, paired::K) where {K<:Signed, V}
    return check_present(dict, paired)
end

function Base.getindex(dict::CardinalPairDict{V,K}, idx1::K, idx2::K) where {K<:Signed, V}
    guarded_index(dict, idx1, idx2)
    paired = pair(idx1, idx2)
    guarded_present(dict, paired)
    return getindex(dict.content, paired)
end
@inline function Base.getindex(dict::CardinalPairDict{V,K}, idx1::J, idx2::J) where {J<:Signed, K<:Signed, V}
    return getindex(dict, idx1%K, idx2%K)
end
@inline function Base.getindex(dict::CardinalPairDict{V,K}, idx12::Tuple{K,K}) where {K<:Signed, V}
    return getindex(dict, idx12[1], idx12[2])
end
@inline function Base.getindex(dict::CardinalPairDict{V,K}, idx12::Tuple{J,J}) where {J<:Signed, K<:Signed, V}
    return getindex(dict, idx12[1]%K, idx12[2]%K)
end

function unsafe_getindex(dict::CardinalPairDict{V,K}, idx1::K, idx2::K) where {K<:Signed, V}
    paired = pair(idx1, idx2)
    guarded_present(dict, idx1, idx2)
    return getindex(dict.content, idx1, idx2)
end
@inline function unsafe_getindex(dict::CardinalPairDict{V,K}, idx1::J, idx2::J) where {J<:Signed, K<:Signed, V}
    return unsafe_getindex(dict, idx1%K, idx2%K)
end

function Base.get(dict::CardinalPairDict{V,K}, idx1::K, idx2::K, default::V) where {K<:Signed, V}
    guarded_index(dict, idx1, idx2)
    return check_present(dict, idx1, idx2) ? getindex(dict.content, idx1, idx2) : default
end
@inline function Base.get(dict::CardinalPairDict{V,K}, idx1::J, idx2::J, default::V) where {J<:Signed, K<:Signed, V}
    return get(dict, idx1%K, idx2%K, default)
end
@inline function Base.get(dict::CardinalPairDict{V,K}, idx12::Tuple{K,K}, default::V) where {K<:Signed, V}
    return get(dict, idx12[1], idx12[2], default)
end
@inline function Base.get(dict::CardinalPairDict{V,K}, idx12::Tuple{J,J}, default::V) where {J<:Signed, K<:Signed, V}
    return get(dict, idx12[1]%K, idx12[2]%K, default)
end

function Base.get!(dict::CardinalPairDict{V,K}, idx1::K, idx2::K, default::V) where {K<:Signed, V}
    guarded_index(dict, idx1, idx2)
    result = default
    if check_present(dict, idx1, idx2)
       result = getindex(dict, idx1, idx2)
    else
       setindex!(dict, default, idx1, idx2)
    end
    return result
end
@inline function Base.get!(dict::CardinalPairDict{V,K}, idx1::J, idx2::J, default::V) where {J<:Signed, K<:Signed, V}
    return get!(dict, idx%K, idx2%K, default)
end
@inline function Base.get!(dict::CardinalPairDict{V,K}, idx12::Tuple{K,K}, default::V) where {K<:Signed, V}
    return get!(dict, idx12[1], idx12[2], default)
end
@inline function Base.get!(dict::CardinalPairDict{V,K}, idx12::Tuple{J,J}, default::V) where {J<:Signed, K<:Signed, V}
    return get!(dict, idx12[1]%K, idx12[2]%K, default)
end

function Base.setindex!(dict::CardinalPairDict{V,K}, value::V, idx1::K, idx2::K) where {K<:Signed, V}
     guarded_index(dict, idx1, idx2)
     paired = pair(idx1, idx2)
     setindex!(dict.present, true, paired)
     setindex!(dict.content, value, paired)
     return value
end
@inline function Base.setindex!(dict::CardinalPairDict{V,K}, value::V, idx1::J, idx2::J) where {J<:Signed, K<:Signed, V}
    return setindex!(dict, value, idx1%K, idx2%K)
end
@inline function Base.setindex!(dict::CardinalPairDict{V,K}, value::V, idx12::Tuple{K,K}) where {K<:Signed, V}
    return setindex!(dict, value, idx12[1], idx12[2])
end
@inline function Base.setindex!(dict::CardinalPairDict{V,K}, value::V, idx12::Tuple{J,J}) where {J<:Signed, K<:Signed, V}
    return setindex!(dict, value, idx12[1]%K, idx12[2]%K)
end

function unsafe_setindex!(dict::CardinalPairDict{V,K}, value::V, idx1::K, idx2::K) where {K<:Signed, V}
     paired = pair(idx1, idx2)
     setindex!(dict.present, true, paired)
     setindex!(dict.content, value, paired)
     return value
end
@inline function unsafe_setindex!(dict::CardinalPairDict{V,K}, value::V, idx1::J, idx2::J) where {J<:Signed, K<:Signed, V}
    return unsafe_setindex!(dict, value, idx1%K, idx2%K)
end


text:
using TextSegmentation
using Test
using PyCall

document = [
    # Derek J. de Solla Price. (2009). On the Origin of Clockwork, Perpetual Motion Devices, and the Compass. Urbana, Illinois: Project Gutenberg. Retrieved February 25, 2022, from www.gutenberg.org/ebooks/30001.
    "Once the floodgates of Arabic learning were opened, a stream of mechanized astronomical models poured into Europe.",
    "Astrolabes and equatoria rapidly became very popular, mainly through the reason for which they had been first devised, the avoidance of tedious written computation.",
    "Many medieval astrolabes have survived, and at least three medieval equatoria are known.",
    "Chaucer is well known for his treatise on the astrolabe; a manuscript in Cambridge, containing a companion treatise on the equatorium, has been tentatively suggested by the present author as also being the work of Chaucer and the only piece written in his own hand.",
    # Carroll, Lewis. (2006). Alice in Wonderland. Urbana, Illinois: Project Gutenberg. Retrieved February 25, 2022, from www.gutenberg.org/ebooks/19033.
    "\"You are not attending!\" said the Mouse to Alice, severely.",
    "\"What are you thinking of?\"",
    "\"I beg your pardon,\" said Alice very humbly, \"you had got to the fifth bend, I think?\"",
    "\"You insult me by talking such nonsense!\" said the Mouse, getting up and walking away.",
    "\"Please come back and finish your story!\" Alice called after it.",
    # John Stuart Mill. (2011). On Liberty. Urbana, Illinois: Project Gutenberg. Retrieved February 25, 2022, from www.gutenberg.org/ebooks/34901.
    "The struggle between Liberty and Authority is the most conspicuous feature in the portions of history with which we are earliest familiar, particularly in that of Greece, Rome, and England.",
    "But in old times this contest was between subjects, or some classes of subjects, and the government.",
    "By liberty, was meant protection against the tyranny of the political rulers.",
    "The rulers were conceived (except in some of the popular governments of Greece) as in a necessarily antagonistic position to the people whom they ruled.",
    ]

@testset "TextTiling.jl" begin
    window_size = 2
    smooth_window_size = 1
    num_topics = 3
    tt = TextTiling.SegmentObject(window_size, smooth_window_size, Utls.tokenize)
    result = TextTiling.segment(tt, document, num_topics)
    @test result == "000100010000"
end

@testset "C99.jl" begin
    n = length(document)
    init_matrix = zeros(n, n)
    window_size = 2
    std_coeff = 1.2
    c99 = C99.SegmentObject(window_size, init_matrix, init_matrix, init_matrix, std_coeff, Utls.tokenize)
    result = C99.segment(c99, document, n)
    @test result == "000100001000"
end

@testset "TopicTiling.jl" begin
    # LDA Topic Model
    pygensim = pyimport("gensim")
    tokenized_train_document = [Utls.tokenize(i) for i in document]
    dictionary = pygensim.corpora.Dictionary(tokenized_train_document)
    corpus = [dictionary.doc2bow(text) for text in tokenized_train_document]
    lda_model = pygensim.models.ldamodel.LdaModel(
        corpus = corpus,
        id2word = dictionary,
        minimum_probability = 0.0001,
        num_topics = 3,
        random_state=1234,
    )

    # TopicTiling
    window_size = 2
    smooth_window_size = 1
    num_topics = 3
    to = TopicTiling.SegmentObject(window_size, smooth_window_size, lda_model, dictionary)
    result = TopicTiling.segment(to, document, num_topics)
    @test result == "000110000000"
end


text:
module UnitCellModelling

using gmsh
using Parameters
using StaticArrays
# using

include("ucmBase.jl")
include("gshapes.jl")
include("geometric_modelling.jl")
include("meshing.jl")

"""
    make_unit_cell_model() -> Dict

It performs mainly three tasks
 1. Geometry creation
 2. Mesh generation
 3. Exporting geometry and mesh data

Returns a julia dictionary with the following key-value pairs

- `"mesh_data"` => `Dict{String, Any}` , where possible key-value pairs of dictionary are
    - `"all_node_tags"` => `Vector{Int}`
    - `"all_node_coordinates"` => `Matrix{Float64}`
    - `"matrix_element_connectivity"` => `Matrix{Int}`
    - `"inclusions_element_connectivity"` => `Matrix{Int}`
    - `"total_nodes"` => Int
    - `"num_2D_elements"` => Int
    - `"num_3D_elements"` => Int

"""
function make_unit_cell_model(
    unit_cell::AbstractUnitCell;
    small_parameter::Float64 = 1e-06,
    geom_export_paths::Tuple{Vararg{String}} = (),  # joinpath(homedir(), "unit_cell.step"),
    extr_dir::String = "XY->Z",
    extr_dir_num_ele::Vector{Int64} = Int64[],
    extr_dir_cum_heights::Vector{Float64} = Float64[],
    extr_dir_recombine_ele::Bool = true,
    mesh_periodicity::Bool = true,
    element_types::Tuple{Vararg{Symbol}} = (:DEFAULT,),
    min_ele_size_factor::Float64 = 1.0,   # FIXME make default to take factor, instead of absolute values
    max_ele_size_factor::Float64 = 2.0,
    mesh_opt_algorithm::String = "Netgen",
    show_mesh_stats::Bool = true,
    show_rve::Bool = true,
)::Dict{String,Any}
    unit_cell_dim::Int = dimension(unit_cell)
    small_uc_side_length = minimum(side_lengths(unit_cell))
    global ϵ = small_parameter

    #
    gmsh.initialize()
    gmsh.option.set_number("General.Verbosity", 0)
    gmsh.option.set_number("General.Terminal", 0)
    gmsh.option.set_number("General.NoPopup", 0)
    gmsh.option.set_number("Geometry.OCCBoundsUseStl", 1)
    gmsh.model.add("unit_cell")
    #
    # ---------------------
    #   Geometry creation
    # ---------------------
    ucp_geom_tags = geometric_model(
        unit_cell,
        extr_dir,
        extr_dir_num_ele,
        extr_dir_cum_heights,
        extr_dir_recombine_ele
    )  # unit cell phase tags
    matrix_pg_tag = gmsh.model.add_physical_group(unit_cell_dim, ucp_geom_tags[1:1])
    inclusions_pg_tags = gmsh.model.add_physical_group(unit_cell_dim, ucp_geom_tags[2:end])
    gmsh.model.set_physical_name(unit_cell_dim, matrix_pg_tag, "Matrix")
    gmsh.model.set_physical_name(unit_cell_dim, inclusions_pg_tags, "Inclusions")
    #
    # ---------------------
    #       Meshing
    # ---------------------
    generate_mesh(
        unit_cell,
        mesh_periodicity,
        element_types,
        extr_dir_num_ele,
        min_ele_size_factor*small_uc_side_length,
        max_ele_size_factor*small_uc_side_length,
        mesh_opt_algorithm,
        show_mesh_stats
    )
    #
    # ---------------------
    #       EXPort
    # ---------------------

    # exporting UC model to the specified fomats
    if !isempty(geom_export_paths)
        for a_export_path in geom_export_paths
            gmsh.write(a_export_path)
        end
    end

    #
    model_data::Dict{String,Any} = Dict{String,Any}()
    model_data["mesh_data"] = get_mesh_data(unit_cell_dim, ucp_geom_tags)  # Dict{String, Any}
    #
    if show_rve
        gmsh.fltk.run()
    end
    #
    gmsh.finalize()
    #
    return model_data
end  # of function

export make_unit_cell_model

end  # of UnitCellModelling module


text:
include("firfilt.jl")
export  FIRFilter

include("firinterp.jl")
export  FIRInterp

include("firdecim.jl")
export  FIRDecim,
        gettaps

# Export common methods
export  destroy,
        print,
        reset!,
        execute,
        push!,
        freqresponse,
        groupdelay

max_stars_repo_name: JayKickliter/LiquidDSP

max_issues_repo_name: JayKickliter/LiquidDSP

max_forks_repo_name: JayKickliter/LiquidDSP


text:
println("Hi from external.jl")

function myfunc(x)
    return sqrt(x)
end

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
"""
    lsnormal(A,b)

Solve a linear least squares problem by the normal equations.
Returns the minimizer of ||b-Ax||.
"""
function lsnormal(A,b)
    N = A'*A;  z = A'*b;
    R = cholesky(N).U
    w = forwardsub(R',z)                   # solve R'z=c
    x = backsub(R,w)                       # solve Rx=z
    return x
end

"""
    lsqrfact(A,b)

Solve a linear least squares problem by QR factorization. Returns
the minimizer of ||b-Ax||.
"""
function lsqrfact(A,b)
    Q,R = qr(A)
    c = Q'*b
    x = backsub(R,c)
    return x
end

"""
    qrfact(A)

QR factorization by Householder reflections. Returns Q and R.
"""
function qrfact(A)
    m,n = size(A)
    Qt = diagm(ones(m))
    R = float(copy(A))
    for k in 1:n
        z = R[k:m,k]
        w = [ -sign(z[1])*norm(z) - z[1]; -z[2:end] ]
        nrmw = norm(w)
        if nrmw < eps() continue; end    # skip this iteration
        v = w / nrmw;
        # Apply the reflection to each relevant column of R and Q
        for j in k:n
            R[k:m,j] -= v*( 2*(v'*R[k:m,j]) )
        end
        for j in 1:m
            Qt[k:m,j] -= v*( 2*(v'*Qt[k:m,j]) )
        end
    end
    return Qt',triu(R)
end


text:
module ExponentialUtilities
using LinearAlgebra, SparseArrays, Printf

"""
    @diagview(A,d) -> view of the `d`th diagonal of `A`.
"""
macro diagview(A,d::Integer=0)
    s = d<=0 ? 1+abs(d) : :(m+$d)
    quote
        m = size($(esc(A)),1)
        @view($(esc(A))[($s):m+1:end])
    end
end

include("exp.jl")
include("phi.jl")
include("arnoldi.jl")
include("krylov_phiv.jl")
include("krylov_phiv_adaptive.jl")
include("kiops.jl")
include("StegrWork.jl")
include("krylov_phiv_error_estimate.jl")

export phi, phi!, KrylovSubspace, arnoldi, arnoldi!, lanczos!, ExpvCache, PhivCache,
    expv, expv!, exp_generic, phiv, phiv!, kiops, expv_timestep, expv_timestep!, phiv_timestep, phiv_timestep!,
    StegrCache, get_subspace_cache
end


text:
macro goma_str(anything)
    println("キュー")
end

@show goma"ゴマちゃんの鳴き声は"
@show goma"9="
@show goma"3x3="


text:
using LinearAlgebra
using SparseArrays
using Random

Random.seed!(0)

# Test with full matrices

n, k = 5, 4

A = randn(n, k)
Q = A*A'
q = randn(n)
f = Quadratic(Q, q)
@test typeof(f) <: ProximalOperators.QuadraticDirect

predicates_test(f)

@test ProximalOperators.is_smooth(f) == true
@test ProximalOperators.is_quadratic(f) == true
@test ProximalOperators.is_set(f) == false

x = randn(n)

grad_fx, fx = gradient(f, x)
@test abs(fx - 0.5*dot(x, Q*x) - dot(x, q)) <= 1e-12
@test norm(grad_fx - (Q*x + q), Inf) <= 1e-12

call_test(f, x)
prox_test(f, x)
prox_test(f, x, 1.5)

f = Quadratic(Q, q, iterative=true)
@test typeof(f) <: ProximalOperators.QuadraticIterative

call_test(f, x)
prox_test(f, x)
prox_test(f, x, 2.1)

# Test with sparse matrices

Q = sparse(Q)
f = Quadratic(Q, q)
@test typeof(f) <: ProximalOperators.QuadraticDirect

grad_fx, fx = gradient(f, x)
@test abs(fx - 0.5*dot(x, Q*x) - dot(x, q)) <= 1e-12
@test norm(grad_fx - (Q*x + q), Inf) <= 1e-12

call_test(f, x)
prox_test(f, x)
prox_test(f, x, 0.8)

f = Quadratic(Q, q, iterative=true)
@test typeof(f) <: ProximalOperators.QuadraticIterative

call_test(f, x)
prox_test(f, x)
prox_test(f, x, 1.3)


text:
module Libav

import BinDeps
using Dates
using SHA
using Random
using Libaudio

folder() = normpath(joinpath(@__FILE__, "../../deps"))

function uninstall()
    rm(joinpath(folder(), "FFmpeg"), force=true, recursive=true)
end

function isinstalled()
    ok = isdir(joinpath(folder(), "FFmpeg"))
end

download(x) = run(BinDeps.download_cmd(x, basename(x)))

function fetch(url, file, sha256r, app)
    download("$(url)$(file)")
    sha256c = open("$file") do io
        sha256(io)
    end
    if isequal(bytes2hex(sha256c), lowercase(sha256r))
        printstyled("SHA256 ok, extracting...\n", color=:light_cyan)
        run(`7z x $file -o$app`)
    else
        printstyled("SHA256 mismatch, nothing installed, please try again\n", color=:light_red)
    end
    rm(file)
end

function install()
    dir = folder()
    !isdir(dir) && mkpath(dir)
    uninstall()
    cd(dir) do
      if Sys.iswindows()
        # arch = Int == Int64 ? "x64" : "x86"
        url = "https://ffmpeg.zeranoe.com/builds/win64/static/"
        file = "ffmpeg-4.0.2-win64-static.zip"
        sha256r = "2BF3726D7B489F1FE5609E1F9D93D827FD9E9C9D21B829B254A5BAE0C3F60EC8"
        fetch(url, file, sha256r, "FFmpeg")
      end
    end
end

function listdevices()
    ffmpeg = joinpath(folder(), "FFmpeg/ffmpeg-4.0.2-win64-static/bin/ffmpeg.exe")
    try
        Sys.iswindows() && run(`$ffmpeg -list_devices true -f dshow -i dummy`)
        Sys.islinux() && run(`v4l2-ctl --list-devices`)
    catch
    end
end

function listoptions(device="Logitech HD Webcam C310")
    ffmpeg = joinpath(folder(), "FFmpeg/ffmpeg-4.0.2-win64-static/bin/ffmpeg.exe")
    try
        Sys.iswindows() && run(`$ffmpeg -f dshow -list_options true -i video="$(device)"`)
        Sys.islinux() && run(`ffmpeg -f v4l2 -list_formats all -i /dev/video0`)
    catch
    end
end

function record(t, mp4="foobar.mp4", s="160x120", f=30, v="Logitech HD Webcam C310", a="Microphone (HD Webcam C310)")
    ffmpeg = joinpath(folder(), "FFmpeg/ffmpeg-4.0.2-win64-static/bin/ffmpeg.exe")
    try
        Sys.iswindows() && run(`$ffmpeg -y -f dshow -t $t -video_size $s -framerate $f -pixel_format bgr24 -i video="$v":audio="$a" $mp4`)
        if Sys.islinux()
            wav = mp4[1:end-3] * "wav"
            open("capture.sh", "w") do io
                write(io, "ffmpeg -y -t $t -f alsa -ac 1 -ar 48000 -i hw:1 $wav &\n")
                write(io, "ffmpeg -y -t $t -f v4l2 -framerate $f -video_size $s -i /dev/video0 $mp4\n")
            end
            run(`chmod +x capture.sh`)
            run(`./capture.sh`)
        end
    catch
    end
end

function ripaudio(mp4="foobar.mp4", wav="foobar.wav")
    ffmpeg = joinpath(folder(), "FFmpeg/ffmpeg-4.0.2-win64-static/bin/ffmpeg.exe")
    try
        run(`$ffmpeg -y -i $mp4 $wav`)
    catch
    end
end

function ripvideo(mp4="foobar.mp4", ppm="foobar.ppm")
    ffmpeg = joinpath(folder(), "FFmpeg/ffmpeg-4.0.2-win64-static/bin/ffmpeg.exe")
    try
        run(`$ffmpeg -y -i $mp4 -f image2pipe -vcodec ppm $ppm`)
    catch
    end
end

#P6
#1920 1080
#255
#<frame RGB>
function ppmheader(ppm)
    header = zeros(UInt8, 128)
    open(ppm, "r") do io
        read!(io, header)
    end
    return header
end

function udiff(a::UInt8, b::UInt8)
    a > b ? (a-b) : (b-a)
end

function rgbintense(ppm="foobar.ppm", h=480, v=640, fps=30, tsilent::Tuple{<:Real,<:Real}=(1,2))

    p6 = [0x50, 0x36, 0x0a]
    depth = [0x32, 0x35, 0x35, 0x0a]
    out = randstring()
    noisefloor = out * ".ppm"

    header = zeros(UInt8, length(p6) + length(string(h)) + 1 + length(string(v)) + 1 + length(depth))
    f0 = zeros(UInt8, h*v*3)
    fn = zeros(UInt64, h*v*3)
    fx = zeros(UInt8, h*v*3)

    n = 1
    open(ppm, "r") do io
        while !eof(io)
            read!(io, header)
            read!(io, fx)
            n += 1
        end
    end
    rgb = zeros(Float32, n, 3)

    open(ppm, "r") do io
        i1 = round(Int, tsilent[1] * fps)
        i2 = round(Int, tsilent[2] * fps)
        for i = 1:i1
            read!(io, header)
            read!(io, f0)
        end
        for i = i1+1:i2
            read!(io, header)
            read!(io, f0)
            fn .+= f0
        end
        f0 .= convert(Vector{UInt8},div.(fn,i2-i1))
    end
    open(noisefloor, "w") do io
        write(io, header)
        write(io, f0)
    end

    open(ppm, "r") do io
        i = 1
        while !eof(io)
            read!(io, header)
            read!(io, fx)
            # opt 1) signal[i] = sum(frame) - md
            # opt 2) signal[i] = udiff.(sum(frame), md)
            for k = 1:3
                rgb[i,k] = sum(udiff.(fx[k:3:end], f0[k:3:end])) / (h*v)
            end
            i += 1
        end
        @assert i == n
    end
    # open(out * ".dat", "w") do f
    #     for j = 1:n
    #         write(f, "$(rgb[j,1]), $(rgb[j,2]), $(rgb[j,3])\n")
    #     end
    # end
    return rgb, out
end

function wavrgbintense(t, mp4, tsilent, wav)

    r = randstring()
    d = r * ".ppm"
    w = r * ".wav"

    Libav.record(t, mp4)
    Libav.ripvideo(mp4, d)
    rgb, out = Libav.rgbintense(d, 160, 120, 30, tsilent)
    Sys.iswindows() && Libav.ripaudio(mp4, w)
    x,fs = Libaudio.wavread(w)

    # norm rgb to [0,1]
    # mark x axis of rgb to sound track
    # merge new tracks

    rgb ./= maximum(rgb)
    nx = size(x,1)
    rgbw = zeros(Float32, nx, 3)
    n = size(rgb,1)
    for i = 1:n
        rgbw[max(min(round(Int, i/30 * convert(Float64, fs)), nx),1), :] = rgb[i,:]
    end
    y = [x rgbw]
    Libaudio.wavwrite(wav, y, fs, 16)
    rm(d)
    rm(w)
    return y
end

end # module


text:
using Dates: days, hour, minute, second, millisecond

function Dates.DateTime(lzdt::LaxZonedDateTime)
    if isrepresentable(lzdt)
        return lzdt.local_datetime
    else
        error("Unable to determine local datetime from an unrepresentable LaxZonedDateTime")
    end
end

Dates.DateTime(lzdt::LaxZonedDateTime, ::Type{UTC}) = DateTime(ZonedDateTime(lzdt), UTC)

Dates.Date(lzdt::LaxZonedDateTime) = Date(DateTime(lzdt))
Dates.Date(lzdt::LaxZonedDateTime, ::Type{UTC}) = Date(DateTime(lzdt, UTC))

Dates.Time(lzdt::LaxZonedDateTime) = Time(DateTime(lzdt))
Dates.Time(lzdt::LaxZonedDateTime, ::Type{UTC}) = Time(DateTime(lzdt, UTC))

TimeZones.timezone(lzdt::LaxZonedDateTime) = lzdt.timezone

"""
    isrepresentable(lzdt::LaxZonedDateTime) -> Bool

Indicates whether a `LaxZonedDateTime` can be represented by a local `DateTime`. Both valid
and invalid time are representable where as unknown time is not representable.
"""
isrepresentable(lzdt::LaxZonedDateTime) = lzdt.representable

"""
    isvalid(lzdt::LaxZonedDateTime) -> Bool

Indicates if a `LaxZonedDateTime` is a valid time. Valid time can always be represented as a
UTC instant.
"""
Base.isvalid(lzdt::LaxZonedDateTime) = isrepresentable(lzdt) && !isa(lzdt.zone, InvalidTimeZone)

"""
    isinvalid(lzdt::LaxZonedDateTime) -> Bool

Indicates if a `LaxZonedDateTime` is an invalid time. Invalid time is representable in local
time but is not representable as a UTC instant.
"""
isinvalid(lzdt::LaxZonedDateTime) = isrepresentable(lzdt) && isa(lzdt.zone, InvalidTimeZone)

isambiguous(lzdt::LaxZonedDateTime) = isa(lzdt.zone, Ambiguous)
isnonexistent(lzdt::LaxZonedDateTime) = isa(lzdt.zone, NonExistent)

Dates.days(lzdt::LaxZonedDateTime) = days(DateTime(lzdt))

for period in (:Hour, :Minute, :Second, :Millisecond)
    accessor = Symbol(lowercase(string(period)))
    @eval begin
        Dates.$accessor(lzdt::LaxZonedDateTime) = $accessor(DateTime(lzdt))
        Dates.$period(lzdt::LaxZonedDateTime) = $period($accessor(lzdt))
    end
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct ServiceAssociationLinkPropertiesFormat <: SwaggerModel
    linkedResourceType::Any # spec type: Union{ Nothing, String } # spec name: linkedResourceType
    link::Any # spec type: Union{ Nothing, String } # spec name: link
    provisioningState::Any # spec type: Union{ Nothing, ProvisioningState } # spec name: provisioningState
    allowDelete::Any # spec type: Union{ Nothing, Bool } # spec name: allowDelete
    locations::Any # spec type: Union{ Nothing, Vector{String} } # spec name: locations

    function ServiceAssociationLinkPropertiesFormat(;linkedResourceType=nothing, link=nothing, provisioningState=nothing, allowDelete=nothing, locations=nothing)
        o = new()
        validate_property(ServiceAssociationLinkPropertiesFormat, Symbol("linkedResourceType"), linkedResourceType)
        setfield!(o, Symbol("linkedResourceType"), linkedResourceType)
        validate_property(ServiceAssociationLinkPropertiesFormat, Symbol("link"), link)
        setfield!(o, Symbol("link"), link)
        validate_property(ServiceAssociationLinkPropertiesFormat, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        validate_property(ServiceAssociationLinkPropertiesFormat, Symbol("allowDelete"), allowDelete)
        setfield!(o, Symbol("allowDelete"), allowDelete)
        validate_property(ServiceAssociationLinkPropertiesFormat, Symbol("locations"), locations)
        setfield!(o, Symbol("locations"), locations)
        o
    end
end # type ServiceAssociationLinkPropertiesFormat

const _property_map_ServiceAssociationLinkPropertiesFormat = Dict{Symbol,Symbol}(Symbol("linkedResourceType")=>Symbol("linkedResourceType"), Symbol("link")=>Symbol("link"), Symbol("provisioningState")=>Symbol("provisioningState"), Symbol("allowDelete")=>Symbol("allowDelete"), Symbol("locations")=>Symbol("locations"))
const _property_types_ServiceAssociationLinkPropertiesFormat = Dict{Symbol,String}(Symbol("linkedResourceType")=>"String", Symbol("link")=>"String", Symbol("provisioningState")=>"ProvisioningState", Symbol("allowDelete")=>"Bool", Symbol("locations")=>"Vector{String}")
Base.propertynames(::Type{ ServiceAssociationLinkPropertiesFormat }) = collect(keys(_property_map_ServiceAssociationLinkPropertiesFormat))
Swagger.property_type(::Type{ ServiceAssociationLinkPropertiesFormat }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ServiceAssociationLinkPropertiesFormat[name]))}
Swagger.field_name(::Type{ ServiceAssociationLinkPropertiesFormat }, property_name::Symbol) =  _property_map_ServiceAssociationLinkPropertiesFormat[property_name]

function check_required(o::ServiceAssociationLinkPropertiesFormat)
    true
end

function validate_property(::Type{ ServiceAssociationLinkPropertiesFormat }, name::Symbol, val)
end


text:
using CausalDiscovery.Autumn
using CausalDiscovery.Autumn.Transform: Statement
# using CausalDiscovery.Autumn.Transform: findnonterminal

function test1()
  p = AExpr(:program, Statement())
  ϕ = Phi()
  subexpr_ = subexpr(p, 1)
  sub(ϕ, subexpr_)
end

function test2()
  aex = AExpr(:program, Statement(), Statement(), Statement(), Statement())
  ϕ = Phi()
  recursub(ϕ, aex)
end

@testset "transform" begin
  test1()
  test2()
end


text:
using Test
using SixelTerm


text:
"""
    Sequential linear programming with line search
"""
mutable struct SlpLS{T,Tv,Tt} <: AbstractSlpOptimizer
    problem::Model{T,Tv,Tt} # problem data

    x::Tv # primal solution
    p::Tv # search direction
    p_slack::Dict{Int,Tv} # search direction at feasibility restoration phase
    lambda::Tv # Lagrangian dual multiplier
    mult_x_L::Tv # reduced cost for lower bound
    mult_x_U::Tv # reduced cost for upper bound

    # Evaluations at `x`
    f::T # objective function
    df::Tv # gradient
    E::Tv # constraint
    dE::Tv # Jacobian

    phi::T # merit function value
    ν::Tv # penalty parameters for the merit function
    directional_derivative::T # directional derivative
    alpha::T # stepsize

    prim_infeas::T # primal infeasibility at `x`
    dual_infeas::T # dual (approximate?) infeasibility
    compl::T # complementary slackness

    optimizer::Union{Nothing,AbstractSubOptimizer} # Subproblem optimizer

    options::Parameters

    feasibility_restoration::Bool # indicator for feasibility restoration
    iter::Int # iteration counter
    ret::Int # solution status
    start_time::Float64 # solution start time
    start_iter_time::Float64 # iteration start time

    function SlpLS(problem::Model{T,Tv,Tt}) where {T,Tv<:AbstractArray{T},Tt}
        slp = new{T,Tv,Tt}()
        slp.problem = problem
        slp.x = Tv(undef, problem.n)
        slp.p = zeros(problem.n)
        slp.p_slack = Dict()
        slp.lambda = zeros(problem.m)
        slp.mult_x_L = zeros(problem.n)
        slp.mult_x_U = zeros(problem.n)
        slp.df = Tv(undef, problem.n)
        slp.E = Tv(undef, problem.m)
        slp.dE = Tv(undef, length(problem.j_str))
        slp.phi = Inf
        slp.ν = Tv(undef, problem.m)

        slp.alpha = 1.0

        slp.prim_infeas = Inf
        slp.dual_infeas = Inf
        slp.compl = Inf

        slp.options = problem.parameters
        slp.optimizer = nothing

        slp.feasibility_restoration = false
        slp.iter = 1
        slp.ret = -5
        slp.start_time = 0.0
        slp.start_iter_time = 0.0

        return slp
    end
end

"""
    run!

Run the line-search SLP algorithm
"""
function run!(slp::SlpLS)

    slp.start_time = time()

    if slp.options.OutputFlag == 1
        sparsity_val = ifelse(
            slp.problem.m > 0,
            length(slp.problem.j_str) / (slp.problem.m * slp.problem.n),
            0.0,
        )
        @printf("LP subproblem sparsity: %e\n", sparsity_val)
        add_statistics(slp.problem, "sparsity", sparsity_val)
    else
        Logging.disable_logging(Logging.Info)
    end

    # Set initial point from MOI
    @assert length(slp.x) == length(slp.problem.x)
    slp.x .= slp.problem.x
    # Adjust the initial point to satisfy the column bounds
    for i = 1:slp.problem.n
        if slp.problem.x_L[i] > -Inf
            slp.x[i] = max(slp.x[i], slp.problem.x_L[i])
        end
        if slp.problem.x_U[i] > -Inf
            slp.x[i] = min(slp.x[i], slp.problem.x_U[i])
        end
    end

    # @show slp.problem.x_U
    # @show slp.x
    # @show slp.problem.x_L
    # if slp.options.ConvexFeasibleInitialization > 0
    #     slp.x = convex_optimize!(slp)
    # end

    if slp.options.ConvexFeasibleInitialization > 0
        slp.x = try
                    convex_initialization!(slp)
                 catch;
                    @info "Convex Feasible Initilization Failed"
                    slp.x
                 end
    end

    slp.iter = 1
    is_valid_step = true
    while true

        slp.start_iter_time = time()

        # evaluate function, constraints, gradient, Jacobian
        eval_functions!(slp)
        slp.alpha = 0.0
        slp.prim_infeas = norm_violations(slp, Inf)
        slp.dual_infeas = KT_residuals(slp)
        slp.compl = norm_complementarity(slp)

        LP_time_start = time()
        # solve LP subproblem (to initialize dual multipliers)
        slp.p, slp.lambda, slp.mult_x_U, slp.mult_x_L, slp.p_slack, status =
            sub_optimize!(slp)

        add_statistics(slp.problem, "LP_time", time() - LP_time_start)

        if status ∉ [MOI.OPTIMAL, MOI.INFEASIBLE]
            @warn("Unexpected LP subproblem solution status ($status)")
            slp.ret == -3
            if slp.prim_infeas <= slp.options.tol_infeas
                slp.ret = 6
            end
            break
        elseif status == MOI.INFEASIBLE
            if slp.feasibility_restoration == true
                @info "Failed to find a feasible direction"
                if slp.prim_infeas <= slp.options.tol_infeas
                    slp.ret = 6
                else
                    slp.ret = 2
                end
                break
            else
                # println("Feasibility restoration ($(status), |p| = $(norm(slp.p, Inf))) begins.")
                slp.feasibility_restoration = true
                continue
            end
        end

        compute_nu!(slp)
        slp.phi = compute_phi(slp)
        slp.directional_derivative = compute_derivative(slp)

        # step size computation
        is_valid_step = compute_alpha(slp)

        print_header(slp)
        print(slp)
        collect_statistics(slp)

        # Iteration counter limit
        if slp.iter >= slp.options.max_iter
            slp.ret = -1
            if slp.prim_infeas <= slp.options.tol_infeas
                slp.ret = 6
            end
            break
        end

        if (
            slp.prim_infeas <= slp.options.tol_infeas &&
            slp.compl <= slp.options.tol_residual
        ) || norm(slp.p, Inf) <= slp.options.tol_direction
            if slp.feasibility_restoration
                slp.feasibility_restoration = false
                slp.iter += 1
                continue
            elseif slp.dual_infeas <= slp.options.tol_residual
                slp.ret = 0
                break
            end
        end

        # Failed to find a step size
        if !is_valid_step
            @info "Failed to find a step size"
            if slp.ret == -3
                if slp.prim_infeas <= slp.options.tol_infeas
                    slp.ret = 6
                else
                    slp.ret = 2
                end
                break
            else
                slp.feasibility_restoration = true
            end

            slp.iter += 1
            continue
        end

        # update primal points
        slp.x += slp.alpha .* slp.p

        slp.iter += 1
    end
    slp.problem.obj_val = slp.problem.eval_f(slp.x)
    slp.problem.status = Int(slp.ret)
    slp.problem.x .= slp.x
    slp.problem.g .= slp.E
    slp.problem.mult_g .= slp.lambda
    slp.problem.mult_x_U .= slp.mult_x_U
    slp.problem.mult_x_L .= slp.mult_x_L
    add_statistic(slp.problem, "iter", slp.iter)
end

"""
    compute_alpha

Compute step size for line search
"""
function compute_alpha(slp::SlpLS)::Bool
    is_valid = true
    slp.alpha = 1.0
    phi_x_p = compute_phi(slp, slp.x, slp.alpha, slp.p)
    eta = slp.options.eta

    while phi_x_p > slp.phi + eta * slp.alpha * slp.directional_derivative
        # The step size can become too small.
        if slp.alpha < slp.options.min_alpha
            # @printf("Descent step cannot be computed.\n")
            # @printf("Feasibility restoration is required but not implemented yet.\n")
            if slp.feasibility_restoration
                slp.ret = -3
            end
            is_valid = false
            break
        end
        slp.alpha *= slp.options.tau
        phi_x_p = compute_phi(slp, slp.x, slp.alpha, slp.p)
    end
    # @show phi_x_p, slp.phi, slp.alpha, slp.directional_derivative, is_valid
    return is_valid
end

"""
    compute_nu!

Compute the penalty parameter for the merit function.
"""
function compute_nu!(slp::SlpLS)
    if slp.iter == 1
        for i = 1:slp.problem.m
            slp.ν[i] = abs(slp.lambda[i])
        end
    else
        for i = 1:slp.problem.m
            slp.ν[i] = max(slp.ν[i], abs(slp.lambda[i]))
        end
    end
end

"""
    compute_phi

Evaluate and return the merit function value at the current point x
"""
compute_phi(slp::SlpLS) = compute_phi(slp, slp.x, 0.0, slp.p)

"""
    print_header

Print the header of iteration information.
"""
function print_header(slp::SlpLS)
    if slp.options.OutputFlag == 0
        return
    end
    if (slp.iter - 1) % 25 == 0
        @printf("  %6s", "iter")
        @printf("  %15s", "f(x_k)")
        # @printf("  %15s", "ϕ(x_k)")
        @printf("  %15s", "D(ϕ,p)")
        # @printf("  %15s", "∇f^Tp")
        @printf("  %14s", "α")
        @printf("  %14s", "|p|")
        @printf("  %14s", "α|p|")
        # @printf("  %14s", "|∇f|")
        @printf("  %14s", "inf_pr")
        @printf("  %14s", "inf_du")
        @printf("  %14s", "compl")
        @printf("  %10s", "time")
        @printf("\n")
    end
end

"""
    print

Print iteration information.
"""
function print(slp::SlpLS)
    if slp.options.OutputFlag == 0
        return
    end
    st = ifelse(slp.feasibility_restoration, "FR", "  ")
    @printf("%2s%6d", st, slp.iter)
    @printf("  %+6.8e", slp.f)
    # @printf("  %+6.8e", slp.phi)
    @printf("  %+.8e", slp.directional_derivative)
    # @printf("  %+.8e", slp.df' * slp.p)
    @printf("  %6.8e", slp.alpha)
    @printf("  %6.8e", norm(slp.p, Inf))
    @printf("  %6.8e", slp.alpha * norm(slp.p, Inf))
    # @printf("  %.8e", norm(slp.df))
    @printf("  %6.8e", slp.prim_infeas)
    @printf("  %.8e", slp.dual_infeas)
    @printf("  %6.8e", slp.compl)
    @printf("  %10.2f", time() - slp.start_time)
    @printf("\n")
end

"""
    collect_statistics

Collect iteration information.
"""
function collect_statistics(slp::SlpLS)
    if slp.options.StatisticsFlag == 0
        return
    end
    mode = (slp.feasibility_restoration) ? "FR" : "OPT"
    add_statistics(slp.problem, "f(x)", slp.f)
    add_statistics(slp.problem, "ϕ(x_k))", slp.phi)
    add_statistics(slp.problem, "D(ϕ,p)", slp.directional_derivative)
    add_statistics(slp.problem, "|p|", norm(slp.p, Inf))
    add_statistics(slp.problem, "|J|2", norm(slp.dE, 2))
    add_statistics(slp.problem, "|J|inf", norm(slp.dE, Inf))
    add_statistics(slp.problem, "inf_pr", slp.prim_infeas)
    add_statistics(slp.problem, "mode", mode)
    # add_statistics(slp.problem, "inf_du", dual_infeas)
    add_statistics(slp.problem, "compl", slp.compl)
    add_statistics(slp.problem, "alpha", slp.alpha)
    add_statistics(slp.problem, "iter_time", time() - slp.start_iter_time)
    add_statistics(slp.problem, "time_elapsed", time() - slp.start_time)
end


text:
using Test
using QuantumOpticsBase
using Random, SparseArrays, LinearAlgebra

@testset "embed" begin

Random.seed!(0)

# Set up operators
spinbasis = SpinBasis(1//2)

b1 = NLevelBasis(3)
b2 = SpinBasis(1//2)
b3 = FockBasis(2)

I1 = dense(identityoperator(b1))
I2 = dense(identityoperator(b2))
I3 = dense(identityoperator(b3))

b = b1 ⊗ b2 ⊗ b3

op1 = DenseOperator(b1, b1, rand(ComplexF64, length(b1), length(b1)))
op2 = DenseOperator(b2, b2, rand(ComplexF64, length(b2), length(b2)))
op3 = DenseOperator(b3, b3, rand(ComplexF64, length(b3), length(b3)))

# Test Vector{Int}, Vector{AbstractOperator}
x = embed(b, [1,2], [op1, op2])
y = op1 ⊗ op2 ⊗ I3
@test 0 ≈ abs(tracedistance_nh(x, y))

x = embed(b, [1,2], [sparse(op1), sparse(op2)])
y = op1 ⊗ op2 ⊗ I3
@test 0 ≈ abs(tracedistance_nh(dense(x), y))

x = embed(b, 1, op1)
y = op1 ⊗ I2 ⊗ I3
@test 0 ≈ abs(tracedistance_nh(x, y))

x = embed(b, 2, op2)
y = I1 ⊗ op2 ⊗ I3
@test 0 ≈ abs(tracedistance_nh(x, y))

x = embed(b, 3, op3)
y = I1 ⊗ I2 ⊗ op3
@test 0 ≈ abs(tracedistance_nh(x, y))

# Test Dict(Int=>AbstractOperator)
x = embed(b, Dict(1 => sparse(op1), 2 => sparse(op2)))
y = op1 ⊗ op2 ⊗ I3
@test 0 ≈ abs(tracedistance_nh(dense(x), y))

x = embed(b, Dict(1 => op1, 2 => op2))
y = op1 ⊗ op2 ⊗ I3
@test 0 ≈ abs(tracedistance_nh(x, y))

x = embed(b, Dict([1,3] => sparse(op1⊗op3)))
y = op1 ⊗ I2 ⊗ op3
@test 0 ≈ abs(tracedistance_nh(dense(x), y))

x = embed(b, Dict([1,3] => op1⊗op3))
y = op1 ⊗ I2 ⊗ op3
@test 0 ≈ abs(tracedistance_nh(x, y))

x = embed(b, Dict([3,1] => sparse(op3⊗op1)))
y = op1 ⊗ I2 ⊗ op3
@test 0 ≈ abs(tracedistance_nh(dense(x), y))

x = embed(b, Dict([3,1] => op3⊗op1))
y = op1 ⊗ I2 ⊗ op3
@test 0 ≈ abs(tracedistance_nh(x, y))

end # testset


text:
using RPRMakie
using Test

RPRMakie.activate!(resource=RPR.RPR_CREATION_FLAGS_ENABLE_CPU, iterations=50)
f, ax, pl = meshscatter(rand(Point3f, 100), color=:blue)
out = joinpath(@__DIR__, "recorded")
isdir(out) && rm(out)
mkdir("recorded")
save(joinpath(out, "test.png"), ax.scene);


text:
using LinearAlgebra
using ForwardDiff
using Ipopt
using Plots
using MathOptInterface
const MOI = MathOptInterface

include("kuka_particle_model.jl")

N = 21
NN = (N-1)*(nq+nu+nγ+ns) + 2*nq
#
# q0 = [0.;0.]
# q1 = [0.;0.]
# qN = [1.;1.]

# q0 = [0.;1.;0.;0.0]
# q1 = [0.;1.;0.;0.0]
# qN = [0.;1.;0.;1.1]

q0 = zeros(nq)
q0[1] = 0
q0[3] = 0
q0[4] = -pi/2
q0[5] = 0.
ee_pos0 = Array(end_effector_function(q0[1:7]))
p_pos0 = ee_pos0
p_pos0[3] -= 0.5
q0[8:10] = p_pos0

q1 = copy(q0)

qN = zeros(nq)
qN[1] = 0
qN[3] = 0
qN[4] = -pi/2
qN[5] = 0.

ee_posN = Array(end_effector_function(qN[1:7]))
p_posN = ee_posN
p_posN[3] += 0.1
# p_posN[2] += 0.5
qN[8:10] = p_posN

qD = zeros(nq)
qD[1] = 0
qD[3] = 0
qD[4] = -pi/2
qD[5] = 0.

ee_posD = Array(end_effector_function(qD[1:7]))
p_posD = ee_posD
qD[8:10] = p_posD

Q = [rand(nq) for k = 1:N+1]
U = [rand(nu) for k = 1:N-1]
Γ = [rand(nγ) for k = 1:N-1]
S = [rand(ns) for k = 1:N-1]

function hold_trajectory(n,m,N, mech::Mechanism, q)
    state = MechanismState(mech)
    nn = num_positions(state)
    set_configuration!(state, q[1:nn])
    vd = zero(state.q)
    u0 = dynamics_bias(state)

    if length(q) > m
        throw(ArgumentError("system must be fully actuated to hold an arbitrary position ($(length(q)) should be > $m)"))
    end
    U0 = zeros(m,N)
    for k = 1:N
        U0[:,k] = u0
    end
    return U0
end

U0_hold = hold_trajectory(7,7,N,mechanismkuka,qN[1:7])
U_hold = [U0_hold[:,k]*h for k = 1:N-1]
uu = zeros(7)
uu = 0.01*rand(7)
# uu[1] = 0.1
# uu[3] = 0.1
# uu[5] = 0.1
U_hold_rand = [U0_hold[:,k]*h .+ uu for k = 1:N-1]

function add_rows_cols!(row,col,_r,_c)
    for cc in _c
        for rr in _r
            push!(row,convert(Int,rr))
            push!(col,convert(Int,cc))
        end
    end
end

function packZ(Q,U,Γ,S)
    Z = copy(Q[1])
    for k = 1:N-1
        append!(Z,[Q[k+1];U[k];Γ[k];S[k]])
    end
    append!(Z,Q[N+1])

    return Z
end

unpackZ(Z) = let nq=nq, nu=nu, nγ=nγ, ns=ns, N=N, NN=NN
    shift = nq+nu+nγ+ns
    Q = [k == 1 ? Z[1:nq] : Z[((k-2)*shift+nq) .+ (1:nq)] for k = 1:N+1]
    U = [Z[((k-1)*shift+2nq) .+ (1:nu)] for k = 1:N-1]
    Γ = [Z[((k-1)*shift+2nq+nu) .+ (1:nγ)] for k = 1:N-1]
    S = [Z[((k-1)*shift+2nq+nu+nγ) .+ (1:ns)] for k = 1:N-1]

    return Q,U,Γ,S
end

ZZ = packZ(Q,U,Γ,S)
# Q_,U_,Γ_,S_ = unpackZ(ZZ)
# Q_ == Q
# U_ == U
# Γ_ == Γ
# S_ == S

W = 1.0*Diagonal([ones(7);0.;0.;0.])
WN = 100.0*Diagonal(ones(nq))#Diagonal([10.0*ones(7);1.;1.;1.])
Ru = 1.0*Diagonal(ones(nu))
# Ru[1,1] = 1.0e-3
# Ru[3,3] = 1.0e-3
# Ru[5,5] = 1.0e-3
Rγ = 0.0*Diagonal(ones(nγ))
Rs = 10000.0*Diagonal(ones(1))

Rs*ones(ns)
g(q,u,γ,s,k,qd=qN,scale=diag(W)) = let W=W, Ru=Ru, Rγ=Rγ, Rs=Rs, qN=qN
    (q - qd)'*Diagonal(scale)*(q - qd) + (u - U_hold[k])'*Ru*(u - U_hold[k]) + γ'*Rγ*γ + sum(Rs*s)
end

dgdq(q,u,γ,s,qd=qN,scale=diag(W)) = let W=W, qN=qN
    2*Diagonal(scale)*(q - qd)
end

dgdu(q,u,γ,s,k) = let R=Ru
    2*R*(u - U_hold[k])
end

dgdγ(q,u,γ,s) = let R=Rγ
    2*R*γ
end

dgds(q,u,γ,s) = let R=Rs, ns=ns
    R*ones(ns)
end

gN(q) = let WN=WN, qN=qN
    (q - qN)'*WN*(q - qN)
end

dgNdq(q) = let WN=WN
    2*WN*(q - qN)
end

cost(Z) = let N=N, g=g, gN=gN, W=W, qN=qN
    J = 0.
    Q,U,Γ,S = unpackZ(Z)

    J += (Q[1] - qN)'*W*(Q[1] - qN)

    for k = 1:N-1
        q = Q[k+1]
        u = U[k]
        γ = Γ[k]
        s = S[k]

        qd = qN
        scale = diag(W)
        if k == 11
            qd=copy(qD)
            qd[8] += 0.5
            # scale = zeros(nq)
            scale[8] = 1000.0
        end
        # elseif k == 16
        #     qd=copy(qD)
        #     qd[10] += 0.5
        #     scale = zeros(nq)
        #     scale[10] = 1000.0
        # elseif k == 26
        #     qd=copy(qD)
        #     qd[10] += 0.5
        #     scale = zeros(nq)
        #     scale[8:10] .= 1000.0
        # else
        #     qd=qN
        #     scale = diag(W)
        # end
        J += g(q,u,γ,s,k,qd,scale)
    end

    qN = Q[N+1]
    J += gN(qN)

    return J
end

cost(ZZ)

∇cost!(∇J,Z) = let nq=nq, nu=nu, nγ=nγ, ns=ns, N=N, NN=NN, qN=qN
    J = 0.
    shift = nq+nu+nγ+ns

    Q,U,Γ,S = unpackZ(Z)

    ∇J[1:nq] = 2*W*(Q[1] - qN)

    for k = 1:N-1
        q = Q[k+1]
        u = U[k]
        γ = Γ[k]
        s = S[k]

        qd = qN
        scale = diag(W)
        if k == 11
            qd=copy(qD)
            qd[8] += 0.5
            # scale = zeros(nq)
            scale[8] = 1000.0
        end
        # elseif k == 16
        #     qd=copy(qD)
        #     qd[10] += 0.5
        #     scale = zeros(nq)
        #     scale[10] = 1000.0
        # elseif k == 26
        #     qd=copy(qD)
        #     qd[10] += 0.5
        #     scale = zeros(nq)
        #     scale[8:10] .= 1000.0
        # else
        #     qd=qN
        #     scale = diag(W)
        # end

        ∇J[((k-1)*shift+nq) .+ (1:nq)] = dgdq(q,u,γ,s,qd,scale)
        ∇J[((k-1)*shift + 2nq) .+ (1:nu)] = dgdu(q,u,γ,s,k)
        ∇J[((k-1)*shift + 2nq + nu) .+ (1:nγ)] = dgdγ(q,u,γ,s)
        ∇J[((k-1)*shift + 2nq + nu + nγ) .+ (1:ns)] = dgds(q,u,γ,s)
    end
    ∇J[((N-1)*shift+nq) .+ (1:nq)] = dgNdq(Q[N+1])

    return nothing
end

JJ = zeros(NN)
∇cost!(JJ,ZZ)
JJ

dynamics_constraints!(g,Z) = let nq=nq, N=N
    Q,U,Γ,S = unpackZ(Z)
    shift = nq

    for k = 1:N-1
        q⁺ = Q[k+2]
        q = Q[k+1]
        q⁻ = Q[k]
        u = U[k]
        γ = Γ[k]
        s = S[k]

        g[(k-1)*shift .+ (1:shift)] = var_int(q⁻,q,u,γ,s,q⁺)
    end
    return nothing
end

gg = rand(nq*(N-1))
dynamics_constraints!(gg,ZZ)
gg

∇dynamics_constraints!(jj,Z) = let nq=nq, nu=nu, nγ=nγ, ns=ns, N=N, NN=NN
    Q,U,Γ,S = unpackZ(Z)
    shift = nq*(3nq + nu + nγ + ns)

    for k = 1:N-1
        q⁺ = Q[k+2]
        q = Q[k+1]
        q⁻ = Q[k]
        u = U[k]
        γ = Γ[k]
        s = S[k]
        copyto!(view(jj,(k-1)*shift .+ (1:shift)),vec(∇var_int(q⁻,q,u,γ,s,q⁺)))
    end
    return nothing
end

jj = zeros((N-1)*nq*(3nq + nu + nγ + ns))
∇dynamics_constraints!(jj,ZZ)
jj

sparsity_dynamics(r_shift=0) = let nq=nq, nu=nu, nγ=nγ, ns=ns, N=N, NN=NN
    row = []
    col = []
    r_shift = r_shift
    c_shift = 0
    shift = nq + nu + nγ + ns

    for k = 1:N-1
        r_idx = r_shift .+ (1:nq)

        if k == 1
            c_idx = 1:(shift + 2*nq)
        else
            c_idx = (c_shift) .+ [(1:nq)...,(shift .+ (1:(shift+nq)))...]
        end

        add_rows_cols!(row,col,r_idx,c_idx)

        r_shift += nq
        if k == 1
            c_shift += nq
        else
            c_shift += shift
        end
    end

    collect(zip(row,col))
end

sparsity_dynamics()

(N-1)*(3nq + nu + nγ + ns)*nq

signed_distance_constraints!(g,Z) = let nγ=nγ, N=N
    Q,U,Γ,S = unpackZ(Z)

    for k = 1:N+1
        copyto!(view(g,(k-1)*nγ .+ (1:nγ)),-ϕ(Q[k]))
    end
end

∇signed_distance_constraints!(jj,Z) = let nq=nq, nγ=nγ, N=N
    Q,U,Γ,S = unpackZ(Z)
    shift = nγ*nq

    for k = 1:N+1
        copyto!(view(jj,(k-1)*shift .+ (1:shift)),-vec(∇ϕ(Q[k])))
    end
end

gg = zeros((N+1)*nγ)
signed_distance_constraints!(gg,ZZ)
gg

sparsity_signed_distance(r_shift=0) = let nq=nq, nu=nu, nγ=nγ, ns=ns, N=N
    row = []
    col = []
    r_shift = r_shift
    c_shift = 0
    shift = nq + nu + nγ + ns

    for k = 1:N+1
        r_idx = r_shift .+ (1:nγ)
        c_idx = c_shift .+ (1:nq)

        add_rows_cols!(row,col,r_idx,c_idx)

        r_shift += nγ
        if k == 1
            c_shift += nq
        else
            c_shift += shift
        end
    end

    collect(zip(row,col))
end

sparsity_signed_distance()
(N+1)*(nγ*nq)

complementary_con(γ,s,q⁺) = dot(γ,ϕ(q⁺)) - s[1]
complementary_con(z) = complementary_con(z[1:nγ],z[nγ .+ (1:ns)],z[(nγ+ns) .+ (1:nq)])
∇complementary_con(z) = ForwardDiff.gradient(complementary_con,z)
∇complementary_con(γ,s,q⁺) = ∇complementary_con([γ;s;q⁺])
∇complementary_con(rand(nγ),rand(ns),rand(nq))

complementarity_constraints!(g,Z) = let N=N
    Q,U,Γ,S = unpackZ(Z)

    for k = 1:N-1
        g[k] = complementary_con(Γ[k],S[k],Q[k+2])
    end
end

gg = zeros(N-1)
complementarity_constraints!(gg,ZZ)
gg

∇complementarity_constraints!(jj,Z) = let nq=nq, nγ=nγ, ns=ns, N=N
    Q,U,Γ,S = unpackZ(Z)
    shift = nq+nγ+ns

    for k = 1:N-1
        copyto!(view(jj,(k-1)*shift .+ (1:shift)),∇complementary_con(Γ[k],S[k],Q[k+2]))
    end
end

jj = zeros((N-1)*(nq+nγ+ns))
∇complementarity_constraints!(jj,ZZ)
jj

sparsity_complementarity(r_shift=0) = let nq=nq, nu=nu, nγ=nγ, ns=ns, N=N
    row = []
    col = []
    r_shift = r_shift
    c_shift = 0
    shift = nq + nu + nγ + ns
    idx = nq .+ [((nq+nu) .+ (1:(nγ+ns+nq)))...]

    for k = 1:N-1

        r_idx = r_shift + 1
        c_idx = c_shift .+ idx

        add_rows_cols!(row,col,r_idx,c_idx)

        r_shift += 1
        c_shift += shift
    end
    collect(zip(row,col))
end

sparsity_complementarity()
(N-1)*(nq+nγ+ns)

primal_bounds(q0,q1,qN) = let nq=nq, nu=nu, nγ=nγ, ns=ns, N=N, NN=NN
    Z_low = -Inf*ones(NN)
    Z_up = Inf*ones(NN)

    shift = nq + nu + nγ + ns

    # set γ, s >= 0
    for k = 1:N-1
        Z_low[((k-1)*shift + 2nq + nu) .+ (1:(nγ+ns))] .= 0.
        # Z_up[((k-1)*shift + 2nq + nu) .+ (1:(nγ+ns))] .= 0.
    end

    # # control limits
    # for k = 1:N-1
    #     Z_low[((k-1)*shift + 2nq) .+ (1:(nu))] = [-1.0;0.;-1.0]
    #     Z_up[((k-1)*shift + 2nq) .+ (1:(nu))] = [1.0;0.;1.0]
    # end

    # initial condition
    Z_low[1:nq] = q0
    Z_up[1:nq] = q0
    Z_low[nq .+ (1:nq)] = q1
    Z_up[nq .+ (1:nq)] = q1

    # terminal condition
    Z_low[NN-(nq-1):NN] = qN
    Z_up[NN-(nq-1):NN] = qN

    Z_low, Z_up
end

# primal_bounds(rand(nq),rand(nq))

function constraint_bounds(con)
    con_low = Float64[]
    con_up = Float64[]

    for c in con
        append!(con_up,zeros(c[1]))

        if c[2] == :inequality
            append!(con_low,-Inf*ones(c[1]))
        else
            append!(con_low,zeros(c[1]))
        end
    end
    con_low, con_up
end

# constraint_bounds([(5,:inequality),(6,:eq)])

function line_trajectory(x0::Vector,xf::Vector,N::Int)
    t = range(0,stop=N,length=N)
    slope = (xf .- x0)./N
    x_traj = [slope*t[k] + x0 for k = 1:N]
    x_traj[1] = x0
    x_traj[end] = xf
    x_traj
end


text:
# bound by a lower and upper quantiles
function trim!(x, percentage = 0.01)
    for j = 1:size(x,2)
        q = quantile(x[:,j], percentage)
        x[:,j] = max.(q,x[:,j])
        q = quantile(x[:,j], 1.0 - percentage)
        x[:,j] = min.(q,x[:,j])
    end
end

max_stars_repo_name: nilshg/Econometrics

max_issues_repo_name: nilshg/Econometrics

max_forks_repo_name: nilshg/Econometrics


text:
LIBINT_VERSION = "2.2.0-beta1"
lib_path = joinpath(dirname(@__FILE__),"usr","lib")

@static if is_linux()
	try begin
			run(`grep --ignore-case --silent 'avx' /proc/cpuinfo`)
			SIMD = "avx"
			cd(joinpath(dirname(@__FILE__),"usr","lib"))
			run(`ln -sf libint2jl-$LIBINT_VERSION-$SIMD.so libint2jl.so`)
			run(`ln -sf libint2-$LIBINT_VERSION-$SIMD.so libint2-QuantumLab.so`)
		end
	catch
		warn("Your system (CPU) does not support AVX.\n
							The Libint2 Library for Linux is only provided as AVX enabled version.\n
							Therefore a Libint2 stub is used.\n
							You should be able to start QuantumLab and perform calculation but
							with a drastically decreased performance.\n
							You can try to build your own version of libint2.so and libint2jl.so.
							As a reference take a look at the Makefile in QuantumLab/deps/usr/src.")
	end
end

@static if is_apple()
	try begin
			run(pipeline(`sysctl machdep.cpu.features`,`grep --ignore-case  --silent 'avx\|sse'`))
			SIMD = "avx"
			try run(pipeline(`sysctl machdep.cpu.features`,`grep --ignore-case --silent 'avx'`))
			catch
				SIMD = "sse"
			end
			cd(joinpath(dirname(@__FILE__),"usr","lib"))
			run(`ln -sf libint2jl-$LIBINT_VERSION-$SIMD.dylib libint2jl.so`)
			run(`ln -sf libint2-$LIBINT_VERSION-$SIMD.dylib libint2-QuantumLab.so`)
		end
	catch
		warn("Your CPU does support neither SSE nor AVX
							The Libint2 Library for Linux is only provided as AVX enabled version.\n
							Therefore a Libint2 stub is used.\n
							You should be able to start QuantumLab and perform calculation but
							with a drastically decreased performance.\n
							You can try to build your own version of libint2.dylib and libint2jl.dylib.
							As a reference take a look at the Makefile in QuantumLab/deps/usr/src.")
	end
end


text:
using LyapunovExponents

function get_angles(solver)
    return [acos(abs(dot(C[:, 1], C[:, 2]))) * 2 / π for C
            in backward_dynamics!(solver; progress=1)]
end

henon_demo = LyapunovExponents.henon_map(t_attr=1000000)
henon_demo.prob :: LEProblem
henon_prob = CLVProblem(henon_demo.prob)  # convert it to CLVProblem
henon_angles = @time get_angles(init(henon_prob))

lozi_prob = CLVProblem(LyapunovExponents.lozi_map(t_attr=1000000).prob)
lozi_angles = @time get_angles(init(lozi_prob))

using Plots
plt = plot(xlabel="Angle [pi/2 rad]", ylabel="Density", legend=:topleft)
stephist!(plt, henon_angles,
          bins=1000, normalize=true, linecolor=1, label="Henon")
stephist!(twinx(plt), lozi_angles,
          bins=1000, normalize=true, linecolor=2, label="Lozi")
plt


text:
for a = 1:1000
  for b = (a+1):1000
    c = sqrt(a^2 + b^2)
    if c == floor(c) && a + b + c == 1000
      c = Int(floor(c))
      println(a * b * c)
    end
  end
end

max_stars_repo_name: sujimodern/project-euler-jl

max_issues_repo_name: sujimodern/project-euler-jl

max_forks_repo_name: sujimodern/project-euler-jl


text:
module Misfits

using ForwardDiff
using Distances
using StatsBase
using LinearAlgebra
using LinearMaps

"""
Distance between x[:,i] and x[:,j]
cosine_dist(x, y)=1 - dot(x, y) / (norm(x) * norm(y))
corr_dist(x, y)=cosine_dist(x - mean(x), y - mean(y))
"""
function error_pairwise_corr_dist(dfdx,x)
	R=pairwise(CorrDist(), x)
	J=sum(R)

	if(!(dfdx===nothing))
		f(x)=error_pairwise_corr_dist(nothing,x)
		ForwardDiff.gradient!(dfdx, f, x);
	end
	return J
end

function error_corr_dist!(dfdx, x, y)
	xn=norm(x)
	yn=norm(y)

	dotxy=0.0
	for i in eachindex(x)
		dotxy+=x[i]*y[i]
	end

	J = 1. - (dotxy/xn/yn)

	if(!(dfdx===nothing))
		for i in eachindex(x)
			dfdx[i]=(dotxy*yn*x[i]/xn - y[i]*yn*xn)/(xn*yn)/(xn*yn)
		end
	end
	return J
end

function error_autocorr_pairwise_corr_dist(dfdx,x)
	X=fft(x, [1]);
	# remove phase
	#XC=real(X .* conj(X))
	J=error_pairwise_corr_dist(nothing, X)
	J=norm(J)

	if(!(dfdx===nothing))
		GXC=similar(X)
		error_pairwise_corr_dist(GXC,X)
		#GX = GXC .* conj(X)
		dfdx[:]=real.(ifft(GXC, [1]))
	end

	return J
end

function error_weighted_norm(dfdx,x,w)
	J=0.0
	for i in eachindex(x)
		J += (x[i]*w[i])*(x[i]*w[i])
	end
	J=sqrt(J)
	if(!(dfdx===nothing))
		f(x)=error_weighted_norm(nothing,x,w)
		ForwardDiff.gradient!(dfdx, f, x);
	end
	return J
end

"""
Compute the error with a normalized vector ̂̂x
̂x=x/|x|

function
* `ghat` : gradient with respect to ̂x
* `g`	 : output gradient with respect to x
* `X`	 : preallocated matrix of dimension (nx,nx) if necessary
"""
function derivative_vector_magnitude!(g,ghat,x,X=nothing)
	xn=norm(x)
	nx=length(x)
	rmul!(x, inv(xn))

	# compute the outer product of
	if(!(X===nothing))
		mul!(X,x,transpose(x))
	else
		X=x*transpose(x)
	end
	for i in 1:nx
		X[i,i]=X[i,i]-1.
	end
	rmul!(X,-inv(xn))
	mul!(g,X,ghat)
	rmul!(x, xn)
end

"""
Normalized least-squares error between two arrays after
estimating a scalar that best fits on to another.
Return misfit and α such that αx-y is minimum.
Normalization is done with respect to the
norm of y.
"""
function error_after_scaling(
			     x::AbstractArray{T},
			     y::AbstractArray{T}
			     ) where T
	any(size(x) ≠ size(y)) && error("x and y different sizes")
	sxx=T(0.0)
	sxy=T(0.0)
	for i in eachindex(x)
		sxx += x[i]*x[i]
		sxy += x[i]*y[i]
	end
	α = sxy * inv(sxx)

	if(!(iszero(α)))
		rmul!(x, α)
		J = error_squared_euclidean!(nothing,  x,   y,   nothing, norm_flag=true)
		rmul!(x, inv(α))
	else
		J = zero(T)
	end

	return J, α
end

"""
Measure the least-squares distance between auto correlations of x and y,
"""
function error_after_normalized_autocor(x::AbstractArray{Float64}, y::AbstractArray{Float64})
	nt=size(x,1)

	ax=autocor(x, 0:nt-1, demean=true) # normalizd autocorr from StatsBase
	ay=autocor(y, 0:nt-1, demean=true)

	return error_squared_euclidean!(nothing,  ax,   ay, nothing, norm_flag=true)
end

"""
Compute the L2 distance between two arrays: ``\\sqrt{\\sum_{i=1}^n |a_i - b|^2}``.
* `norm_flag` : optional; normalizes the distance with the norm of `y` (don't use in the inner loops)
"""
function error_squared_euclidean!(dfdx,  x,   y,   w; norm_flag=false)
	J=zero(eltype(x))
	if(w===nothing)
		for i in eachindex(x)
			@inbounds J += abs2(x[i]-y[i])
		end
	else
		for i in eachindex(x)
			@inbounds J += w[i] * abs2(x[i]-y[i])
		end
	end
	if(!(dfdx === nothing))
		if(w===nothing)
			for i in eachindex(x)
				@inbounds dfdx[i] = 2.0 * (x[i]-y[i])
			end
		else
			for i in eachindex(x)
				@inbounds dfdx[i] = 2.0 * w[i] * (x[i]-y[i])
			end
		end
	end
	if(norm_flag) # second norm of y
		ynorm=zero(Float64)
		for i in eachindex(y)
			@inbounds ynorm += abs2(y[i])
		end
		# divide the functional with ynorm
		J /= ynorm
		if(!(dfdx === nothing))
			rmul!(dfdx, inv(ynorm)) # divide dfdx with ynorm too
		end
	end
	return J
end

"""
Computed the weighted norm of
"""
function error_weighted_norm!(dfdx,  x,   w)
	J=zero(Float64)
	if(w===nothing)
		for i in eachindex(x)
			J += (x[i]) * (x[i])
		end
	else
		for i in eachindex(x)
			J += w[i] * (x[i]) * (x[i])
		end
	end
	if(!(dfdx === nothing))
		if(w===nothing)
			for i in eachindex(x)
				dfdx[i] = 2.0 * (x[i])
			end
		else
			for i in eachindex(x)
				dfdx[i] = 2.0 * w[i] * (x[i])
			end
		end
	end
	return J
end

"""
Calculate the front load of dfdx
"""
function front_load!(dfdx,  x::AbstractMatrix)
	nt=size(x,1)
	nr=size(x,2)
	J=zero(Float64)
	for ir in 1:nr
		for it in 1:nt
			J += abs2(x[it,ir] * inv(nt-1) * (it-1))
		end
	end
	if(!(dfdx === nothing))
		for ir in 1:nr
			for it in 1:nt
				dfdx[it,ir] = 2.0 * (x[it,ir]) *
					abs2(inv(nt-1) * (it-1))
			end
		end
	end
	return J
end

"""
Compute error b/w g1 and g2 after ingnoring time translation
allocates a lot of memory, dont use for big data
"""
function error_after_translation(g1,g2)
	err=[]
	for is in 1:size(g1,1)
		g11=circshift(g1,(is,0))
		push!(err,Misfits.error_after_scaling(g11,g2)[1])
	end
	return minimum(err)
end

"""
type for computing the generalized least-squares error
* `Q` : must be symmetric, inverse covariance matrix
"""
mutable struct P_gls{T}
	Q::LinearMaps.LinearMap{T} # inverse covariance matrix, must be symmetric
	r::Vector{T} # store residual temporarily
	Qr::Vector{T} # store Q*r temporarily
end

"""
# Q is a diagonal matrix with α
* `n` : size
"""
function P_gls(n,α)
	Q=LinearMap(
	     (y,x)->LinearAlgebra.mul!(y,x,α),
	     (y,x)->LinearAlgebra.mul!(y,x,α),
		      n, n, ismutating=true, issymmetric=true)
	r=zeros(n)
	Qr=zeros(n)
	return P_gls(Q,r,Qr)
end

function P_gls(Q)
	n=size(Q,1)
	r=zeros(n)
	Qr=zeros(n)
	return P_gls(Q,r,Qr)
end

function func_grad!(dfdx,  x,   y, pa::P_gls)
	# compute the difference
	for i in eachindex(x)
		pa.r[i]=x[i]-y[i]
	end

	mul!(pa.Qr, pa.Q, pa.r)

	J=pa.r' * pa.Qr

	if(!(dfdx === nothing))
		mul!(dfdx,pa.Q,pa.r)
		rmul!(dfdx,2.0)
	end
	return J
end

end # module


text:
# Repeated-measures Designs
#
# author: yizhan miao
# email: yzmiao@protonmail.com
# last update: Oct 26 2018

"""One-way Repeated Measures ANOVA

Syntax:
- `RESULT, _F, _p = rm_anova(Y::Array{T, 1}, S::Array{D, 1}, F::Array{D, 1};name::Symbol=:F, preview::Bool=true) where {T <: Number, D <: Integer}`

Arguments:
- `Y`: data set in 1d vector
- `S`: corresponding subject marker, as Array{Integer, 1}
- `F`: corresponding factor marker, as Array{Integer, 1}

Keywords:
- `name`: the name of the factor [default: `:F`]
- `preview`: print anova table [default: `true`]

Returns:
- ANOVA Result
- F statistic
- p value

Example:
```julia
julia> rm_anova_Y = [
    -4.4297,  4.7513,  3.2971, -1.4606,
     4.8458,  5.1163,  6.0739, -1.9225,
     6.1542, 10.4794, 12.4438,  9.2150];
julia> rm_anova_S = [1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4];
julia> rm_anova_F = [1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3];
julia> ANOVA.rm_anova(rm_anova_Y, rm_anova_S, rm_anova_F, name=:factor);
```
```
One-way repeated measures ANOVA
-----------------------------------------------
source  SS      df    MS      F       Prob
-----------------------------------------------
factor  169.437 2     84.718  12.364  0.00744
error   41.112  6     6.852   NaN     NaN
S       74.253  3     24.751  NaN     NaN
total   284.802 11    NaN     NaN     NaN
-----------------------------------------------
```
"""
function rm_anova(Y::Array{T, 1}, S::Array{D, 1}, F::Array{D, 1};name::Symbol=:F, preview::Bool=true) where {T <: Number, D <: Integer}
    # check balance
    _factor = unique(F)
    _n_factor = zeros(Int64, length(_factor))
    for (idx, each) in enumerate(_factor)
        _n_factor[idx] = F[F .== each] |> length
    end

    _subject = unique(S)
    _n_subject = zeros(Int64, length(_subject))
    for (idx, each) in enumerate(_subject)
        _n_subject[idx] = S[S .== each] |> length
    end

    @assert(length(unique(_n_factor)) == 1, "input data is not balanced across factors.")
    @assert(length(unique(_n_subject)) == 1, "input data is not balanced across subject: $(unique(_n_subject)).")
    # TODO: support unbalanced ANOVA

    _mean_total = mean(Y)
    _mean_treat = [mean(Y[F .== _F]) for _F in _factor]
    _mean_subject = [mean(Y[S .== _S]) for _S in _subject]

    _SS_total = (Y .- _mean_total).^2 |> sum
    _SS_treat = sum((_mean_treat .- _mean_total).^2) * length(_subject)
    _SS_subject = sum((_mean_subject .- _mean_total) .^ 2) * length(_factor)
    _SS_SF = _SS_total - _SS_treat - _SS_subject

    _df_total = length(Y) - 1
    _df_treat = length(_factor) - 1
    _df_subject = length(_subject) - 1
    _df_SF = _df_total - _df_treat - _df_subject

    ##### rm anova #####
    _MS_treat = _SS_treat / _df_treat
    _MS_subject = _SS_subject / _df_subject
    _MS_SF = _SS_SF / _df_SF

    _F = _MS_treat / _MS_SF
    _F_dist = FDist(_df_treat, _df_SF)

    _p = ccdf(_F_dist, _F)

    TOTAL = anova_entry(:total, _SS_total, _df_total, NaN, NaN, NaN)
    TREAT = anova_entry(name, _SS_treat, _df_treat, _MS_treat, _F, _p)
    SF    = anova_entry(:error, _SS_SF, _df_SF, _MS_SF, NaN, NaN)
    SUBJECT = anova_entry(:S, _SS_subject, _df_subject, _MS_subject, NaN, NaN)

    RESULT = anova_result("One-way repeated measures ANOVA", [TREAT, SF, SUBJECT, TOTAL])

    preview&&preview_anova_result(RESULT)

    RESULT, _F, _p

end

"""Two-way Repeated Measures ANOVA

Syntax:
- `Result = rm_anova2(Y::Array{T, 1}, S::Array{D, 1}, A::Array{D, 1}, B::Array{D, 1};names::Tuple{Symbol, Symbol}=(:A, :B), preview::Bool=true) where {T <: Number, D <: Integer}`

Arguments:
- `Y`: data set in 1d vector
- `S`: corresponding subject marker, as Array{Integer, 1}
- `A`: corresponding factor A marker, as Array{Integer, 1}
- `B`: corresponding factor B marker, as Array{Integer, 1}

Keywords:
- `names`: the name of the factor [default: `(:A, :B)`]
- `preview`: print anova table [default: `true`]

Return:
- ANOVA Result

"""
function rm_anova2(Y::Array{T, 1}, S::Array{D, 1}, A::Array{D, 1}, B::Array{D, 1};names::Tuple{Symbol, Symbol}=(:A, :B), preview::Bool=true) where {T <: Number, D <: Integer}
    # check balance
    _factor_A = unique(A)
    _n_factor_A = zeros(Int64, length(_factor_A))
    for (idx, each) in enumerate(_factor_A)
        _n_factor_A[idx] = A[A .== each] |> length
    end

    _factor_B = unique(B)
    _n_factor_B = zeros(Int64, length(_factor_B))
    for (idx, each) in enumerate(_factor_B)
        _n_factor_B[idx] = B[B .== each] |> length
    end

    _subject = unique(S)
    _n_subject = zeros(Int64, length(_subject))
    for (idx, each) in enumerate(_subject)
        _n_subject[idx] = S[S .== each] |> length
    end

    @assert(length(unique(_n_factor_A)) == 1, "input data is not balanced across factor A.")
    @assert(length(unique(_n_factor_B)) == 1, "input data is not balanced across factor B.")
    @assert(length(unique(_n_subject)) == 1,  "input data is not balanced across subject.")
    #TODO: support unbalanced ANOVA

    _mean_total = mean(Y)
    _mean_A = [mean(Y[A .== _A]) for _A in _factor_A]
    _mean_B = [mean(Y[B .== _B]) for _B in _factor_B]
    _mean_BA = [mean(Y[(A .== _A) .& (B .== _B)]) for _A in _factor_A, _B in _factor_B]
    _mean_S = [mean(Y[S .== _S]) for _S in _subject]

    _mean_AS = [mean(Y[(A .== _A) .& (S .== _S)]) for _A in _factor_A, _S in _subject]
    _mean_BS = [mean(Y[(B .== _B) .& (S .== _S)]) for _B in _factor_B, _S in _subject]

    _SS_total = (Y .- _mean_total).^2 |> sum
    _SS_A = sum((_mean_A .- _mean_total) .^ 2) * length(_subject) * length(_factor_B)
    _SS_B = sum((_mean_B .- _mean_total) .^ 2) * length(_subject) * length(_factor_A)
    _SS_BA = sum((_mean_BA .- _mean_total) .^ 2) * length(_subject) - _SS_A - _SS_B
    _SS_S = sum((_mean_S .- _mean_total).^2) * length(_factor_A) * length(_factor_B)
    _SS_AS = sum((_mean_AS .- _mean_total) .^ 2) * length(_factor_B) - _SS_A - _SS_S
    _SS_BS = sum((_mean_BS .- _mean_total) .^ 2) * length(_factor_A) - _SS_B - _SS_S
    _SS_BAS = _SS_total - _SS_A - _SS_B - _SS_BA - _SS_S - _SS_AS - _SS_BS

    _df_total = length(Y) - 1
    _df_A = length(_factor_A) - 1
    _df_B = length(_factor_B) - 1
    _df_S = length(_subject) - 1
    _df_BA = _df_A * _df_B
    _df_AS = _df_A * _df_S
    _df_BS = _df_B * _df_S
    _df_BAS = _df_BA * _df_S

    _MS_A = _SS_A / _df_A
    _MS_B = _SS_B / _df_B
    _MS_BA = _SS_BA / _df_BA
    _MS_S = _SS_S / _df_S
    _MS_AS = _SS_AS / _df_AS
    _MS_BS = _SS_BS / _df_BS
    _MS_BAS = _SS_BAS / _df_BAS

    ##### ANOVA #####
    _F_A = _MS_A / _MS_AS
    _F_B = _MS_B / _MS_BS
    _F_BA = _MS_BA / _MS_BAS

    _p_A = ccdf(FDist(_df_A, _df_AS), _F_A)
    _p_B = ccdf(FDist(_df_B, _df_BS), _F_B)
    _p_BA = ccdf(FDist(_df_BA, _df_BAS), _F_BA)

    TOTAL = anova_entry(:total, _SS_total, _df_total, NaN, NaN, NaN)
    _S    = anova_entry(:S, _SS_S, _df_S, _MS_S, NaN, NaN)
    _A    = anova_entry(names[1], _SS_A, _df_A, _MS_A, _F_A, _p_A)
    _B    = anova_entry(names[2], _SS_B, _df_B, _MS_B, _F_B, _p_B)
    _BA   = anova_entry(Symbol("$(names[2])*$(names[1])"), _SS_BA, _df_BA, _MS_BA, _F_BA, _p_BA)
    _AS   = anova_entry(Symbol("S*$(names[1])"), _SS_AS, _df_AS, _MS_AS, NaN, NaN)
    _BS   = anova_entry(Symbol("S*$(names[2])"), _SS_BS, _df_BS, _MS_BS, NaN, NaN)
    _BAS   = anova_entry(Symbol("S*$(names[2])*$(names[1])"), _SS_BAS, _df_BAS, _MS_BAS, NaN, NaN)

    RESULT = anova_result("Two-way repeated measures ANOVA", [_S, _A, _B, _BA, _AS, _BS, _BAS, TOTAL])

    preview&&preview_anova_result(RESULT)

    RESULT
end


text:
module ExNonpar
using SDE
using SDE.Schauder
require(Pkg.dir("SDE","src", "SDEPlot.jl"))
using SDEPlot

bex1(s,x) = 4(sin(2pi*x) - x)
bex2(s,x) = -4min(x+.5, 0) - 4max(x-.5, 0)
si(s,x) = 1.
range(a) = (minimum(a), maximum(a))

function example(L, logN,T, f, B; VERBOSE = 0)

	global y, t
	srand(42)
	beta = 0.5
	a = -2
	b = 2
	K = 0
	if (B=="B1"); K = 1;end # Schauder basis with constant
	if (B=="B2"); K = 2;end # Schauder basis plus two affine component
				# otherwise Schauder basis

	VERBOSE > 0 && println("L $L, beta $beta")
	N = 10^logN
	dt = T/N
	P = Diffusion{1}(f, si, ())
    tt = linspace(0., T, N+1)
    W = sample(tt, Wiener())
    u = 0.5(b+a)
    y = euler(u, W, P).yy

	if (B == "B1")
		truec = fe_transfB1(x-> f(0., x), a,b, L)
	elseif (B == "B2")
		truec = fe_transfB2(x-> f(0., x), a,b, L)
	else
		truec = fe_transf(x-> f(0., x), a,b, L)

	end
    VERBOSE > 0 && println("true drift f.e. coefficients",truec')

	println("range y $(range(y))")

	post = bayes_drift(y, dt, a, b, L, 10ones(K), beta , B)

	pl = visualize_posterior(post,x -> f(0.,x), 1.96)

	(post, pl)
end
println("Try: post = ExNonpar.example(5,5, 100, ExNonpar.bex1, \"B2\"); Winston.display(post[2])")

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 = "rocPRIM"
version = v"4.2.0"

# Collection of sources required to build
sources = [
    ArchiveSource("https://github.com/ROCmSoftwarePlatform/rocPRIM/archive/rocm-$(version).tar.gz",
                  "3932cd3a532eea0d227186febc56747dd95841732734d9c751c656de9dd770c8"),
]

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

amdgpu_targets="gfx900,gfx906,gfx908,gfx1010,gfx1011,gfx1012"

ln -s ${prefix}/bin/clang ${prefix}/tools/clang
mkdir build && cd build
export ROCM_PATH=${prefix}
export HIP_CLANG_PATH=$WORKSPACE/destdir/tools
export HIP_CLANG_HCC_COMPAT_MODE=1
export HIP_RUNTIME=rocclr
export HIP_COMPILER=clang
export HIP_PLATFORM=amd
export HIPCC_VERBOSE=1
cmake -DCMAKE_PREFIX_PATH=${prefix} \
      -DCMAKE_INSTALL_PREFIX=${prefix} \
      -DCMAKE_BUILD_TYPE=Release \
      -DCMAKE_CXX_COMPILER=${prefix}/hip/bin/hipcc \
      -DAMDGPU_TARGETS=$(echo $amdgpu_targets | tr ',' ';') \
      ..
make -j${nproc}
make install
rm ${prefix}/tools/clang
"""

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

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

# Dependencies that must be installed before this package can be built
dependencies = [
    BuildDependency(PackageSpec(; name="ROCmLLVM_jll", version=v"4.2.0")),
    Dependency("hsa_rocr_jll"; compat="4.2.0"),
    Dependency("ROCmCompilerSupport_jll"; compat="4.2.0"),
    Dependency("ROCmOpenCLRuntime_jll"; compat="4.2.0"),
    Dependency("HIP_jll"; compat="4.2.0"),
    Dependency("rocm_cmake_jll"),
    Dependency("rocminfo_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies,
               julia_compat="1.7",
               preferred_gcc_version=v"8", preferred_llvm_version=v"11")

max_stars_repo_name: jakebolewski/Yggdrasil

max_issues_repo_name: jakebolewski/Yggdrasil

max_forks_repo_name: jakebolewski/Yggdrasil


text:
@adjoint (::Type{T})(args...) where T<:Array = T(args...), Δ -> nothing

@nograd size, length, eachindex, Colon(), findfirst, randn, ones, zeros, one, zero

@adjoint Base.vect(xs...) = Base.vect(xs...), Δ -> (Δ...,)

@adjoint copy(x::AbstractArray) = copy(x), ȳ -> (ȳ,)

Base.zero(xs::AbstractArray{Any}) = fill!(similar(xs), nothing)

@adjoint function getindex(xs::Array, i...)
  xs[i...], function (Δ)
    Δ′ = zero(xs)
    Δ′[i...] = Δ
    (Δ′, map(_ -> nothing, i)...)
  end
end

@adjoint! setindex!(xs::AbstractArray, x...) = setindex!(xs, x...),
  _ -> error("Mutating arrays is not supported")

# General

@adjoint collect(x) = collect(x), Δ -> (Δ,)

@adjoint permutedims(xs, dims) = permutedims(xs, dims),
  Δ -> (permutedims(Δ, invperm(dims)), nothing)

@adjoint reshape(xs, dims...) = reshape(xs, dims...),
  Δ -> (reshape(Δ, size(xs)),map(_->nothing,dims)...)

@adjoint function hvcat(rows::Tuple{Vararg{Int}}, xs::T...) where T<:Number
  hvcat(rows, xs...), ȳ -> (nothing, ȳ...)
end

pull_block_vert(sz, Δ, A::AbstractVector) = Δ[sz-length(A)+1:sz]
pull_block_vert(sz, Δ, A::AbstractMatrix) = Δ[sz-size(A, 1)+1:sz, :]
@adjoint function vcat(A::Union{AbstractVector, AbstractMatrix}...)
  sz = cumsum([size.(A, 1)...])
  return vcat(A...), Δ->(map(n->pull_block_vert(sz[n], Δ, A[n]), eachindex(A))...,)
end

pull_block_horz(sz, Δ, A::AbstractVector) = Δ[:, sz]
pull_block_horz(sz, Δ, A::AbstractMatrix) = Δ[:, sz-size(A, 2)+1:sz]
@adjoint function hcat(A::Union{AbstractVector, AbstractMatrix}...)
  sz = cumsum([size.(A, 2)...])
  return hcat(A...), Δ->(map(n->pull_block_horz(sz[n], Δ, A[n]), eachindex(A))...,)
end

@adjoint function repeat(xs; inner=ntuple(_->1, ndims(xs)), outer=ntuple(_->1, ndims(xs)))
  repeat(xs, inner = inner, outer = outer), function (Δ)
    Δ′ = zero(xs)
    S = size(xs)

    # Loop through each element of Δ, calculate source dimensions, accumulate into Δ′
    for (dest_idx, val) in pairs(IndexCartesian(), Δ)
        # First, round dest_idx[dim] to nearest gridpoint defined by inner[dim], then
        # wrap around based on original size S.
        src_idx = [mod1(div(dest_idx[dim] - 1, inner[dim]) + 1, S[dim]) for dim in 1:length(S)]
        Δ′[src_idx...] += val
    end
    return (Δ′,)
  end
end

# Reductions

@adjoint function sum(xs::AbstractArray; dims = :)
  if dims === (:)
    sum(xs), Δ -> (FillArray(Δ, size(xs)),)
  else
    sum(xs, dims = dims), Δ -> (similar(xs) .= Δ,)
  end
end

function _forward(cx::Context, ::typeof(sum), f, xs::AbstractArray)
  y, back = forward(cx, (xs -> sum(f.(xs))), xs)
  y, ȳ -> (nothing, nothing, back(ȳ)...)
end

@adjoint function prod(xs; dims = :)
  if dims === (:)
    prod(xs), Δ -> (prod(xs) ./ xs .* Δ,)
  else
    prod(xs, dims = dims),
      Δ -> (reshape(.*(circshift.([reshape(xs, length(xs))], 1:length(xs)-1)...), size(xs)) .* Δ,)
  end
end

@adjoint function maximum(xs; dims = :)
  max, i = findmax(xs, dims = dims)
  max, function (Δ)
    Δ isa Real && Δ <= sqrt(eps(float(Δ))) && return nothing
    Δ′ = zero(xs)
    Δ′[i] = Δ
    return (Δ′,)
  end
end

@adjoint function minimum(xs; dims = :)
  min, i = findmin(xs, dims = dims)
  min, function (Δ)
    Δ′ = zero(xs)
    Δ′[i] = Δ
    return (Δ′,)
  end
end

# LinAlg

@adjoint function(a::AbstractVecOrMat * b::AbstractVecOrMat)
  return a * b, function(Δ)
    return (reshape(Δ * transpose(b), size(a)), reshape(transpose(a) * Δ, size(b)))
  end
end

@adjoint transpose(x) = transpose(x), Δ -> (transpose(Δ),)
@adjoint Base.adjoint(x) = x', Δ -> (Δ',)
@adjoint parent(x::LinearAlgebra.Adjoint) = parent(x), ȳ -> (LinearAlgebra.Adjoint(ȳ),)

function _kron(mat1::AbstractMatrix,mat2::AbstractMatrix)
    m1, n1 = size(mat1)
    mat1_rsh = reshape(mat1,(1,m1,1,n1))

    m2, n2 = size(mat2)
    mat2_rsh = reshape(mat2,(m2,1,n2,1))

    return reshape(mat1_rsh.*mat2_rsh, (m1*m2,n1*n2))
end

@adjoint kron(a::AbstractMatrix, b::AbstractMatrix) = forward(_kron, a, b)

@adjoint iterate(r::UnitRange, i...) = iterate(r, i...), _ -> nothing

@adjoint diag(A::AbstractMatrix) = diag(A), Δ->(Diagonal(Δ),)

@adjoint function \(A::AbstractMatrix, B::AbstractVecOrMat)
  Y = A \ B
  return Y, function(Ȳ)
      B̄ = A' \ Ȳ
      return (-B̄ * Y', B̄)
  end
end

@adjoint function /(A::AbstractMatrix, B::AbstractMatrix)
  Y = A / B
  return Y, function(Ȳ)
      Ā = Ȳ / B'
      return (Ā, -Y' * Ā)
  end
end

_symmetric_back(Δ) = UpperTriangular(Δ) + LowerTriangular(Δ)' - Diagonal(Δ)
_symmetric_back(Δ::UpperTriangular) = Δ
@adjoint function Symmetric(A::AbstractMatrix)
  back(Δ::AbstractMatrix) = (_symmetric_back(Δ),)
  back(Δ::NamedTuple) = (_symmetric_back(Δ.data),)
  return Symmetric(A), back
end

# Implementation due to Seeger, Matthias, et al. "Auto-differentiating linear algebra."
@adjoint function cholesky(Σ::Union{StridedMatrix, Symmetric{<:Real, <:StridedMatrix}})
  C = cholesky(Σ)
  return C, function(Δ)
    U, Ū = C.U, Δ.factors
    Σ̄ = Ū * U'
    Σ̄ = copytri!(Σ̄, 'U')
    Σ̄ = ldiv!(U, Σ̄)
    BLAS.trsm!('R', 'U', 'T', 'N', one(eltype(Σ)), U.data, Σ̄)
    @inbounds for n in diagind(Σ̄)
      Σ̄[n] /= 2
    end
    return (UpperTriangular(Σ̄),)
  end
end

@adjoint function cholesky(Σ::Real)
  C = cholesky(Σ)
  return C, Δ::NamedTuple->(Δ.factors[1, 1] / (2 * C.U[1, 1]),)
end

# Various sensitivities for `literal_getproperty`, depending on the 2nd argument.
@adjoint function literal_getproperty(C::Cholesky, ::Val{:uplo})
  return literal_getproperty(C, Val(:uplo)), function(Δ)
    return ((uplo=nothing, info=nothing, factors=nothing),)
  end
end
@adjoint function literal_getproperty(C::Cholesky, ::Val{:info})
  return literal_getproperty(C, Val(:info)), function(Δ)
    return ((uplo=nothing, info=nothing, factors=nothing),)
  end
end
@adjoint function literal_getproperty(C::Cholesky, ::Val{:U})
  return literal_getproperty(C, Val(:U)), function(Δ)
    Δ_factors = C.uplo == 'U' ? UpperTriangular(Δ) : LowerTriangular(copy(Δ'))
    return ((uplo=nothing, info=nothing, factors=Δ_factors),)
  end
end
@adjoint function literal_getproperty(C::Cholesky, ::Val{:L})
  return literal_getproperty(C, Val(:L)), function(Δ)
    Δ_factors = C.uplo == 'L' ? LowerTriangular(Δ) : UpperTriangular(copy(Δ'))
    return ((uplo=nothing, info=nothing, factors=Δ_factors),)
  end
end

@adjoint function logdet(C::Cholesky)
  return logdet(C), function(Δ)
    return ((uplo=nothing, info=nothing, factors=Diagonal(2 .* Δ ./ diag(C.factors))),)
  end
end

@adjoint function +(A::AbstractMatrix, S::UniformScaling)
  return A + S, Δ->(Δ, (λ=sum(view(Δ, diagind(Δ))),))
end


text:
for (elty, ext) in ((:Float32, :s),
                    (:Float64, :d),
                    (:ComplexF32, :c),
                    (:ComplexF64, :z))
    for mattype in ("", "Dist")
        mat = Symbol(mattype, "Matrix")
        @eval begin

            # Hessenberg
            function _hessenberg!(uplo::UpperOrLower, A::$mat{$elty}, t::$mat{$elty})
                ElError(ccall(($(string("ElHessenberg", mattype, "_", ext)), libEl), Cuint,
                    (UpperOrLower, Ptr{Cvoid}, Ptr{Cvoid}),
                    uplo, A.obj, t.obj))
                return A, t
            end
            _hessenberg!(A::$mat{$elty}) = _hessenberg!(UPPER, A, $mat($elty))
        end
    end
end

struct ElHessenberg{T,S<:ElementalMatrix} <: LinearAlgebra.Factorization{T}
    factors::S
    τ::S
    ElHessenberg{T,S}(factors::ElementalMatrix{T}, τ::ElementalMatrix{T}) where {T,S<:ElementalMatrix} = new{T,S}(factors, τ)
end

ElHessenberg(factors::ElementalMatrix{T}, τ::ElementalMatrix{T}) where {T} = ElHessenberg{T,typeof(factors)}(factors, τ)

LinearAlgebra.hessenberg!(A::ElementalMatrix) = ElHessenberg(_hessenberg!(A)...)