text:
include("../examples/solenoid.jl")
using Test
using Zygote
using PyPlot
function test_attractor()
	s = [2.0, 4.0, 0.0]
	u = rand(3)
	n_runup = 2000
	for i = 1:n_runup
		u = next(u,s)
	end
	@show u
	n = 2000000
	u_trj = step(u, s, n)
	u_trj = u_trj[:,1:10:end]
	x, y, z = view(u_trj,1,:), view(u_trj,2,:), view(u_trj,3,:)
	fig, ax = subplots(1,1)
	ax.plot(x,y,".",ms=2)
	#ax.set_xlabel("x",fontsize=28)
	#ax.set_ylabel("y",fontsize=28)
	ax.tick_params(axis="both",labelsize=28)
	ax.axis("scaled")
	fig, ax = subplots(1,1)
	ax.plot(x,z,".",ms=2)
	ax.set_xlabel("x",fontsize=28)
	ax.set_ylabel("z",fontsize=28)
	ax.tick_params(axis="both",labelsize=28)
end
function test_dstep()
	eps = 1.e-6
	s = [1., 4., 0.]
	u = rand(3)
	u1px = u .+ eps*[1,0,0]
	u1mx = u .- eps*[1,0,0]
	u1py = u .+ eps*[0,1,0]
	u1my = u .- eps*[0,1,0]
	u1pz = u .+ eps*[0,0,1]
	u1mz = u .- eps*[0,0,1]

	du_dx_fd = (next(u1px,s) - next(u1mx,s))/(2*eps)
	du_dy_fd = (next(u1py,s) - next(u1my,s))/(2*eps)
	du_dz_fd = (next(u1pz,s) - next(u1mz,s))/(2*eps)

	du = dstep(u, s)
	du_dx = du[:,1]
	du_dy = du[:,2]
	du_dz = du[:,3]

	@show abs.(du_dx_fd .- du_dx)
	@show abs.(du_dy_fd .- du_dy)
	@show abs.(du_dz_fd .- du_dz)
	@test du_dx_fd ≈ du_dx atol=1.e-8
	@test du_dy_fd ≈ du_dy atol=1.e-8
	@test du_dz_fd ≈ du_dz atol=1.e-8
end
function test_d2step()
# d2step has to be defined in the main script as below, as opposed to writing a function in solenoid.jl, because we don't want to do source-to-source AD every time the second derivative is needed.
	eps = 1.e-6
	s = [1.0, 4.0,0.1]
	u = rand(3)
	u1px = u .+ eps*[1,0,0]
	u1mx = u .- eps*[1,0,0]
	u1py = u .+ eps*[0,1,0]
	u1my = u .- eps*[0,1,0]
	u1pz = u .+ eps*[0,0,1]
	u1mz = u .- eps*[0,0,1]

	ddu_dx_fd = (dstep(u1px,s) - dstep(u1mx,s))/(2*eps)
	ddu_dy_fd = (dstep(u1py,s) - dstep(u1my,s))/(2*eps)
	ddu_dz_fd = (dstep(u1pz,s) - dstep(u1mz,s))/(2*eps)

	ddu(x) = jacobian(x -> dstep(x,s), x)
	ddu_ad = ddu(u)[1]
	ddu_dx_ad = reshape(ddu_ad[:,1],3,3)
	ddu_dy_ad = reshape(ddu_ad[:,2],3,3)
	ddu_dz_ad = reshape(ddu_ad[:,3],3,3)

	@show ddu_dx_ad .- ddu_dx_fd
	@show ddu_dy_ad .- ddu_dy_fd
	@show ddu_dz_ad .- ddu_dz_fd

	@test ddu_dx_ad ≈ ddu_dx_fd atol=1.e-8
	@test ddu_dy_ad ≈ ddu_dy_fd atol=1.e-8
	@test ddu_dz_ad ≈ ddu_dz_fd atol=1.e-8

end
function test_pert()
	eps = 1.e-6
	s = [1., 4., 0.]
	u = rand(3)
	for i = 1:3
		eps_arr = zeros(3)
		eps_arr[i] = eps
		sp = s .+ eps_arr
		sm = s .- eps_arr

		pp_fd = (next(u, sp) - next(u, sm))/(2*eps)
		pp = pert(u, s, i)
		@show i
		@show pp_fd, pp
		@test pp_fd ≈ pp atol=1.e-8
	end
end

max_stars_repo_name: nishaChandramoorthy/linearResponse

max_issues_repo_name: nishaChandramoorthy/linearResponse

max_forks_repo_name: nishaChandramoorthy/linearResponse


text:
using Documenter, FWF

makedocs(modules=[FWF],
        doctest=true)

deploydocs(deps   = Deps.pip("mkdocs", "python-markdown-math"),
    repo = "github.com/GITHUBNAME/GITHUBREPO.git",
    julia  = ".6",
    osname = "linux")

ext: jl


text:
export ft_setbaudrate

function ft_setbaudrate(ft_handle::Culong, baud::Integer)
  ft_status = ccall((:FT_SetBaudRate, d2xx),
                     Cuint,
                     (Culong,Cuint),
                     ft_handle,baud)
  checkstatus(ft_status)
  return nothing
end


text:
# v0.11 deprecations
@deprecate poisson poisson_loss false
@deprecate hinge hinge_loss false
@deprecate squared_hinge squared_hinge_loss false
@deprecate normalise(x) normalise(x, dims=1) false

@deprecate binarycrossentropy(ŷ, y) Losses.binarycrossentropy(ŷ, y, agg=identity) false
@deprecate logitbinarycrossentropy(ŷ, y) Losses.logitbinarycrossentropy(ŷ, y, agg=identity) false

function Broadcast.broadcasted(::typeof(binarycrossentropy), ŷ, y)
    @warn "binarycrossentropy.(ŷ, y) is deprecated, use Losses.binarycrossentropy(ŷ, y, agg=identity) instead"
    Losses.binarycrossentropy(ŷ, y, agg=identity)
end

function Broadcast.broadcasted(::typeof(logitbinarycrossentropy), ŷ, y)
    @warn "logitbinarycrossentropy.(ŷ, y) is deprecated, use Losses.logitbinarycrossentropy(ŷ, y, agg=identity) instead"
    Losses.logitbinarycrossentropy(ŷ, y, agg=identity)
end


text:
using Soss

contaminatedNormal = @model begin
    p ~ Uniform()
    x ~ Normal()
    ε ~ Normal(σ=10)
    return p<0.1 ? x : ε
end

m = @model prior begin
    x ~ prior()
    y ~ Binomial(n=10, logitp=x)
    return y
end

y = rand(m(prior=contaminatedNormal))

using SampleChainsDynamicHMC

sample(DynamicHMCChain, m(prior=contaminatedNormal) | (;y=4))


text:
using CImGui
using ImGuiMakieBackend
using GLMakie

x = rand(10)
y = rand(10)
colors = rand(10)
scene = scatter(x, y, color = colors)

set_makie_renderloop()


text:
# SPD-License-Identifier: MIT

using ArgParse
using JLD2
using MarkovModels
using TOML

function parse_commandline()
    s = ArgParseSettings()
    @add_arg_table s begin
        "topology"
            required = true
            help = "hmm topology in TOML format"
        "units"
            required = true
            help = "list of units with their categories"
        "hmms"
            required = true
            help = "output hmms in BSON format"
    end
    s.description = """
    Build a set of HMM. The topology file should be formatted as:

    states = [

    \ua0\ua0{ id = 1, initweight = 1.0, finalweight = 0.0 },

    \ua0\ua0...

    ]

    links = [

    \ua0\ua0{ src = 1, dest = 2, weight = 0.5 },

    \ua0\ua0...

    ]

    tiestates = false | true

    """
    parse_args(s)
end

function loadunits(file)
    units, categories = [], []
    open(file, "r") do f
        for line in eachline(f)
            tokens = split(line)
            push!(units, tokens[1])
            push!(categories, tokens[2:end])
        end
    end
    units, categories
end

function get_unit_topo(topo, category)
    i = 1
    while i <= length(category) && category[i] ∈ keys(topo)
        topo = topo[category[i]]
        i += 1
    end
    topo
end

function makehmm!(pdfid_mapping, unit, topo, pdfid)
    SF = LogSemifield{Float32}
    fsm = VectorFSM{SF}()

    states = Dict()
    for (i, state) in enumerate(topo["states"])
        initweight = SF(log(state["initweight"]))
        finalweight = SF(log(state["finalweight"]))
        s = addstate!(fsm, i; initweight, finalweight)
        states[i] = s
        pdfid_mapping[(unit, i)] = pdfid
        pdfid += 1
    end

    for arc in topo["arcs"]
        addarc!(fsm, states[arc["src"]], states[arc["dest"]],
                SF(log(arc["weight"])))
    end

    fsm |> renormalize, pdfid
end

function main(args)
    topo = TOML.parsefile(args["topology"])

    units, categories = loadunits(args["units"])
    hmms = Dict()
    pdfid_mapping = Dict()
    pdfid = 1
    for (unit, category) in zip(units, categories)
        unit_topo = get_unit_topo(topo, category)
        hmms[unit], pdfid = makehmm!(pdfid_mapping, unit, topo, pdfid)
    end

    data = Dict(
        "units" => hmms,
        "pdfid_mapping" => pdfid_mapping
    )
    save(args["hmms"], data)
end

args = parse_commandline()
main(args)

max_stars_repo_name: MartinKocour/SpeechLab

max_issues_repo_name: MartinKocour/SpeechLab

max_forks_repo_name: MartinKocour/SpeechLab


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

# tests the output of the embedding example is correct
using Test

if Sys.iswindows()
    # libjulia needs to be in the same directory as the embedding executable or in path
    ENV["PATH"] = string(Sys.BINDIR, ";", ENV["PATH"])
end

@test length(ARGS) == 1

@testset "embedding example" begin
    out = Pipe()
    err = Pipe()
    embedded_cmd_path = abspath(ARGS[1])
    p = cd(@__DIR__) do
        run(pipeline(Cmd([embedded_cmd_path]), stdin=devnull, stdout=out, stderr=err), wait=false)
    end
    close(out.in)
    close(err.in)
    out_task = @async readlines(out)
    err = read(err, String)
    @test err == "MethodError: no method matching this_function_has_no_methods()\n"
    @test success(p)
    lines = fetch(out_task)
    @test length(lines) == 10
    @test parse(Float64, lines[1]) ≈ sqrt(2)
    @test lines[8] == "called bar"
    @test lines[9] == "calling new bar"
    @test lines[10] == "      From worker 2:\tTaking over the world..."
end

max_stars_repo_name: TimoLarson/julia

max_issues_repo_name: AsForMe/julia

max_forks_repo_name: AsForMe/julia


text:
using Base.Test

export @customtest,
    @customtest_throws

#=
customtest and customtest_throws are wrappers around the normal
@test and @test_throws macros. They ensure that the test execution
is independent of the test order by resetting the interal package state.
Currently, the internal package state consists only of the next free
symbol.
=#

macro customtest(test)
    quote
    resetfreesymbols()
    @test $(esc(test))
    resetfreesymbols()
    end
end

macro customtest_throws(exception, test)
    quote
    resetfreesymbols()
    @test_throws $(esc(exception)) $(esc(test))
    resetfreesymbols()
    end
end


text:
"
MutationOperator represents the method that mutates the individuals selected to
create a new generation of individuals.

# Fields
- `_method::Function`: method used for mutation.
- `_probability::Float32`: probability of mutation.
- `_varArgs::Array{Any}`: arguments necessary for the mutation method.
"
struct MutationOperator
    _method::Function
    _probability::Float32
    _varArgs::Array{Any}
end # struct

"""
    getMethod(mutationOp::MutationOperator)::Function

Returns the method for mutate the offspring.
"""
getMethod(mutationOp::MutationOperator)::Function = mutationOp._method
# function

"""
    getMutateProbability(mutationOp::MutationOperator)::Float32

Obtains the probability for a set of individuals to be mutated.
"""
getMutateProbability(mutationOp::MutationOperator)::Float32 = mutationOp._probability
# function

"""
    getFunctionArgs(mutationOp::MutationOperator)::Array

Obtains the aditional arguments associated to mutation method.
"""
getFunctionArgs(mutationOp::MutationOperator)::Array = mutationOp._varArgs
#function

"""
    getDefaultMutationOp(individualType::DataType)::Tuple{Function, Array}

Returns a default mutation method for each registered individual type. For
`BinaryGenotype`, it is [`standardMutation`](@ref) with a probability of 0.1;
for `IntegerGenotype`, it is [`uniformMutation`](@ref) of one gene within a range of 0:10;
for both `CGPGenotype` and `STGPGenotype` it is [`pointMutation`](@ref) with a
probability of 0.2; for `GEPGenotype` it is [`geneMutation`](@ref).
"""
function getDefaultMutationOp(individualType::DataType)::Tuple{Function, Array}
    method = nothing
    varArgs = []
    if individualType <: BinaryGenotype
        method = standardMutation
        varArgs = [0.1]
    elseif individualType <: IntegerGenotype
        method = uniformMutation
        varArgs = [0, 10]
    elseif individualType <: CGPGenotype || individualType <: STGPGenotype
        method = pointMutation
        varArgs = [0.2]
    elseif individualType <: GEPGenotype
        method = geneMutation
    end

    return method, varArgs
end # function

"""
    mutate_(mutationOp::MutationOperator, offspring::Array{Individual},
            experimentInfo::ExperimentInfo)::Array{Individual}

Performs a mutation on a set of individuals.
"""
function mutate_(mutationOp::MutationOperator, offspring::Array{Individual},
                 experimentInfo::ExperimentInfo)::Array{Individual}

    nSelected = length(offspring)
    offspringRep = getGenotype(offspring)
    rng = getRNG(experimentInfo)
    mutateProbability = getMutateProbability(mutationOp)
    method = getMethod(mutationOp)
    functionArgs = getFunctionArgs(mutationOp)

    if experimentInfo._individualType <: GAGenotype
        for i in eachindex(offspring)

            random = rand(rng)

            if random < mutateProbability
                offspringRep[i] = method(offspringRep[i], rng,  functionArgs...)
                offspring[i] = Individual(offspringRep[i])
            end
        end
    else
        for i in eachindex(offspring)

            random = rand(rng)

            if random < mutateProbability
                offspringRep[i] = method(offspringRep[i], experimentInfo._GPExperimentInfo, rng, functionArgs...)
                offspring[i] = Individual(offspringRep[i])
            end
        end
    end

    return offspring
end # function

max_stars_repo_name: SergioGmezM/GeneticJulia

max_issues_repo_name: SergioGmezM/GeneticJulia

max_forks_repo_name: SergioGmezM/GeneticJulia


text:
# ---
# title: 1222. Queens That Can Attack the King
# id: problem1222
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array
# link: <https://leetcode.com/problems/queens-that-can-attack-the-king/description/>
# hidden: true
# ---
#
# On an **8x8** chessboard, there can be multiple Black Queens and one White
# King.
#
# Given an array of integer coordinates `queens` that represents the positions
# of the Black Queens, and a pair of coordinates `king` that represent the
# position of the White King, return the coordinates of all the queens (in any
# order) that can attack the King.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2019/10/01/untitled-diagram.jpg)
#
#
#
#     Input: queens = [[0,1],[1,0],[4,0],[0,4],[3,3],[2,4]], king = [0,0]
#     Output: [[0,1],[1,0],[3,3]]
#     Explanation:
#     The queen at [0,1] can attack the king cause they're in the same row.
#     The queen at [1,0] can attack the king cause they're in the same column.
#     The queen at [3,3] can attack the king cause they're in the same diagnal.
#     The queen at [0,4] can't attack the king cause it's blocked by the queen at [0,1].
#     The queen at [4,0] can't attack the king cause it's blocked by the queen at [1,0].
#     The queen at [2,4] can't attack the king cause it's not in the same row/column/diagnal as the king.
#
#
# **Example 2:**
#
# **![](https://assets.leetcode.com/uploads/2019/10/01/untitled-diagram-1.jpg)**
#
#
#
#     Input: queens = [[0,0],[1,1],[2,2],[3,4],[3,5],[4,4],[4,5]], king = [3,3]
#     Output: [[2,2],[3,4],[4,4]]
#
#
# **Example 3:**
#
# **![](https://assets.leetcode.com/uploads/2019/10/01/untitled-diagram-2.jpg)**
#
#
#
#     Input: queens = [[5,6],[7,7],[2,1],[0,7],[1,6],[5,1],[3,7],[0,3],[4,0],[1,2],[6,3],[5,0],[0,4],[2,2],[1,1],[6,4],[5,4],[0,0],[2,6],[4,5],[5,2],[1,4],[7,5],[2,3],[0,5],[4,2],[1,0],[2,7],[0,1],[4,6],[6,1],[0,6],[4,3],[1,7]], king = [3,4]
#     Output: [[2,3],[1,4],[1,6],[3,7],[4,3],[5,4],[4,5]]
#
#
#
#
# **Constraints:**
#
#   * `1 <= queens.length <= 63`
#   * `queens[0].length == 2`
#   * `0 <= queens[i][j] < 8`
#   * `king.length == 2`
#   * `0 <= king[0], king[1] < 8`
#   * At most one piece is allowed in a cell.
#
#
## @lc code=start
using LeetCode

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


text:
push!(LOAD_PATH, "../src")

using HDF5
using StaticArrays
using Statistics
using LinearAlgebra
using .Threads
using Serialization
using Random
using Healpix
using Parameters

#include("GadgetReader.jl")
#include("OctTree.jl")
using GadgetReader
using OctTree

using AstroChemistry

const N = 3

function vec2svec(vec::Matrix{T}) where {T}
    svec = [SVector{3,T}(vec[:,i]) for i in 1:size(vec,2)]
end
function mat2smat(mat::Array{T,3}) where {T}
    smat = [SMatrix{3,3,T}(mat[:,:,i]) for i in 1:size(mat,3)]
end

const ANGLE = 0.7
const ShieldingLength = 1e10

const XH = 0.71
const BOLTZMANN=1.3806e-16
const PROTONMASS=1.6726e-24
const GRAVCON=6.67e-8
const UnitMass_in_g = 1.989e43
const UnitLength_in_cm    = 3.085678e21
const UnitTime_in_s = 3.08568e+16

const UnitDensity_in_cgs = UnitMass_in_g / UnitLength_in_cm^3
const UnitDensity_in_pccm = UnitDensity_in_cgs/PROTONMASS
const Year_in_s = 31556926.

const fac_col = (UnitMass_in_g/UnitLength_in_cm^2)/PROTONMASS

const facNHtoAv = 5.35e-22
const gamma_eff = 3.51

function solve_chem_all_particles(file_path, Zp)

    T = Float64

    all_species::Vector{String} =
    [
    "H", "H-", "H2", "H+", "H2+", "H3+", "e-", "He", "He+", "HeH+",
    "C", "C+", "CO", "HCO+",
    "O", "O+", "OH", "OH+", "H2O+", "H3O+", "H2O",
    "O2", "CO+", "O2+",
    "CH2", "CH2+", "CH", "CH+", "CH3+",
    "Si+", "Si"
    #"SiO", "SiO+", "SiOH+",
    #"S", "S+",
    ]
    N_spec = length(all_species)
    N_neq = 0

    net, dict = initialize_chemistry_network(all_species)
    @unpack iH2, iCO, iC, fac_H, fac_C, fac_O, charge = net
    abtot = AbundTotal() * Zp
    @unpack abC, abO, abSi = abtot

    fname = file_path * "/snapshot.hdf5"
    Npart, pos, vel, rho, u, mass, hsml, id_gas,
    boxsize, time = read_snap(fname);
    println("boxsize = ", boxsize)

    boxsizes = SVector(boxsize, boxsize, boxsize)
    center = 0.5 .* boxsizes
    topnode_length = boxsizes

    X = vec2svec(pos);

    ξ = 1.3e-16  #H2
    IUV = 1.0

    mu = 2.3 #mean molecular weight
    nH = rho .* (XH * UnitDensity_in_cgs / PROTONMASS)
    temp = u ./ (1.5 * BOLTZMANN /mu /PROTONMASS / 1e10)

    mass_H2 = zeros(Npart)
    mass_CO = zeros(Npart)

    #NH = NH2 = NCO = NC = 0.0 #don't declear here!! otherwise they won't be thread-safe

    abund_all = [zeros(N_spec) for _ in 1:Npart]

    for i in 1:Npart
        xneq = SVector{N_neq,T}()
        #make sure the abC is consistent between simulations & postprocessing
        abund_all[i][dict["C+"]] = abC
        #initial guess for H2 & H+ (useful for calcultating steady state H2)
        #abund_all[i][dict["H2"]] = 0.49
        init_abund(abund_all[i], xneq, abtot, net)
    end

    NH_eff = zeros(Npart)
    NH2_eff = zeros(Npart)
    NCO_eff = zeros(Npart)
    ga_out = Vector{TreeGather{Float64}}(undef,Npart)
    tree_out = nothing
    status = zeros(Int, Npart)

    Nstep = 3
    time_max = 1e9 #unit = yr
    dt = time_max / Nstep
    NCpix = zeros(NPIX)

    for j in 1:Nstep
        @. mass_H2 = mass * getindex(abund_all,iH2) * XH * 2.0
        @. mass_CO = mass * getindex(abund_all,iCO) * XH * 28.0

        println("iteration ", j)

        println("buildtree...")
        @time tree = buildtree(X, hsml, mass, mass_H2, mass_CO, center, topnode_length);

        println("loop over particles and solve the chemistry network...")
        ichem = findall((nH.>0.1).&(temp.<5e3))
        println(length(ichem), " cold and dense particles found...")
        @time @threads for i in ichem #better load balance
            if i%1000 == 0
                print("i=", i, " ")
            end
            ga = TreeGather{T}()
            treewalk(ga,X[i],tree,ANGLE,ShieldingLength,boxsizes)
            ga_out[i] = ga
            NH = ga.column_all .* (fac_col*XH)
            NH2 = ga.column_H2 .* (fac_col/2)
            NCO = ga.column_CO .* (fac_col/28)
            NC = 0.0
            facNHtoAv_eff = gamma_eff*facNHtoAv*Zp
            NH_eff[i]  = -log(mean(exp.(-NH .*facNHtoAv_eff))) / facNHtoAv_eff
            NH2_eff[i] = -log(mean(exp.(-NH2.*facNHtoAv_eff))) / facNHtoAv_eff
            NCO_eff[i] = -log(mean(exp.(-NCO.*facNHtoAv_eff))) / facNHtoAv_eff

            xneq = SVector{N_neq,T}()
            if nH[i] > 1e5 nH[i] = 1e5 end
            par = Par{NPIX,N_neq,T}(nH[i], temp[i], ξ, IUV, Zp, Zp,
                SVector{NPIX,T}(NH),
                SVector{NPIX,T}(NH2),
                SVector{NPIX,T}(NCO),
                SVector{NPIX,T}(NCpix), xneq)
            #@show NH,NH2,NCO,NCpix,xneq
            #dt = 1e9 / (Zp * nH[i]) #in years
            retcode, _  = solve_equilibrium_abundances(abund_all[i], dt, par, abtot, net)
            if retcode != :Success
                status[i] = 1  #1 means failed
                @show retcode
                @show i, nH[i], temp[i], ξ, IUV, Zp
                @show NH_eff[i], NH2_eff[i], NCO_eff[i], xneq
                #@show NH, NH2, NCO
            end
        end
        println("chemsitry done!")
        @show sum(status)

        xCO=getindex.(abund_all, dict["CO"])
        println("CO mass = ",  1e10*sum(mass.*xCO) * XH * 28 )
        xH2=getindex.(abund_all, dict["H2"])
        println("H2 mass = ",  1e10*sum(mass.*xH2) * XH * 2  )

        #this may seem redundant as it'll be computed in the next iteration
        #however, treewalk is so fast (compared to chemistry) that it doesn't hurt
        #and we get the updated column densities for little overhead
        println("get column densities for the outputs...")
        @. mass_H2 = mass * getindex(abund_all,iH2) * XH * 2.0
        @. mass_CO = mass * getindex(abund_all,iCO) * XH * 28.0
        println("buildtree...")
        @time tree = buildtree(X, hsml, mass, mass_H2, mass_CO, center, topnode_length);
        println("loop over particles and get column densities...")
        @time @threads for i in eachindex(NH_eff)
            ga = TreeGather{T}()
            treewalk(ga,X[i],tree,ANGLE,ShieldingLength,boxsizes)
            ga_out[i] = ga
            NH = ga.column_all .* (fac_col*XH)
            NH2 = ga.column_H2 .* (fac_col/2)
            NCO = ga.column_CO .* (fac_col/28)
            facNHtoAv_eff = gamma_eff*facNHtoAv*Zp
            NH_eff[i]  = -log(mean(exp.(-NH .*facNHtoAv_eff))) / facNHtoAv_eff
            NH2_eff[i] = -log(mean(exp.(-NH2.*facNHtoAv_eff))) / facNHtoAv_eff
            NCO_eff[i] = -log(mean(exp.(-NCO.*facNHtoAv_eff))) / facNHtoAv_eff
        end
        println("done!")
        tree_out = tree

        # write to file
        abund_all_arr = zeros(N_spec,Npart)
        for i in eachindex(abund_all)
            abund_all_arr[:,i] = abund_all[i]
        end
        #T = Float64
        fnamebase = "/data-chem-"
        if net.grRec == true
            fnamebase *= "GrRec-"
        end
        fname = file_path * fnamebase * string(j) *".hdf5"
        fid=h5open(fname,"w")
        grp_head = create_group(fid,"Header");
        attributes(fid["Header"])["all_species"] = all_species
        attributes(fid["Header"])["time"] = time
        grp_part = create_group(fid,"Chemistry");
        h5write(fname, "Gas/Position"        , pos)
        h5write(fname, "Gas/Density"         , nH)
        h5write(fname, "Gas/Temperature"     , temp)
        h5write(fname, "Gas/ID"              , id_gas)
        h5write(fname, "Chemistry/Abundances"     , abund_all_arr)
        h5write(fname, "Chemistry/NH_eff"         , NH_eff)
        h5write(fname, "Chemistry/NH2_eff"        , NH2_eff)
        h5write(fname, "Chemistry/NCO_eff"        , NCO_eff)
        h5write(fname, "Chemistry/Status"         , status)
        close(fid)
    end #iteration loop
    return pos, nH, temp, id_gas, abund_all, ga_out, X, NH_eff, NH2_eff, NCO_eff, tree_out
end

file_path = "./"
Zp = 1.0
pos, nH, temp, id_gas, abund_all, ga_out, X, NH_eff, NH2_eff, NCO_eff, tree_out = solve_chem_all_particles(file_path, Zp)


text:
# Some quick and dirty broadcasting AD benchmarks, 25 April 2020

julia> using Zygote, Tracker, ReverseDiff, Tullio

julia> x1 = randn(10); f1(x) = 1+tanh(x); # A simple function which isn't hard-coded

julia> Tracker.gradient(x -> sum(f1.(x)), x1)[1]'
Tracked 1×10 LinearAlgebra.Adjoint{Float64,Array{Float64,1}}:
 0.333674  0.174402  0.786548  0.222385  …  0.143003  0.361015  0.910887  0.940078

julia> Zygote.gradient(x -> sum(f1.(x)), x1)[1]'
1×10 LinearAlgebra.Adjoint{Float64,Array{Float64,1}}:
 0.333674  0.174402  0.786548  0.222385  …  0.143003  0.361015  0.910887  0.940078

julia> ReverseDiff.gradient(x -> sum(f1.(x)), (x1,))[1]'
1×10 LinearAlgebra.Adjoint{Float64,Array{Float64,1}}:
 0.333674  0.174402  0.786548  0.222385  …  0.143003  0.361015  0.910887  0.940078

julia> Tracker.gradient(x -> (@tullio s = 1+tanh(x[i])), x1)[1]'
Tracked 1×10 LinearAlgebra.Adjoint{Float64,Array{Float64,1}}:
 0.333674  0.174402  0.786548  0.222385  …  0.143003  0.361015  0.910887  0.940078

julia> Zygote.gradient(x -> sum(@tullio s[i] := 1+tanh(x[i])), x1)[1]'
1×10 LinearAlgebra.Adjoint{Float64,Array{Float64,1}}:
 0.333674  0.174402  0.786548  0.222385  …  0.143003  0.361015  0.910887  0.940078

# Time that on a bigger array:
# right now Zygote is very slow, and ReverseDiff worse.

julia> x2 = randn(1000,1000);

julia> @btime sum(f1.($x2));
  25.181 ms (2 allocations: 7.63 MiB)

julia> @btime sum(f1, $x2);
  24.414 ms (0 allocations: 0 bytes)

julia> @btime sum(Broadcast.broadcasted(f1, $x2)) # on Julia 1.5
  29.678 ms (0 allocations: 0 bytes)

julia> @btime @tullio s = f1($x2[i,j])
  13.195 ms (41 allocations: 2.66 KiB)

julia> @btime Tracker.gradient(x -> sum(f1.(x)), $x2);
  58.780 ms (44 allocations: 53.41 MiB)

julia> @btime Zygote.gradient(x -> sum(f1.(x)), $x2);
  119.821 ms (3000046 allocations: 122.07 MiB)

julia> @btime ReverseDiff.gradient(x -> sum(f1.(x)), ($x2,));
  941.047 ms (15000027 allocations: 567.06 MiB)

julia> @btime Tracker.gradient(x -> (@tullio s = 1+tanh(x[i,j])), $x2);
  33.391 ms (88 allocations: 30.52 MiB)

julia> @btime Zygote.gradient(x -> sum(@tullio s[i] := 1+tanh(x[i,j])), $x2);
  29.928 ms (84 allocations: 7.64 MiB)

julia> Tullio.@printgrad 1+tanh(x) x
δx = 1 - tanh(x) ^ 2

# Compare to tanh, which is a special case of Zygote's:

julia> @btime Zygote.gradient(x -> sum(tanh.(x)), $x2);
  25.633 ms (20 allocations: 15.26 MiB)

# Fancier ways to use @tullio:

julia> using LoopVectorization

julia> @btime @tullio s = f1($x2[i,j])
  3.194 ms (39 allocations: 2.59 KiB)

julia> @btime Zygote.gradient(x -> sum(@tullio s[i] := 1+tanh(x[i,j])), $x2);
  10.663 ms (81 allocations: 7.64 MiB)

julia> @btime Zygote.gradient(x -> (@tullio s := 1+tanh(x[i,j])), $x2); # skip making s[i]
  9.406 ms (79 allocations: 7.63 MiB)

julia> using ForwardDiff

julia> @btime Zygote.gradient(x -> sum(@tullio s[i] := 1+tanh(x[i,j]) grad=Dual), $x2);
  10.512 ms (81 allocations: 7.64 MiB)

# Another problem:

julia> @btime sum($x2 .+ $x2' ./ 2);
  2.968 ms (2 allocations: 7.63 MiB)

julia> @btime Tracker.gradient(x -> sum(x .+ x' ./ 2), $x2);
  25.699 ms (210 allocations: 68.67 MiB)

julia> @btime Zygote.gradient(x -> sum(x .+ x' ./ 2), $x2);
  9.006 ms (15 allocations: 38.15 MiB)

julia> @btime ReverseDiff.gradient(x -> sum(x .+ x' ./ 2), $x2);
  1.219 s (19000027 allocations: 719.65 MiB)

julia> @btime Zygote.gradient(x -> (@tullio s := x[i,j] + x[j,i]/2), $x2);
  1.756 ms (164 allocations: 7.64 MiB)

# And without @avx magic...

julia> @btime Zygote.gradient(x -> (@tullio s := x[i,j] + x[j,i]/2  avx=false), $x2);
  2.368 ms (165 allocations: 7.64 MiB)

julia> @btime Zygote.gradient(x -> (@tullio s := x[i,j] + x[j,i]/2  avx=false threads=false), $x2);
  5.146 ms (25 allocations: 7.63 MiB)

julia> @btime Tracker.gradient(x -> (@tullio s := x[i,j] + x[j,i]/2  avx=false), $x2);
  7.169 ms (172 allocations: 30.53 MiB)

julia> @btime ReverseDiff.gradient(x -> (@tullio s := x[i,j] + x[j,i]/2  avx=false), $x2);
  7.955 ms (164 allocations: 30.53 MiB)

# ReverseDiff should soon have forward-mode broadcasting from this package:

julia> using DistributionsAD

julia> @btime ReverseDiff.gradient(x -> sum(f1.(x)), ($x2,));
  35.920 ms (34 allocations: 53.41 MiB)

julia> @btime ReverseDiff.gradient(x -> sum(x .+ x' ./ 2), $x2);
  1.365 s (19000087 allocations: 727.28 MiB)

julia> @btime ReverseDiff.gradient(x -> sum(x .+ x'), $x2); # without the ./2 bit!
  48.920 ms (1000030 allocations: 99.18 MiB)

# Yota only supports broadcasting of primitives, i.e. things with an explicit scalar rule

julia> using Yota

julia> @diffrule f1(x)  x  dy*(1-tanh(x)^2)

julia> @btime Yota.grad(x -> sum(f1.(x)), $x2);
  77.449 ms (86 allocations: 45.78 MiB)

julia> @btime Yota.grad(x -> sum(x .+ x' ./ 2), $x2);
  8.423 ms (89 allocations: 38.15 MiB)


text:
####
#### [x|y|x2|y2]title[!]
####

"""
    $SIGNATURES

Internal function to set the title of axes (`el==:title`) or axis objects (`el==:xtitle`...).
"""
function _title(el::Symbol, text::String=""; axes=nothing, opts...)
    axes = check_axes(axes)
    # retrieve on what object to apply the title either the current axes or a sub-axis
    obj = (el == :axis) ? axes : getfield(axes, el)
    # create a new title object with the text and apply properties
    obj.title = Title(text=text)
    set_properties!(obj.title; defer_preview=true, opts...)
    return PreviewFigure(gcf())
end

# Generate xtitle, and associated for each axis
for axs ∈ ("", "x", "y", "x2", "y2")
    f   = Symbol(axs * "title")
    f2  = Symbol(axs * "label")
    ex = quote
        # mutate
        $f(t::String=""; o...) = _title(Symbol($axs * "axis"), t; axes=gca(), o...)
        # synonyms xlabel, ylabel, etc
        !isempty($axs) && ($f2 = $f)
    end
    eval(ex)
end

####
#### [x|y|x2|y2]ticks
####

function _ticks(axis_sym::Symbol, loc::Vector{Float64}=Float64[],
                lab::Vector{String}=String[]; axes=nothing, opts...)
    axes = check_axes(axes)
    # retrieve the appropriate axis
    axis = getfield(axes, axis_sym)
    # create a new ticks object
    axis.ticks = Ticks()
    # if locs are empty, just pass options and return
    if isempty(loc)
        isempty(lab) || throw(ArgumentError("Cannot pass ticks labels without specifying " *
                                            "ticks locations"))
        set_properties!(axis.ticks; defer_preview=true, opts...)
        return PreviewFigure(gcf())
    end
    axis.ticks.places = loc
    # process labels if any are passed
    if !isempty(lab)
        length(lab) == length(loc) || throw(ArgumentError("Ticks locations and labels must " *
                                                          "have the same length."))
        axis.ticks.labels = TicksLabels(names=lab)
    end
    # set remaining properties and return
    set_properties!(axis.ticks; defer_preview=true, opts...)
    return PreviewFigure(gcf())
end

# Generate xticks!, xticks, and associated for each axis
for axs ∈ ("x", "y", "x2", "y2")
    f  = Symbol(axs * "ticks")
    a  = Symbol(axs * "axis")
    ex = quote
        # xticks!([1, 2], ["a", "b"]; o...)
        $f(loc::AVR, lab=String[]; o...) = _ticks(Symbol($axs * "axis"), fl(loc), lab; o...)
        # xticks!("off"; o...)
        function $f(s::String=""; o...)
            isempty(s) && return _ticks(Symbol($axs * "axis"); o...)
            ax = getfield(gca(), Symbol($axs * "axis"))
            s_lc = lowercase(s)
            if s_lc == "off"
                reset!(ax.ticks)
                ax.ticks.off = true
                ax.ticks.labels.off = true
            elseif s_lc == "on"
                ax.ticks.off = false
                ax.ticks.labels.off = false
                set_properties!(ax.ticks; defer_preview=true, o...)
            else
                throw(ArgumentError("Unrecognised shorthand: $s"))
            end
            return PreviewFigure(gcf())
        end
    end
    eval(ex)
end


text:
using Logging, LoggingExtras, Dates

const EMAIL = # Bot email
const API_KEY = # Bot api key
const ZULIP_EP = "https://julialang.zulipchat.com/"

const SODB = "db/sodb.sqlite"
const LOGS_DIR = "<path to logs>"

const date_format = "yyyy-mm-dd HH:MM:SS"

timestamp_logger(logger) = TransformerLogger(logger) do log
    merge(log, (; message = "[$(Dates.format(now(), date_format))] $(log.message)"))
end

ConsoleLogger(stdout, show_limited = false) |> timestamp_logger |> global_logger


text:
# Autogenerated wrapper script for CSDP_jll for i686-w64-mingw32
export libcsdp

using OpenBLAS32_jll
using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("CSDP")
JLLWrappers.@declare_library_product(libcsdp, "libcsdp.dll")
function __init__()
    JLLWrappers.@generate_init_header(OpenBLAS32_jll, CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libcsdp,
        "bin\\libcsdp.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
"""
    fit_adc(bvalues, signals)

Fit the apparent diffusion coefficient (ADC) using a polynomial logfit of the normalized signals
`signals` against the b-values `bvalues`.
"""
function fit_adc(bvalues, signals)
    # Fit log of signal
    p = fit(bvalues, log.(abs.(signals)))

    # ADC is the negative linear term
    -p[1]
end


text:
using Primes: nextprime

################################################################
function fed_tau_by_process_time(T::Float64, M::Int64, tau_max::Float64,
                            reordering::Bool)
    ## All cycles have the same fraction of the stopping time
    fed_tau_by_cycle_time(T/M, tau_max, reordering)
end

################################################################
function fed_tau_by_cycle_time(t::Float64, tau_max::Float64,
                               reordering::Bool)

    n = 0          ## Number of time steps
    scale = 0.0    ## Ratio of t we search to maximal t

    ## Compute necessary number of time steps
    # n = ceil(Int64, sqrt(3.0*t/tau_max+0.25)-0.5-1.0e-8)

    n = floor(Int64, ceil(sqrt(3.0 * t / tau_max + 0.25) - 0.5 - 1.0e-8)+ 0.5)

    scale = 3.0*t/(tau_max*Float64(n*(n+1)))

    ## Call internal FED time step creation routine
    fed_tau_internal(n,scale,tau_max,reordering)
end

################################################################
function fed_tau_internal(n::Int64,  scale::Float64, tau_max::Float64,
                          reordering::Bool)

    if (n <= 0)
        return 0, nothing
    end

    #### Allocate memory for the time step size
    tau = zeros(n)

    if (reordering)
        tauh = zeros(n)
    end

    ## Compute time saver
    c = 1.0 / (4.0 * n + 2.0)
    d = scale * tau_max / 2.0

    ## Set up originally ordered tau vector
    for k in 1:n
        h = cos(pi * (2.0 * k - 1.0) * c)

        if reordering
            tauh[k] = d / (h * h);
        else
            tau[k] = d / (h * h);
        end
    end

    ## Permute list of time steps according to chosen reordering function
    if reordering
        ## Choose kappa cycle with k = n/2
        ## This is a heuristic. We can use Leja ordering instead!!
        kappa = n ÷ 2

        ## Get modulus for permutation
        prime = nextprime(n + 1)

        ## Perform permutation
        k = 0
        l = 0
        while l < n
            index = ((k+1)*kappa) % prime - 1
            while (index >= n)
                k+=1
                index = ((k+1)*kappa) % prime - 1
            end

            tau[l+1] = tauh[index+1]

            k += 1
            l += 1
        end
    end

    return n, tau
end


text:
@testset "Chimera-like instance has the correct low energy spectrum" begin
    m = 3
    n = 4
    t = 3

    β = 1.

    L = n * m * t
    num_states = 22

    # energies
    exact_energies =
    [
        -16.4, -16.4, -16.4, -16.4, -16.1, -16.1, -16.1, -16.1, -15.9,
        -15.9, -15.9, -15.9, -15.9, -15.9, -15.6, -15.6, -15.6, -15.6,
        -15.6, -15.6, -15.4, -15.4
    ]

    # degenerate fg solutions
    exact_states =
    [   # E =-16.4
        [
        [1, 4, 5, 1, 2, 2, 1, 1, 1, 4, 2, 1], [1, 4, 7, 1, 2, 2, 1, 1, 1, 4, 2, 1],
        [1, 4, 5, 1, 2, 2, 1, 1, 1, 4, 6, 1], [1, 4, 7, 1, 2, 2, 1, 1, 1, 4, 6, 1]
        ],
        # E =-16.1
        [
        [2, 5, 4, 1, 1, 3, 1, 1, 1, 5, 7, 1], [2, 5, 2, 1, 1, 3, 1, 1, 1, 5, 3, 1],
        [2, 5, 4, 1, 1, 3, 1, 1, 1, 5, 3, 1], [2, 5, 2, 1, 1, 3, 1, 1, 1, 5, 7, 1]
        ],
        # E = -15.9
        [
        [1, 4, 1, 1, 2, 2, 1, 1, 1, 4, 2, 1], [1, 4, 3, 1, 2, 2, 1, 1, 1, 4, 2, 1],
        [1, 4, 6, 1, 2, 2, 1, 1, 1, 4, 2, 1], [1, 4, 3, 1, 2, 2, 1, 1, 1, 4, 6, 1],
        [1, 4, 1, 1, 2, 2, 1, 1, 1, 4, 6, 1], [1, 4, 6, 1, 2, 2, 1, 1, 1, 4, 6, 1]
        ],
        # E = -15.6
        [
        [2, 5, 3, 1, 1, 3, 1, 1, 1, 5, 3, 1], [2, 5, 3, 1, 1, 3, 1, 1, 1, 5, 7, 1],
        [2, 5, 8, 1, 1, 3, 1, 1, 1, 5, 3, 1], [2, 5, 6, 1, 1, 3, 1, 1, 1, 5, 7, 1],
        [2, 5, 6, 1, 1, 3, 1, 1, 1, 5, 3, 1], [2, 5, 8, 1, 1, 3, 1, 1, 1, 5, 7, 1]
        ],
        # E = -15.4
        [
        [1, 4, 7, 1, 2, 2, 1, 1, 1, 2, 6, 1], [1, 4, 5, 1, 2, 2, 1, 1, 1, 2, 6, 1]
        ],
    ]

    deg = Dict(
        1 => 1, 2 => 1, 3 => 1, 4 => 1,
        5 => 2, 6 => 2, 7 => 2, 8 => 2,
        9 => 3, 10 => 3, 11 => 3, 12 => 3, 13 => 3, 14 => 3,
        15 => 4, 16 => 4, 17 => 4, 18 => 4, 19 => 4, 20 => 4,
        21 => 5, 22 => 5,
    )

    instance = "$(@__DIR__)/instances/pathological/chim_$(m)_$(n)_$(t).txt"

    ig = ising_graph(instance)

    fg = factor_graph(
        ig,
        spectrum=full_spectrum,
        cluster_assignment_rule=super_square_lattice((m, n, t))
    )

    for transform ∈ all_lattice_transformations
        peps = PEPSNetwork{false}(m, n, fg, transform, β=β)
        #update_gauges!(peps, :id)
        sol = low_energy_spectrum(peps, num_states)

        @test sol.energies ≈ exact_energies
        for (i, σ) ∈ enumerate(sol.states)
            @test σ ∈ exact_states[deg[i]]
        end
    end
end


text:
"""
    pedsort(src, unknown, out)
---

# Description
Given a input stream `src` contains columns:
1. id
2. sire
3. dam
4. other information

And a `Set` of `unknown` strings, this program will sort the pedigree, such that `Sire` and `Dam` appear in
the `ID` column first.  In other words, an offspring appears after its parents in the sorted pedigree.

Results will be written to `out`, which has columns:
1. recoded `sire` name
2. recoded `dam` name
3. original name of current `ID`, which is the current row number (starts from 1)
4. `other information` in the orginal stream.

All former `unknown`s are coded as `0`.
"""
function pedsort(src, unknown::Set{String}, out)
    title("Sort pedigree in $src")
    #for line
end


text:
push!(LOAD_PATH, "/mnt/g/Python Programme/BALL.jl/src" )
using BALL
using BenchmarkTools
using StaticArrays
pdb_filepath = "1ake.pdb"

comp = System(pdb_filepath)
dfs = DataFrameSystem(comp)

dfs2 = DataFrameSystem(comp)
for x in dfs.atoms.position_
    x = SVector{3,Float64}(x[1]/sqrt(3),x[2]/sqrt(3),x[3]/sqrt(3))
end

rmsd(dfs, dfs2, calphaonly=false)

elapsed = @elapsed rmsd(dfs,dfs2)

println(elapsed)
@show @benchmark rmsd(dfs,dfs2)


text:
include("CreateMesh.jl")

#=== Case 1 layered medium ===#
vp = [5000.0,5300,5800,6000]
vs = [3500.0,3200,4000,4500]
rho = 2.8
depth = [700,200,400,600]
Horizon = 1000
pkf = 10
T = 2
ext = 10
iflag = 2
dx = nothing
dz = nothing
dt = nothing
nT = nothing

Model = model(vp,vs,rho,depth,Horizon,pkf,T,ext,iflag,dx,dz,dt,nT)

nwf = initnelwf(Model.nDZ,Model.nHX,Model.ext,Model.iflag)

sn = 2
loc = [20 20;30 30]
ot = [0.1, 0.2]
mt = zeros(2,3)
waveform = Ricker(30,Model.dt)
waveform = [waveform waveform]

sou = initsource(sn,loc,ot,mt,waveform,Model,nwf)


text:
using Test, LinearAlgebra, MPStates

@testset "Operations with Mps" begin
rtest1 = MPStates.testMps("rtest1")
ctest1 = MPStates.testMps("ctest1")
rtest2 = MPStates.testMps("rtest2")
ctest2 = MPStates.testMps("ctest2")
rtest3 = MPStates.testMps("rtest3")
@testset "measure occupation at one site" begin
    @test expected(rtest1, "n", 1) ≈ 1/9
    @test expected(rtest1, "n", 2) ≈ 4/9
    @test expected(rtest1, "n", 3) ≈ 4/9
    @test expected(rtest1, "n", 4) ≈ 1.
    @test expected(rtest1, "n", 5) ≈ 0.64
    @test expected(rtest1, "n", 6) ≈ 1.
    @test expected(ctest1, "n", 1) ≈ 1/9
    @test expected(ctest1, "n", 2) ≈ 4/9
    @test expected(ctest1, "n", 3) ≈ 4/9
    @test expected(ctest1, "n", 4) ≈ 1.
    @test expected(ctest1, "n", 5) ≈ 0.64
    @test expected(ctest1, "n", 6) ≈ 1.
    @test expected(rtest2, "n", 1) ≈ 5/9
    @test expected(rtest2, "n", 2) ≈ 5/9
    @test expected(rtest2, "n", 3) ≈ 1.
    @test expected(rtest2, "n", 4) ≈ 1.
    @test expected(rtest2, "n", 5) ≈ 1.
    @test expected(rtest2, "n", 6) ≈ 0.64
    @test expected(ctest2, "n", 1) ≈ 5/9
    @test expected(ctest2, "n", 2) ≈ 5/9
    @test expected(ctest2, "n", 3) ≈ 1.
    @test expected(ctest2, "n", 4) ≈ 1.
    @test expected(ctest2, "n", 5) ≈ 1.
    @test expected(ctest2, "n", 6) ≈ 0.64
    @test expected(rtest3, "Sz", 1) ≈ -1/6
    @test expected(rtest3, "Sz", 2) ≈ 1/2
    @test expected(rtest3, "Sz", 3) ≈ 1/6
    @test expected(rtest3, "Sz", 4) ≈ -1/6
    @test expected(rtest3, "Sz", 5) ≈ 1/2
    @test expected(rtest3, "Sz", 6) ≈ 1/6
    @test expected(rtest3, "S+", 1) ≈ -sqrt(2)/6
    @test expected(rtest3, "S+", 2) ≈ 0. atol = 1e-15
    @test expected(rtest3, "S+", 3) ≈ 0. atol = 1e-15
    @test expected(rtest3, "S+", 4) ≈ -sqrt(2)/6
    @test expected(rtest3, "S+", 5) ≈ 0. atol = 1e-15
    @test expected(rtest3, "S+", 6) ≈ 0. atol = 1e-15
    @test expected(rtest3, "S-", 1) ≈ -sqrt(2)/6
    @test expected(rtest3, "S-", 2) ≈ 0. atol = 1e-15
    @test expected(rtest3, "S-", 3) ≈ 0. atol = 1e-15
    @test expected(rtest3, "S-", 4) ≈ -sqrt(2)/6
    @test expected(rtest3, "S-", 5) ≈ 0. atol = 1e-15
    @test expected(rtest3, "S-", 6) ≈ 0. atol = 1e-15
end

@testset "measure fermionic correlations" begin
    @test expected(rtest1, "c+", 1, "c", 2, ferm_op="Z") ≈ -2/9
    @test expected(rtest1, "c+", 1, "c", 4, ferm_op="Z") ≈ 0. atol=1e-15
    @test expected(rtest1, "c+", 3, "c", 2, ferm_op="Z") ≈ -4/9
    @test expected(rtest1, "c+", 2, "c", 6, ferm_op="Z") ≈ 0. atol=1e-15
    @test expected(rtest1, "c+", 6, "c", 1, ferm_op="Z") ≈ 0.  atol=1e-15
    # On-site operations.
    @test expected(rtest1, "c+", 1, "c+", 1) ≈ 0.  atol=1e-15
    @test expected(rtest1, "c+", 1, "c", 1) ≈ 1/9
    @test expected(rtest1, "c+", 3, "c", 3) ≈ 4/9
    @test expected(ctest1, "c+", 1, "c", 2, ferm_op="Z") ≈ complex(0., 2/9)
    @test expected(ctest1, "c+", 1, "c", 4, ferm_op="Z") ≈ 0. atol=1e-15
    @test expected(ctest1, "c+", 3, "c", 2, ferm_op="Z") ≈ complex(0., 4/9)
    @test expected(ctest1, "c+", 2, "c", 6, ferm_op="Z") ≈ 0. atol=1e-15
    @test expected(ctest1, "c+", 6, "c", 1, ferm_op="Z") ≈ 0.  atol=1e-15
    @test expected(ctest1, "c", 2, "c", 6) ≈ 0. atol=1e-15
    @test expected(ctest1, "c+", 6, "c", 6) ≈ 1.
    @test expected(rtest2, "c+", 1, "c", 2, ferm_op="Z") ≈ -4/9
    @test expected(rtest2, "c+", 1, "c", 4, ferm_op="Z") ≈ 0. atol=1e-15
    @test expected(rtest2, "c+", 2, "c", 6, ferm_op="Z") ≈ 2/9*0.6*0.8
    @test expected(rtest2, "c+", 6, "c", 1, ferm_op="Z") ≈ 2/9*0.6*0.8
    # On-site operations.
    @test expected(rtest2, "c+", 2, "c", 2) ≈ 5/9
    @test expected(ctest2, "c+", 1, "c", 2, ferm_op="Z") ≈ complex(0., -4/9)
    @test expected(ctest2, "c+", 1, "c", 4, ferm_op="Z") ≈ 0. atol=1e-15
    @test expected(ctest2, "c+", 2, "c", 6, ferm_op="Z") ≈ 2/9*0.6*0.8
    @test expected(ctest2, "c+", 6, "c", 1, ferm_op="Z") ≈ complex(0., -2/9*0.6*0.8)
end

@testset "measure correlations" begin
    @test expected(rtest1, "b+", 1, "b", 2) ≈ -2/9
    @test expected(rtest1, "b+", 1, "b", 4) ≈ 0. atol=1e-15
    @test expected(rtest1, "b+", 3, "b", 2) ≈ -4/9
    @test expected(rtest1, "b+", 2, "b", 6) ≈ 0. atol=1e-15
    @test expected(rtest1, "b+", 6, "b", 1) ≈ 0.  atol=1e-15
    # On-site operations.
    @test expected(rtest1, "b+", 1, "b+", 1) ≈ 0.  atol=1e-15
    @test expected(rtest1, "b+", 1, "b", 1) ≈ 1/9
    @test expected(rtest1, "b+", 3, "b", 3) ≈ 4/9
    @test expected(ctest1, "b+", 1, "b", 2) ≈ complex(0., 2/9)
    @test expected(ctest1, "b+", 1, "b", 4) ≈ 0. atol=1e-15
    @test expected(ctest1, "b+", 3, "b", 2) ≈ complex(0., 4/9)
    @test expected(ctest1, "b+", 2, "b", 6) ≈ 0. atol=1e-15
    @test expected(ctest1, "b+", 6, "b", 1) ≈ 0.  atol=1e-15
    # On-site operations.
    @test expected(ctest1, "b", 2, "b", 6) ≈ 0. atol=1e-15
    @test expected(ctest1, "b+", 6, "b", 6) ≈ 1.
    @test expected(rtest2, "b+", 1, "b", 2) ≈ -4/9
    @test expected(rtest2, "b+", 1, "b", 4) ≈ 0. atol=1e-15
    @test expected(rtest2, "b+", 2, "b", 6) ≈ -2/9*0.6*0.8
    @test expected(rtest2, "b+", 6, "b", 1) ≈ 2/9*0.6*0.8
    # On-site operations.
    @test expected(rtest2, "b+", 2, "b", 2) ≈ 5/9
    @test expected(ctest2, "b+", 1, "b", 2) ≈ complex(0., -4/9)
    @test expected(ctest2, "b+", 1, "b", 4) ≈ 0. atol=1e-15
    @test expected(ctest2, "b+", 2, "b", 6) ≈ -2/9*0.6*0.8
    @test expected(ctest2, "b+", 6, "b", 1) ≈ complex(0., -2/9*0.6*0.8)
end

@testset "measure spin correlations" begin
    @test expected(rtest3, "Sz", 1, "Sz", 2) ≈ 1/6
    @test expected(rtest3, "Sz", 1, "Sz", 3) ≈ 1/6
    @test expected(rtest3, "Sz", 1, "Sz", 4) ≈ 1/36
    @test expected(rtest3, "Sz", 1, "Sz", 1) ≈ 1/2
    @test expected(rtest3, "Sz", 2, "Sz", 5) ≈ 1/4
    @test expected(rtest3, "Sz", 2, "Sz", 2) ≈ 5/6
    @test expected(rtest3, "S+", 2, "Sz", 1) ≈ 0. atol=1e-15
    @test expected(rtest3, "S+", 2, "Sz", 3) ≈ 0. atol=1e-15
    @test expected(rtest3, "S+", 1, "Sz", 1) ≈ 0. atol=1e-15
    @test expected(rtest3, "S+", 1, "Sz", 2) ≈ -sqrt(2)/6
    @test expected(rtest3, "S+", 1, "Sz", 3) ≈ -sqrt(2)/6
    @test expected(rtest3, "S+", 4, "Sz", 6) ≈ -sqrt(2)/6
    @test expected(rtest3, "S+", 1, "Sz", 5) ≈ -sqrt(2)/12
    @test expected(rtest3, "S+", 1, "S-", 2) ≈ -1/3
    @test expected(rtest3, "S+", 2, "S-", 1) ≈ -1/3
    @test expected(rtest3, "S+", 1, "S-", 3) ≈ 1/3
    @test expected(rtest3, "S+", 3, "S-", 1) ≈ 1/3
    @test expected(rtest3, "S+", 2, "S-", 3) ≈ -1/3
    @test expected(rtest3, "S+", 3, "S-", 2) ≈ -1/3
    @test expected(rtest3, "S+", 1, "S-", 1) ≈ 4/3
    @test expected(rtest3, "S+", 2, "S-", 2) ≈ 5/3
    @test expected(rtest3, "S+", 3, "S-", 3) ≈ 5/3
    @test expected(rtest3, "S-", 4, "S+", 4) ≈ 5/3
    @test expected(rtest3, "S-", 5, "S+", 5) ≈ 2/3
    @test expected(rtest3, "S-", 6, "S+", 6) ≈ 4/3
end

@testset "measure 2 point occupations" begin
    @test expected(rtest1, "n", 1, "n", 2) ≈ 0. atol=1e-15
    @test expected(rtest1, "n", 1, "n", 4) ≈ 1/9
    @test expected(rtest1, "n", 3, "n", 2) ≈ 0. atol=1e-15
    @test expected(rtest1, "n", 2, "n", 5) ≈ 4/9*0.64
    @test expected(rtest1, "n", 6, "n", 1) ≈ 1/9
    @test expected(ctest1, "n", 1, "n", 2) ≈ 0. atol=1e-15
    @test expected(ctest1, "n", 1, "n", 4) ≈ 1/9
    @test expected(ctest1, "n", 3, "n", 2) ≈ 0. atol=1e-15
    @test expected(ctest1, "n", 2, "n", 5) ≈ 4/9*0.64
    @test expected(ctest1, "n", 6, "n", 1) ≈ 1/9
    @test expected(rtest2, "n", 1, "n", 2) ≈ 1/9
    @test expected(rtest2, "n", 1, "n", 4) ≈ 5/9
    @test expected(rtest2, "n", 2, "n", 6) ≈ 5/9*0.64
    @test expected(rtest2, "n", 6, "n", 1) ≈ 5/9*0.64
    @test expected(ctest2, "n", 1, "n", 2) ≈ 1/9
    @test expected(ctest2, "n", 1, "n", 4) ≈ 5/9
    @test expected(ctest2, "n", 2, "n", 6) ≈ 5/9*0.64
    @test expected(ctest2, "n", 6, "n", 1) ≈ 5/9*0.64
end

@testset "contraction of two MPS" begin
    L = 6
    GHZ = Mps(L, "GHZ")
    W = Mps(L, "W")
    full = Mps(L, "full")
    product = Mps(L, "product")
    cGHZ = Mps(ComplexF64, L, "GHZ")
    cW = Mps(ComplexF64, L, "W")
    cfull = Mps(ComplexF64, L, "full")
    cproduct = Mps(ComplexF64, L, "product")
    @test contract(GHZ, W) ≈ 0. atol=1e-15
    @test contract(GHZ, full) ≈ 1/sqrt(2^(L-1))
    @test contract(GHZ, product) ≈ 1/sqrt(2)
    @test contract(W, full) ≈ sqrt(L/2^L)
    @test contract(W, product) ≈ 0. atol=1e-15
    @test contract(full, product) ≈ 1/sqrt(2^L)
    @test contract(rtest1, W) ≈ 0. atol=1e-15
    @test contract(rtest1, GHZ) ≈ 0. atol=1e-15
    @test contract(rtest1, full) ≈ (1/3 - 2/3 + 2/3)*(-0.6 + 0.8)/sqrt(2^L)
    @test contract(ctest1, cW) ≈ 0. atol=1e-15
    @test contract(ctest1, cGHZ) ≈ 0. atol=1e-15
    @test contract(ctest1, cfull) ≈ (-1im/3 - 2im/3 - 2/3)*(0.6im + 0.8)/sqrt(2^L)
    @test contract(cfull, ctest1) ≈ (1im/3 + 2im/3 - 2/3)*(-0.6im + 0.8)/sqrt(2^L)
    @test contract(rtest2, W) ≈ 0. atol=1e-15
    @test contract(rtest2, GHZ) ≈ 0.8*1/3/sqrt(2)
    @test contract(rtest2, full) ≈ (1/3 - 2/3 + 2/3)*(0.6 + 0.8)/sqrt(2^L)
    @test contract(ctest2, cW) ≈ 0. atol=1e-15
    @test contract(ctest2, cGHZ) ≈ -0.8*1/3/sqrt(2)
    @test contract(cGHZ, ctest2) ≈ -0.8*1/3/sqrt(2)
    @test contract(ctest2, cfull) ≈ (-1im/3 - 2/3 - 2im/3)*(0.6 - 0.8im)/sqrt(2^L)
    @test contract(cfull, ctest2) ≈ (1im/3 - 2/3 + 2im/3)*(0.6 + 0.8im)/sqrt(2^L)
end

@testset "norm of a MPS" begin
    L = 6
    GHZ = Mps(L, "GHZ")
    W = Mps(L, "W")
    full = Mps(L, "full")
    product = Mps(L, "product")
    AKLT = Mps(L, "AKLT")
    @test norm(GHZ) ≈ 1.
    @test norm(full) ≈ 1.
    @test norm(product) ≈ 1.
    @test norm(W) ≈ 1.
    @test norm(AKLT) ≈ 1.
    @test norm(rtest1) ≈ 1.
    @test norm(rtest2) ≈ 1.
    @test norm(ctest1) ≈ 1.
    @test norm(ctest2) ≈ 1.
    @test norm(rtest3) ≈ 1.
end

@testset "schmidt decomposition" begin
    L = 4
    GHZ = Mps(L, "GHZ")
    @test MPStates.schmidt_decomp(GHZ, 1) ≈ [1/sqrt(2), 1/sqrt(2)]
    @test MPStates.schmidt_decomp(GHZ, 2) ≈ [1/sqrt(2), 1/sqrt(2)]
    W = Mps(L, "W")
    @test MPStates.schmidt_decomp(W, 1) ≈ [sqrt(3)/2, 0.5]
    @test MPStates.schmidt_decomp(W, 2) ≈ [1/sqrt(2), 1/sqrt(2)]
    @test MPStates.schmidt_decomp(rtest1, 1) ≈ [2sqrt(2)/3, 1/3]
    @test MPStates.schmidt_decomp(rtest1, 2) ≈ [sqrt(5)/3, 2/3]
    @test MPStates.schmidt_decomp(rtest1, 3) ≈ [1., 0.]
    @test MPStates.schmidt_decomp(rtest1, 4) ≈ [1., 0.]
    @test MPStates.schmidt_decomp(rtest1, 5) ≈ [1., 0.]
    @test MPStates.schmidt_decomp(rtest1, 6) ≈ [1.]
    @test MPStates.schmidt_decomp(ctest1, 1) ≈ [2sqrt(2)/3, 1/3]
    @test MPStates.schmidt_decomp(ctest1, 2) ≈ [sqrt(5)/3, 2/3]
    @test MPStates.schmidt_decomp(ctest1, 3) ≈ [1., 0.]
    @test MPStates.schmidt_decomp(ctest1, 4) ≈ [1., 0.]
    @test MPStates.schmidt_decomp(ctest1, 5) ≈ [1., 0.]
    @test MPStates.schmidt_decomp(ctest1, 6) ≈ [1.]
    @test MPStates.schmidt_decomp(rtest2, 1) ≈ [sqrt((1+sqrt(17)/9)/2),
                                                sqrt((1-sqrt(17)/9)/2)]
    @test MPStates.schmidt_decomp(rtest2, 2) ≈ [1., 0.]
    @test MPStates.schmidt_decomp(rtest2, 3) ≈ [1., 0.]
    @test MPStates.schmidt_decomp(rtest2, 4) ≈ [1., 0.]
    @test MPStates.schmidt_decomp(rtest2, 5) ≈ [1., 0.]
    @test MPStates.schmidt_decomp(rtest2, 6) ≈ [1.]
    @test MPStates.schmidt_decomp(ctest2, 1) ≈ [sqrt((1+sqrt(17)/9)/2),
                                                sqrt((1-sqrt(17)/9)/2)]
    @test MPStates.schmidt_decomp(ctest2, 2) ≈ [1., 0.]
    @test MPStates.schmidt_decomp(ctest2, 3) ≈ [1., 0.]
    @test MPStates.schmidt_decomp(ctest2, 4) ≈ [1., 0.]
    @test MPStates.schmidt_decomp(ctest2, 5) ≈ [1., 0.]
    @test MPStates.schmidt_decomp(ctest2, 6) ≈ [1.]
    A1 = 1/sqrt(6)*[[0. -1. 1. 0. 0.];
                    [1. 0. 0. 1. -1.];
                    [0. 0. 0. -1. 0.]]
    A2 = 1/sqrt(6)*[[0. -1. 0. 0. 0.];
                    [1. 0. 0. -1. 1.];
                    [0. 1. -1. 0. 0.]]
    @test MPStates.schmidt_decomp(rtest3, 1) ≈ svdvals(A1)
    @test MPStates.schmidt_decomp(rtest3, 2) ≈ svdvals(A2)
    @test MPStates.schmidt_decomp(rtest3, 3) ≈ [1.]
    @test MPStates.schmidt_decomp(rtest3, 4) ≈ svdvals(A1)
    @test MPStates.schmidt_decomp(rtest3, 5) ≈ svdvals(A2)
    @test MPStates.schmidt_decomp(rtest3, 6) ≈ [1.]
end

@testset "entanglement entropy" begin
    L = 4
    GHZ = Mps(L, "GHZ")
    @test ent_entropy(GHZ, 1) ≈ 1/sqrt(2)
    @test ent_entropy(GHZ, 2) ≈ 1/sqrt(2)
    W = Mps(L, "W")
    @test ent_entropy(W, 1) ≈ 0.5 - sqrt(3)/2*log2(sqrt(3)/2)
    @test ent_entropy(W, 2) ≈ 1/sqrt(2)
    @test ent_entropy(rtest1, 1) ≈ -2sqrt(2)/3*log2(2sqrt(2)/3) - 1/3*log2(1/3)
    @test ent_entropy(rtest1, 2) ≈ -sqrt(5)/3*log2(sqrt(5)/3) - 2/3*log2(2/3)
    @test ent_entropy(rtest1, 3) ≈ 0. atol=1e-15
    @test ent_entropy(rtest1, 4) ≈ 0. atol=1e-15
    @test ent_entropy(ctest1, 1) ≈ -2sqrt(2)/3*log2(2sqrt(2)/3) - 1/3*log2(1/3)
    @test ent_entropy(ctest1, 2) ≈ -sqrt(5)/3*log2(sqrt(5)/3) - 2/3*log2(2/3)
    @test ent_entropy(ctest1, 3) ≈ 0. atol=1e-15
    @test ent_entropy(ctest1, 4) ≈ 0. atol=1e-15
    @test ent_entropy(rtest2, 1) ≈ (
        -sqrt((1+sqrt(17)/9)/2)*log2(sqrt((1+sqrt(17)/9)/2))
        -sqrt((1-sqrt(17)/9)/2)*log2(sqrt((1-sqrt(17)/9)/2))
        )
    @test ent_entropy(rtest2, 2) ≈ 0. atol=1e-15
    @test ent_entropy(rtest2, 3) ≈ 0. atol=1e-15
    @test ent_entropy(rtest2, 4) ≈ 0. atol=1e-15
    @test ent_entropy(ctest2, 1) ≈ (
        -sqrt((1+sqrt(17)/9)/2)*log2(sqrt((1+sqrt(17)/9)/2))
        -sqrt((1-sqrt(17)/9)/2)*log2(sqrt((1-sqrt(17)/9)/2))
        )
    @test ent_entropy(ctest2, 2) ≈ 0. atol=1e-15
    @test ent_entropy(ctest2, 3) ≈ 0. atol=1e-15
    @test ent_entropy(ctest2, 4) ≈ 0. atol=1e-15
    A1 = 1/sqrt(6)*[[0. -1. 1. 0. 0.];
                    [1. 0. 0. 1. -1.];
                    [0. 0. 0. -1. 0.]]
    svals1 = svdvals(A1)
    A2 = 1/sqrt(6)*[[0. -1. 0. 0. 0.];
                    [1. 0. 0. -1. 1.];
                    [0. 1. -1. 0. 0.]]
    svals2 = svdvals(A2)
    @test ent_entropy(rtest3, 1) ≈ sum(-svals1.*log2.(svals1))
    @test ent_entropy(rtest3, 2) ≈ sum(-svals2.*log2.(svals2))
    @test ent_entropy(rtest3, 3) ≈ 0. atol=1e-15
    @test ent_entropy(rtest3, 4) ≈ sum(-svals1.*log2.(svals1))
    @test ent_entropy(rtest3, 5) ≈ sum(-svals2.*log2.(svals2))
    @test ent_entropy(rtest3, 6) ≈ 0. atol=1e-15
end

@testset "enlargement of bond dimension of MPS" begin
    L = 10
    GHZ = Mps(L, "GHZ")
    enlarge_bond_dimension!(GHZ, 5)
    @test size(GHZ.M[1], 1) == 1
    @test size(GHZ.M[1], 3) == 2
    @test size(GHZ.M[3], 1) == 4
    @test size(GHZ.M[3], 3) == 5
    for i=4:7
        @test size(GHZ.M[i], 1) == 5
        @test size(GHZ.M[i], 3) == 5
    end
    @test size(GHZ.M[8], 1) == 5
    @test size(GHZ.M[8], 3) == 4
    @test size(GHZ.M[L], 1) == 2
    @test size(GHZ.M[L], 3) == 1
    # Check that the properties of the Mps are left intact.
    full = Mps(L, "full")
    @test contract(GHZ, full) ≈ 1/sqrt(2^(L-1))
    enlarge_bond_dimension!(full, 11)
    @test contract(GHZ, full) ≈ 1/sqrt(2^(L-1))
end

@testset "SVD truncation of MPS" begin
    svd_truncate!(rtest1, 1)
    @test expected(rtest1, "n", 1) ≈ 0.  atol=1e-15
    @test (expected(rtest1, "n", 2)^2 + expected(rtest1, "n", 3)^2) ≈ 1.
    @test expected(rtest1, "n", 4) ≈ 1.
    @test expected(rtest1, "n", 5) ≈ 0.64
    @test expected(rtest1, "n", 6) ≈ 1.
    svd_truncate!(ctest1, 1)
    @test expected(ctest1, "n", 1) ≈ 0.  atol=1e-15
    @test (expected(ctest1, "n", 2)^2 + expected(ctest1, "n", 3)^2) ≈ 1.
    @test expected(ctest1, "n", 4) ≈ 1.
    @test expected(ctest1, "n", 5) ≈ 0.64
    @test expected(ctest1, "n", 6) ≈ 1.
    svd_truncate!(rtest2, 1)
    @test expected(rtest2, "n", 1) ≈ (0.7882054380161092)^2
    @test expected(rtest2, "n", 2) ≈ (0.7882054380161092)^2
    @test expected(rtest2, "n", 3) ≈ 1.
    @test expected(rtest2, "n", 4) ≈ 1.
    @test expected(rtest2, "n", 5) ≈ 1.
    @test expected(rtest2, "n", 6) ≈ 0.64
    svd_truncate!(ctest2, 1)
    @test expected(ctest2, "n", 1) ≈ (0.7882054380161092)^2
    @test expected(ctest2, "n", 2) ≈ (0.7882054380161092)^2
    @test expected(ctest2, "n", 3) ≈ 1.
    @test expected(ctest2, "n", 4) ≈ 1.
    @test expected(ctest2, "n", 5) ≈ 1.
    @test expected(ctest2, "n", 6) ≈ 0.64
    crtest3 = deepcopy(rtest3)
    A2 = 1/sqrt(6)*[[0. -1. 0. 0. 0.];
                    [1. 0. 0. -1. 1.];
                    [0. 1. -1. 0. 0.]]
    u, s, vt = svd(A2)
    svd_truncate!(crtest3, 1)
    @test expected(crtest3, "Sz", 3) ≈ -u[1, 1] + u[3, 1]
end

@testset "save and read Mps in hdf5 format" for T in [Float64, ComplexF64]
    L = 10
    GHZ = Mps(T, L, "GHZ")

    filename = "foo.h5"
    # Remove file if previous test crashed and file was not removed.
    isfile(filename) && rm(filename)

    save_mps(filename, GHZ)
    # Read Mps.
    psi = read_mps(filename)
    @test psi.L == L
    @test psi.d == 2
    @test eltype(psi.M[1]) == T
    for i=1:L
        @test psi.M[i] ≈ GHZ.M[i]
    end
    # Remove hdf5 testing file.
    isfile(filename) && rm(filename)
end
end # @testset "Operations with Mps"


text:
include("./solve_ipopt.jl")
using QPnorm
using Random
using DataFrames, CSV

rng = MersenneTwister(123)

function compute_metrics(P::Matrix{T}, q::Vector{T},
    A::Matrix{T}, b::Vector{T}, r::T,
    x::Vector{T}, λ::Vector{T}) where {T}

    m, n = size(A)

    f = dot(x, P*x)/2 + dot(x, q)

    grad_residual = norm(P*x + q + A'*λ[1:end-1] + λ[end]*x, Inf)
    infeasibility = max(maximum(A*x - b), abs(norm(x)^2 - r^2), 0)
    complementarity = maximum(minimum([abs.(λ[1:end-1]) abs.(A*x - b)], dims=2))
    dual_infeasibility = max(-minimum(λ[1:end-1]), 0.0)

    active_set = λ .>= 1e-8
    active_set[end] = true # Norm constraint is an equality constraint; thus it should always be included
    V = nullspace([A; x'][active_set, :])
    if length(V) > 0
        min_eig = minimum(eigvals(Symmetric(V'*(P + I*λ[end])*V)))
    else
        min_eig = 0.0
    end

    return f, grad_residual, infeasibility, dual_infeasibility, complementarity, min_eig
end

df = DataFrame(n = Int[], m = Int[],
    time = Float64[], time_ipopt = Float64[],
    objective = Float64[], objective_ipopt = Float64[],
    infeasibility = Float64[], infeasibility_ipopt = Float64[],
    dual_infeasibility = Float64[], dual_infeasibility_ipopt = Float64[],
    complementarity = Float64[], complementarity_ipopt = Float64[],
    grad_residual = Float64[], grad_residual_ipopt = Float64[],
    min_eig = Float64[], min_eig_ipopt = Float64[]
)

dimensions = Vector(range(100, 2000, length=10))
prepend!(dimensions, 10) # Small problem to warm up solvers
r = 100.0

for dimension in dimensions
    n = Int(floor(dimension))
    m = Int(floor(1.5*n))
    P = randn(rng, n, n); P = (P + P')/2
    q = randn(rng, n)
    A = randn(rng, m, n); b = randn(rng, m)

    x_ipopt = Float64[]; f_ipopt = NaN; infeasibility_ipopt = NaN; t_ipopt = NaN
    grad_residual_ipopt = NaN; min_eig_ipopt = NaN; multipliers_ipopt = Float64[];
    complementarity_ipopt = NaN; dual_infeasibility_ipopt = NaN
    try
        x_ipopt, multipliers_ipopt, t_ipopt = solve_ipopt(P, q, A, b, r, r; print_level=0)
        f_ipopt, grad_residual_ipopt, infeasibility_ipopt, dual_infeasibility_ipopt, complementarity_ipopt, min_eig_ipopt = compute_metrics(Matrix(P), q, Matrix(A), b, r, x_ipopt, multipliers_ipopt)
    catch e
        nothing
    end

    x = Float64[]; f = NaN; infeasibility = NaN; t = NaN
    grad_residual = NaN; min_eig = NaN; multipliers = Float64[];
    complementarity = NaN; dual_infeasibility = NaN
    try
        t = @elapsed x_init = find_feasible_point(A, b, r, r)
        t += @elapsed x, multipliers = QPnorm.solve_boundary(P, q, A, b, r, x_init, verbosity=1, printing_interval=5000, max_iter=5000)
        t = @elapsed x_init = find_feasible_point(A, b, r, r)
        t += @elapsed x, multipliers = QPnorm.solve_boundary(P, q, A, b, r, x_init, verbosity=1, printing_interval=5000, max_iter=5000)
        f, grad_residual, infeasibility, dual_infeasibility, complementarity, min_eig = compute_metrics(P, q, A, b, r, x, multipliers)
    catch e
        nothing
    end

    push!(df, [n, m,
            t, t_ipopt,
            f, f_ipopt,
            infeasibility, infeasibility_ipopt,
            dual_infeasibility, dual_infeasibility_ipopt,
            complementarity, complementarity_ipopt,
            grad_residual, grad_residual_ipopt,
            min_eig, min_eig_ipopt])

    df |> CSV.write(string("results_random.csv"))
end


text:
using GVFHordes
import BSON

struct LearnedPolicy{M, D} <: GVFHordes.GVFParamFuncs.AbstractPolicy
    policy::M
    info::D
end

function LearnedPolicy(loc::AbstractString)
    d = BSON.load(loc, @__MODULE__)
    if "LearnedPolicy" ∈ keys(d)
        d["LearnedPolicy"]
    elseif "learner" ∈ keys(d)
        LearnedPolicy(Flux.Chain(d["fc"], d["learner"]), d["info"])
    end
end

Base.get(π::LearnedPolicy, state_t, action_t) = π(state_t, action_t)
# Base.get(π::LearnedPolicy; kwargs...) = function(;kwargs...)
#     return π(kwargs[:state_t], kwargs[:action_t])
# end

# Base.get(π::LearnedPolicy; kwargs...) = π(kwargs[:state_t], kwargs[:action_t])

function Base.get(π::LearnedPolicy; kwargs...)
    # @show π(kwargs[:state_t])
    # @show kwargs[:action_t]
    # @show kwargs[:state_t]
    # @show (kwargs[:action_t], findmax(π(kwargs[:state_t]))[2])
    return Int(kwargs[:action_t] == findmax(π(kwargs[:state_t]))[2])
end

(lp::LearnedPolicy)(s) = lp.policy(s)
(lp::LearnedPolicy)(s, a) = lp.policy(s)[a]

# Policy learning struct
mutable struct PolicyLearner{L<:Learner, FC, E, D, C, Φ, DEMON} <: MinimalRLCore.AbstractAgent
    learner::L
    fc::FC
    exploration::E

    discount::D
    cumulant::C
    o_t::Φ
    a_t::Int
    BehaviourDemon::DEMON
end

function MinimalRLCore.start!(pl::PolicyLearner, o)
    s_t = pl.fc(o)
    pl.o_t = copy(o)
    pl.a_t = μ(pl, s_t)

    pl.a_t
end

function μ_dist(pl::PolicyLearner, state)
    qs = pl.learner(state)
    pl.exploration(qs)
end

function μ(pl::PolicyLearner, state)
    action_probs = μ_dist(pl, state)
    action = sample(1:length(action_probs), Weights(action_probs))
end

function MinimalRLCore.step!(pl::PolicyLearner, o_tp1, rew, term)
    s_tp1 = pl.fc(o_tp1)
    next_action = μ(pl, s_tp1)

    # rew = pl.cumulant(;r=rew, o_t=pl.o_t, a_t=pl.a_t, o_tp1=o_tp1)
    # get(pl.cumulant;r=rew, state_t=pl.o_t, action_t=pl.a_t, state_tp1=o_tp1)
    # γ = if pl.discount isa Number
    #     pl.discount
    # else
    #     pl.discount(;r=rew, o_t=pl.o_t, a_t=pl.a_t, o_tp1=o_tp1)
    # end

    update!(pl.learner,
            pl.BehaviourDemon,
            pl.o_t,
            o_tp1,
            pl.fc(pl.o_t),
            pl.a_t,
            pl.fc(o_tp1),
            next_action,
            term,
            (state, obs) -> μ_dist(pl, state),
            rew)

    pl.o_t = copy(o_tp1)
    pl.a_t = next_action

    pl.a_t
end

save(pl::PolicyLearner, loc, info::AbstractDict) =
    BSON.bson(loc,
              Dict("learner"=>pl.learner,
                   "fc"=>pl.fc,
                   "exp"=>pl.exploration,
                   "info"=>info))

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
# *******************************************
#  Checker
# *******************************************
function check_config_consistency(c1::ConfigFFT, c2::ConfigFFT)
    if !(c1 === c2)
        errmsg = "ConfigFFTInconsistencyError"
        throw(ErrorException(errmsg))
    end
end

function check_size_consistency(
        vals::Array, config::ConfigFFT)

    if size(vals) != config.ngrids
        errmsg = (
            "size(vals) and config.ngrids must match:\n"
            * "size(vals)=$(size(vals)), and "
            * "config.ngrids=$(config.ngrids)"
        )
        throw(ErrorException(errmsg))
    end
end

function check_length_consistency(
        array1::Array, array2::Array
    )

    if length(array1) != length(array2)
        errmsg = (
            "lengths of arrays must be equal"
        )
        throw(ErrorException(errmsg))
    end

end

"""
warn_if_ngrid_is_not_power_of_2 displays warning
if each ngrid in ngrids is not power of 2,
for FFTW efficiency.
"""
function warn_if_ngrid_is_not_power_of_2(ngrids)
    for (axis, ngrid) in enumerate(ngrids)
        if !is_power_of_2(ngrid)
            println(
                "WARNING:",
                " ngrids[", axis, "]=", ngrid,
                " should be power of 2 for FFTW efficiency."
            )
        end
    end
end

function is_power_of_2(num::Int)

    if num<=0
        errmsg="num must be >0"
        throw(DomainError(num, errmsg))
    elseif num%2!=0
        return num==1
    else
        return is_power_of_2(num÷2)
    end

end

max_stars_repo_name: mat-der-D/SpectralMethodFFTW

max_issues_repo_name: mat-der-D/SpectralMethodFFTW

max_forks_repo_name: mat-der-D/SpectralMethodFFTW


text:
using Gnuplot

x = collect(1.:100);

#-----------------------------------------------------------------
for i in 1:10
    @gp :gp1 "plot sin($i*x)"
    @gp :gp2 "plot sin($i*x)"
    @gp :gp3 "plot sin($i*x)"
    sleep(0.3)
end
quitall()

#-----------------------------------------------------------------
@gp "plot sin(x)"
@gp "plot sin(x)" "pl cos(x)"
@gp "plo sin(x)" "s cos(x)"

@gp "plot sin(x)" :-
@gp :- "plot cos(x)"

@gp "plot sin(x)" 2 xr=(-2pi,2pi) "pause 2" "plot cos(4*x)"

x = range(-2pi, stop=2pi, length=100);
y = 1.5 * sin.(0.3 .+ 0.7x) ;
noise = randn(length(x))./2;
e = 0.5 * fill(1, size(x));

@gp hist(noise, nbins=10)

@gp x y
@gp x y "w l"

name = "\$MyDataSet1"
@gp x y name "plot $name w l" "pl $name u 1:(2*\$2) w l"

@gsp randn(Float64, 30, 50)
@gp randn(Float64, 30, 50) "w image"
@gsp x y y

@gp("set key horizontal", "set grid",
    xrange=(-7,7), ylabel="Y label",
    x, y, "w l t 'Real model' dt 2 lw 2 lc rgb 'red'",
    x, y+noise, e, "w errorbars t 'Data'")

@gp "f(x) = a * sin(b + c*x); a = 1; b = 1; c = 1;"   :-
@gp :- x y+noise e name                               :-
@gp :- "fit f(x) $name u 1:2:3 via a, b, c;"          :-
@gp :- "set multiplot layout 2,1"                     :-
@gp :- "plot $name w points" ylab="Data and model"    :-
@gp :- "plot $name u 1:(f(\$1)) w lines"              :-
@gp :- 2 xlab="X label" ylab="Residuals"              :-
@gp :- "plot $name u 1:((f(\$1)-\$2) / \$3):(1) w errorbars notit"

# Retrieve values for a, b and c
a = 0; b = 0; c = 0;
try
    a = parse(Float64, exec("print a"))
    b = parse(Float64, exec("print b"))
    c = parse(Float64, exec("print c"))
catch
end
gnuplot(:dry)
@gp    :dry "f(x) = a * sin(b + c*x); a = 1; b = 1; c = 1;"  :-
@gp :- :dry "a = $a; b = $b; c = $c"                         :-
@gp :- :dry "set multiplot layout 2,1" ylab="Data and model" :-
name = "\$MyDataSet1"
@gp :- :dry x y+noise e name                                 :-
@gp :- :dry "plot $name w points"                            :-
@gp :- :dry "plot $name u 1:(f(\$1)) w lines"                :-
@gp :- :dry 2 xlab="X label" ylab="Residuals"                :-
@gp :- :dry "plot $name u 1:((f(\$1)-\$2) / \$3):(1) w errorbars notit" :-
@gp :- :dry
save("test.gp")        # write on file test.gp
quitall()
exec("load 'test.gp'") # load file test.gp

#-----------------------------------------------------------------
@gp("""
        approx_1(x) = x - x**3/6
        approx_2(x) = x - x**3/6 + x**5/120
        approx_3(x) = x - x**3/6 + x**5/120 - x**7/5040
        label1 = "x - {x^3}/3!"
        label2 = "x - {x^3}/3! + {x^5}/5!"
        label3 = "x - {x^3}/3! + {x^5}/5! - {x^7}/7!"
        #
        set termoption enhanced
        save_encoding = GPVAL_ENCODING
        set encoding utf8
        #
        set title "Polynomial approximation of sin(x)"
        set key Left center top reverse
        set xrange [ -3.2 : 3.2 ]
        set xtics ("-π" -pi, "-π/2" -pi/2, 0, "π/2" pi/2, "π" pi)
        set format y "%.1f"
        set samples 500
        set style fill solid 0.4 noborder""",
    "plot '+' using 1:(sin(\$1)):(approx_1(\$1)) with filledcurve title label1 lt 3",
    "plot '+' using 1:(sin(\$1)):(approx_2(\$1)) with filledcurve title label2 lt 2",
    "plot '+' using 1:(sin(\$1)):(approx_3(\$1)) with filledcurve title label3 lt 1",
    "plot sin(x) with lines lw 1 lc rgb 'black'")

#-----------------------------------------------------------------
@gp("""
        set zrange [-1:1]
        unset label
        unset arrow
        sinc(u,v) = sin(sqrt(u**2+v**2)) / sqrt(u**2+v**2)
        set xrange [-5:5]; set yrange [-5:5]
        set arrow from 5,-5,-1.2 to 5,5,-1.2 lt -1
        set label 1 "increasing v" at 6,0,-1
        set arrow from 5,6,-1 to 5,5,-1 lt -1
        set label 2 "u=0" at 5,6.5,-1
        set arrow from 5,6,sinc(5,5) to 5,5,sinc(5,5) lt -1
        set label 3 "u=1" at 5,6.5,sinc(5,5)
        set parametric
        set hidden3d offset 0	# front/back coloring makes no sense for fenceplot #
        set isosamples 2,33
        xx=-5; dx=(4.99-(-4.99))/9
        x0=xx; xx=xx+dx
        x1=xx; xx=xx+dx
        x2=xx; xx=xx+dx
        x3=xx; xx=xx+dx
        x4=xx; xx=xx+dx
        x5=xx; xx=xx+dx
        x6=xx; xx=xx+dx
        x7=xx; xx=xx+dx
        x8=xx; xx=xx+dx
        x9=xx; xx=xx+dx""",
    "splot [u=0:1][v=-4.99:4.99]x0, v, (u<0.5) ? -1 : sinc(x0,v) notitle",
	"splot x1, v, (u<0.5) ? -1 : sinc(x1,v) notitle",
	"splot x2, v, (u<0.5) ? -1 : sinc(x2,v) notitle",
	"splot x3, v, (u<0.5) ? -1 : sinc(x3,v) notitle",
	"splot x4, v, (u<0.5) ? -1 : sinc(x4,v) notitle",
	"splot x5, v, (u<0.5) ? -1 : sinc(x5,v) notitle",
	"splot x6, v, (u<0.5) ? -1 : sinc(x6,v) notitle",
	"splot x7, v, (u<0.5) ? -1 : sinc(x7,v) notitle",
	"splot x8, v, (u<0.5) ? -1 : sinc(x8,v) notitle",
	"splot x9, v, (u<0.5) ? -1 : sinc(x9,v) notitle")

quitall()


text:
# numerically integrate Mₙ(k²)/(4br)ᵐ
function Mn_num(k2::T, m::Int64) where T
    f(x) = sqrt(k2 - sin(x)^2)^m
    if k2 < 1
        kap2 = T(asin(sqrt(big(k2))))
    else
        kap2 = 0.5 * π
    end
    Mn, err = quadgk(f, -kap2, kap2, rtol=1e-15)
    return Mn
end

function test_Mn(r::T, b::T) where T
    n_max = 30
    u = ones(n_max + 2)
    prec_frac = zeros(n_max + 3)
    prec_abs = zeros(n_max + 3)

    ld = PolynomialLimbDark(u)
    ld_b = PolynomialLimbDark(big.(u))

    t = ld(b, r)
    t_b = ld_b(big(b), big(r))

    onembmr2 = (r + 1 - b) * (1 - r + b)
    k2 = max(0, onembmr2 / (4 * b * r))

    sqbr = sqrt(b * r)

    color = ones(Int, n_max + 3)
    for m in 0:ld.n_max
        Mnn = Mn_num(k2, m) * (2 * sqbr)^m
        Mn = ld.Mn[m + 1]
        Mn_b = Float64(ld_b.Mn[m + 1])
        t1 = isapprox(Mnn, Mn; atol=1e-15, rtol=1e-6)
        t2 = isapprox(Mn, Mn_b; atol=1e-15, rtol=1e-6)
        if !t1
            diff = Mnn - Mn
            @warn "Mn discrepancy" b r k2 Mnn Mn diff
        end
        if !t2
            diff = Mn - Mn_b
            @warn "Mn discrepancy" b r k2 Mn Mn_b diff
        end
        @test t1
        @test t2
        t1 || t2 && (color[begin + m] = 2)
        prec_frac[begin + m] = abs(Mn / Mn_b - 1)
        prec_abs[begin + m] = abs(asinh(Mn) - asinh(Mn_b))
    end
    PLOT && mnerror(prec_frac, prec_abs, c=color, title="r=$r, b=$b")
    PLOT && savefig(plotsdir("Mn_error_r=$(r)_b=$b.png"))
    return nothing
end

@testset "Mn compute" begin
    r = 100.0
    ϵ = 1e-8
    bs = [r - 1 + ϵ, r, r + 1 - ϵ]

    test_Mn.(r, bs)

    r = 0.01
    # points close to boundary
    bs = [ϵ, r, 1 - r - ϵ, 1 - r + ϵ, 1, r + 1 - ϵ]
    test_Mn.(r, bs)

    # random points
    ntest = 10
    for i in 1:ntest
        r = 2.0
        b = 0.25
        while r > 1 + b || r < b - 1
            r = 2 * rand(rng)
            b = 2 * rand(rng)
        end
        test_Mn(r, b)
    end
end

@testset "Mn raise" begin
    nothing
end


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

export
    connectedcomponents, connectedmembership,
    findbiconnected, cutvertices, bridges, edgebiconnectedcomponents

"""
    connectedcomponents(graph::UndirectedGraph) -> Vector{Set{Int}}

Compute connectivity and return sets of the connected components.
"""
function connectedcomponents(graph::UndirectedGraph)
    nodes = nodeset(graph)
    components = Set{Int}[]
    while !isempty(nodes)
        root = pop!(nodes)
        tree = Set(keys(dfstree(adjacencies, graph, root)))
        push!(components, tree)
        setdiff!(nodes, tree)
    end
    return components
end

"""
    connectedmembership(graph::UndirectedGraph) -> Vector{Int}

Return connected component membership array.
"""
function connectedmembership(graph::UndirectedGraph)
    mem = zeros(Int, nodecount(graph))
    for (i, conn) in enumerate(connectedcomponents(graph))
        for n in conn
            mem[n] = i
        end
    end
    return mem
end

struct BiconnectedState{T<:UndirectedGraph}
    graph::T

    pred::Dict{Int,Int}
    level::Dict{Int,Int}
    low::Dict{Int,Int}

    cutvertices::Set{Int}
    bridges::Vector{Int}
    biconnected::Vector{Set{Int}} # biconnected edges

    function BiconnectedState{T}(graph) where {T<:UndirectedGraph}
        new(graph, Dict(), Dict(), Dict(), Set(), [], [])
    end
end
BiconnectedState(graph::UndirectedGraph) = BiconnectedState{typeof(graph)}(graph)

function dfs!(state::BiconnectedState, n::Int)
    state.pred[n] = n
    dfs!(state, 1, n)
end

function dfs!(state::BiconnectedState, depth::Int, n::Int)
    state.level[n] = depth
    state.low[n] = depth
    childcnt = 0
    compbuf = Set{Int}()
    for (ninc, nadj) in neighbors(state.graph, n)
        if !haskey(state.level, nadj)
            @debug "new node: $(nadj)"
            childcnt += 1
            state.pred[nadj] = n
            comp = dfs!(state, depth + 1, nadj)
            push!(comp, ninc)
            if state.low[nadj] >= state.level[n]
                @debug "articulation point: $(n)"
                if state.low[nadj] > state.level[n]
                    @debug "bridge $(ninc)"
                    push!(state.bridges, ninc) # except for bridgehead
                end
                push!(state.biconnected, comp)
                push!(state.cutvertices, n)
            else
                union!(compbuf, comp)
            end
            state.low[n] = min(state.low[n], state.low[nadj])
        elseif state.pred[n] != nadj
            @debug "cycle found: $(n)"
            state.low[n] = min(state.low[n], state.level[nadj])
            push!(compbuf, ninc)
        end
    end
    if depth == 1 && childcnt < 2
        @debug "revert: $(n)"
        delete!(state.cutvertices, n)
    end
    return compbuf
end

function findbiconnected(graph::UndirectedGraph)
    state = BiconnectedState(graph)
    nodes = nodeset(graph)
    while !isempty(nodes)
        dfs!(state, pop!(nodes))
        setdiff!(nodes, keys(state.level))
    end
    return state
end

"""
    cutvertices(graph::UndirectedGraph) -> Set{Int}
    cutvertices(state::BiconnectedState) -> Set{Int}

Compute biconnectivity and return cut vertices (articulation points).
"""
cutvertices(graph::UndirectedGraph) = findbiconnected(graph).cutvertices
cutvertices(state::BiconnectedState) = state.cutvertices

"""
    bridges(graph::UndirectedGraph) -> Set{Int}
    bridges(state::BiconnectedState) -> Set{Int}

Compute biconnectivity and return bridges.
"""
bridges(graph::UndirectedGraph) = findbiconnected(graph).bridges
bridges(state::BiconnectedState) = state.bridges

"""
    edgebiconnectedcomponents(graph::UndirectedGraph) -> Vector{Vector{Int}}

Compute sets of biconnected component edges.
"""
edgebiconnectedcomponents(graph::UndirectedGraph) = findbiconnected(graph).biconnected
edgebiconnectedcomponents(state::BiconnectedState) = state.biconnected


text:
const USE_GPU = true
using ImplicitGlobalGrid
import MPI
using ParallelStencil
using ParallelStencil.FiniteDifferences3D
@static if USE_GPU
    @init_parallel_stencil(CUDA, Float64, 3);
else
    @init_parallel_stencil(Threads, Float64, 3);
end

@parallel function diffusion3D_step!(T2, T, Ci, lam, dt, _dx, _dy, _dz)
    @inn(T2) = @inn(T) + dt*(lam*@inn(Ci)*(@d2_xi(T)*_dx^2 + @d2_yi(T)*_dy^2 + @d2_zi(T)*_dz^2));
    return
end

@parallel function diffusion3D_step!(T2, T, Ci, lam, dt, _dx, _dy, _dz)
    @inn(T2) = @inn(T) + dt*(lam*@inn(Ci)*(@d2_xi(T)*_dx^2 + @d2_yi(T)*_dy^2 + @d2_zi(T)*_dz^2));
    return
end

function diffusion3D()
# Physics
lam        = 1.0;                                        # Thermal conductivity
c0         = 2.0;                                        # Heat capacity
lx, ly, lz = 1.0, 1.0, 1.0;                              # Length of computational domain in dimension x, y and z

# Numerics
nx, ny, nz = 512, 512, 512;                              # Number of gridpoints in dimensions x, y and z
nt         = 100;                                        # Number of time steps
me, dims   = init_global_grid(nx, ny, nz);
dx         = lx/(nx_g()-1);                              # Space step in x-dimension
dy         = ly/(ny_g()-1);                              # Space step in y-dimension
dz         = lz/(nz_g()-1);                              # Space step in z-dimension
_dx, _dy, _dz = 1.0/dx, 1.0/dy, 1.0/dz;

# Array initializations
T   = @zeros(nx, ny, nz);
T2  = @zeros(nx, ny, nz);
Ci  = @zeros(nx, ny, nz);

# Initial conditions
Ci .= 1/c0;                                              # 1/Heat capacity
T  .= 1.7;
T2 .= T;                                                 # Assign also T2 to get correct boundary conditions.

# Time loop
dt   = min(dx^2,dy^2,dz^2)/lam/maximum(Ci)/8.1;          # Time step for 3D Heat diffusion
for it = 1:nt
    if (it == 11) tic(); end                             # Start measuring time.
    @hide_communication (16, 2, 2) begin
        @parallel diffusion3D_step!(T2, T, Ci, lam, dt, _dx, _dy, _dz);
        T2[  1,:,:] .= 1.0;
        T2[end,:,:] .= 1.0;
        update_halo!(T2);
    end
    T, T2 = T2, T;
end
time_s = toc()

# Performance
A_eff = (2*1+1)*1/1e9*nx*ny*nz*sizeof(Data.Number);      # Effective main memory access per iteration [GB] (Lower bound of required memory access: T has to be read and written: 2 whole-array memaccess; Ci has to be read: : 1 whole-array memaccess)
t_it  = time_s/(nt-10);                                  # Execution time per iteration [s]
T_eff = A_eff/t_it;                                      # Effective memory throughput [GB/s]
if (me==0) println("time_s=$time_s T_eff=$T_eff"); end

finalize_global_grid();
end

diffusion3D()


text:
abstract type IntegratorCache{DT,D} end

abstract type ODEIntegratorCache{DT,D} <: IntegratorCache{DT,D} end
abstract type DAEIntegratorCache{DT,D} <: IntegratorCache{DT,D} end
abstract type IODEIntegratorCache{DT,D} <: IntegratorCache{DT,D} end
abstract type IDAEIntegratorCache{DT,D} <: IntegratorCache{DT,D} end
abstract type PODEIntegratorCache{DT,D} <: IntegratorCache{DT,D} end
abstract type PDAEIntegratorCache{DT,D} <: IntegratorCache{DT,D} end

create_integrator_cache(integrator::Integrator) = error("create_integrator_cache()! not implemented for ", typeof(integrator))

initialize!(::Integrator, ::IntegratorCache) = nothing

integrate_step!(integrator::Integrator, ::IntegratorCache) = error("integrate_step()! not implemented for ", typeof(integrator))

function CommonFunctions.reset!(cache::Union{ODEIntegratorCache, DAEIntegratorCache}, Δt)
    cache.t̅  = cache.t
    cache.q̅ .= cache.q
    cache.v̅ .= cache.v
    cache.t += Δt
    cache.n += 1
end

function CommonFunctions.reset!(cache::Union{PODEIntegratorCache, IODEIntegratorCache}, Δt)
    cache.t̅  = cache.t
    cache.q̅ .= cache.q
    cache.p̅ .= cache.p
    cache.v̅ .= cache.v
    cache.f̅ .= cache.f
    cache.t += Δt
    cache.n += 1
end

function CommonFunctions.reset!(cache::Union{PDAEIntegratorCache, IDAEIntegratorCache}, Δt)
    cache.t̅  = cache.t
    cache.q̅ .= cache.q
    cache.p̅ .= cache.p
    cache.λ̅ .= cache.λ
    cache.v̅ .= cache.v
    cache.f̅ .= cache.f
    cache.u̅ .= cache.u
    cache.g̅ .= cache.g
    cache.t += Δt
    cache.n += 1
end

function CommonFunctions.set_solution!(cache::ODEIntegratorCache, sol, n=0)
    t, q = sol
    cache.n  = n
    cache.t  = t
    cache.q .= q
    cache.v .= 0
end

function CommonFunctions.set_solution!(cache::Union{PODEIntegratorCache, IODEIntegratorCache}, sol, n=0)
    t, q, p = sol
    cache.n  = n
    cache.t  = t
    cache.q .= q
    cache.p .= p
    cache.v .= 0
    cache.f .= 0
end

function CommonFunctions.set_solution!(cache::Union{PDAEIntegratorCache, IDAEIntegratorCache}, sol, n=0)
    t, q, p, λ = sol
    cache.n  = n
    cache.t  = t
    cache.q .= q
    cache.p .= p
    cache.λ .= λ
    cache.v .= 0
    cache.f .= 0
end

function CommonFunctions.get_solution(cache::ODEIntegratorCache)
    (cache.t, cache.q)
end

function CommonFunctions.get_solution(cache::DAEIntegratorCache)
    (cache.t, cache.q, cache.λ)
end

function CommonFunctions.get_solution(cache::Union{PODEIntegratorCache, IODEIntegratorCache})
    (cache.t, cache.q, cache.p)
end

function CommonFunctions.get_solution(cache::Union{PDAEIntegratorCache, IDAEIntegratorCache})
    (cache.t, cache.q, cache.p, cache.λ)
end

function copy_solution!(sol::Solution, cache::IntegratorCache, n, m)
    copy_solution!(sol, get_solution(cache)..., n, m)
end

function cut_periodic_solution!(cache::IntegratorCache, periodicity::Vector)
    cut_periodic_solution!(cache.q, periodicity, cache.s̃)
    cache.q .+= cache.s̃
    cache.q̅ .+= cache.s̃
end


text:
module NumericalRange

using LinearAlgebra
using Plots

"""
        nrange(A; nk = 1, thmax = 32, noplot = false)

Numerical range (or field of values) of the matrix `A`.
Evaluates and plots the numerical range of the
`nk` largest leading principal submatrices of `A`, using `thmax`
equally spaced angles in the complex plane.
The defaults are `nk = 1` and `thmax = 32`
The eigenvalues of `A` are plotted as `x`.  The function returns `f` and `e`,
where `f` is the numerical range and `e` is a vector of eigenvalues of `A`.
Setting `noplot = true` suppresses the plot.

This function is a direct translation of `fv.m` in Professor Nick Higham's Matrix Computation Toolbox,
https://www.mathworks.com/matlabcentral/fileexchange/2360-the-matrix-computation-toolbox,
based on an original routine by A. Ruhe.

References:
R. A. Horn and C. R. Johnson, Topics in Matrix Analysis, Cambridge
     University Press, 1991; sec. 1.5.
A. S. Householder, The Theory of Matrices in Numerical Analysis,
     Blaisdell, New York, 1964; sec. 3.3.
C. R. Johnson, Numerical determination of the field of values of a
     general complex matrix, SIAM J. Numer. Anal., 15 (1978),
     pp. 595-602.

"""
function nrange(B; nk=1, thmax=32, noplot=false)
        function rq(A, x)
                return x'*A*x/(x'*x)
        end
        thmax -= 1 # the function uses thmax + 1 angles

        (n, p) = size(B)
        if n != p
                DimensionMismatch("Matrix must be square.")
        end

        z = Matrix{ComplexF64}(undef, 2*thmax + 1, 1)
        F = eigen(B)
        e = F.values
        f = []

        # filter out cases where B is Hermitian or skew-Hermitian
        if B == B'
                f = [minimum(e), maximum(e)]
        elseif B == -B'
                e = imag(e)
                f = [minimum(e), maximum(e)]
                e *= im
                f *= im
        else
                for m = 1:nk
                        ns = n + 1 - m
                        A = B[1:ns, 1:ns]
                        for i = 0:thmax
                                th = i/thmax*pi
                                Ath = exp(im*th)*A
                                H = 0.5*(Ath + Ath')
                                X = eigen(H, sortby = x -> real(x))
                                V = X.vectors
                                z[i + 1] = rq(A, V[:, 1])
                                z[i + 1 + thmax] = rq(A, V[:, ns])
                        end
                        f = [f; z]
                end
                # join up the boundary
                f = [f; f[1, :]]
        end
        if thmax == 0
                f = e
        end

        if !noplot
                p = plot(real(f), imag(f), framestyle=:origin, aspect_ratio = 1, legend=false)
                scatter!(p, real(e), imag(e), markershape=:x, legend=false)
                display(p)
        end

        return f, e
end

end

max_stars_repo_name: ThomasChaffey/NumericalRange


text:
@testset "test for errors" begin
    PR.silence!()

    m = Model(cbc_optimizer)
    @variable(m, -1.0 <= x <= 1.0)
    @variable(m, y)

    # Test for too few points in base partition.
    @test_throws ErrorException construct_univariate_relaxation!(
        m,
        a -> a^3,
        x,
        y,
        [0.0],
        true,
    )

    # Test for unbounded partition point.
    @test_throws ErrorException construct_univariate_relaxation!(
        m,
        a -> a^3,
        x,
        y,
        [0.0, 1e16],
        true,
    )

    # # Test for unbounded function value.
    @test_throws ErrorException construct_univariate_relaxation!(
        m,
        tan,
        x,
        y,
        [0.0, π / 2.0],
        true,
    )

    # # Test for unbounded derivative value.
    @test_throws ErrorException construct_univariate_relaxation!(
        m,
        a -> sqrt(1 - a^2),
        x,
        y,
        [0.0, 1.0],
        true,
    )

    # # Test for unbounded derivative value.
    @test_throws ErrorException construct_univariate_relaxation!(
        m,
        a -> sqrt(1 - a^2),
        x,
        y,
        [0.0, 1.0],
        true,
    )

    # # Test for equal derivative values in adjacent points of the base partition.
    @test_throws ErrorException construct_univariate_relaxation!(
        m,
        sin,
        x,
        y,
        [0.0, 2 * π],
        true,
    )

    # # Test for invalid ordering of points in base partition.
    @test_throws ErrorException construct_univariate_relaxation!(
        m,
        a -> a^3,
        x,
        y,
        [0.0, -1.0, 1.0],
        true,
    )

    @test_throws ErrorException construct_univariate_relaxation!(
        m,
        a -> a^2,
        x,
        y,
        [0.0, 1e-16],
        true,
    )

    x_lb, x_ub = 10 * rand(2) .* [-1, 1]
    y_lb, y_ub = 10 * rand(2) .* [-1, 1]
    x_mid = (x_lb + x_ub) / 2.0
    y_mid = (y_lb + y_ub) / 2.0
    m = JuMP.Model(cbc_optimizer)
    JuMP.@variable(m, x_lb <= x <= x_ub)
    JuMP.@variable(m, y_lb <= y <= y_ub)
    JuMP.@variable(m, z)
    @test_throws ErrorException construct_bilinear_relaxation!(
        m,
        x,
        y,
        z,
        [x_lb, x_mid, x_ub],
        [y_lb, y_mid, y_ub],
    )
end


text:
##############################################################################
#
# The hash trick: use a hash function instead of a lexicon to determine the
# columns of a DocumentTermMatrix-like encoding of the data
#
# To produce a DTM for a Corpus for which we do not have an existing lexicon,
# we need someway to map the terms from each document into column indices.
#
# We use the now standard "Hash Trick" (CITATION NEEDED), in which we
# hash strings and then reduce the resulting integers modulo N, which
# defines the numbers of columns we want our DTM to have. This amounts to
# doing a non-linear dimensionality reduction with low probability that similar
# terms hash to the same dimension.
#
# To make things easier, we wrap Julia's hash functions in a new type,
# TextHashFunction, which maintains information about the desired cardinality
# of the hashes.
#
##############################################################################

mutable struct TextHashFunction
    hash_function::Function
    cardinality::Int
end

"""
```
TextHashFunction(cardinality)
TextHashFunction(hash_function, cardinality)
```

The need to create a lexicon before we can construct a document term matrix is often prohibitive.
We can often employ a trick that has come to be called the Hash Trick in which we replace terms
with their hashed valued using a hash function that outputs integers from 1 to N.

Parameters:
	-  cardinality	    = Max index used for hashing (default 100)
 	-  hash_function    = function used for hashing process (default function present, see code-base)

```julia-repl
julia> h = TextHashFunction(10)
TextHashFunction(hash, 10)
```
"""
TextHashFunction(cardinality::Int) = TextHashFunction(hash, cardinality)

TextHashFunction() = TextHashFunction(hash, 100)

cardinality(h::TextHashFunction) = h.cardinality

"""
```
index_hash(str, TextHashFunc)
```

Shows mapping of string to integer.

Parameters:
	-  str		   = Max index used for hashing (default 100)
 	-  TextHashFunc    = TextHashFunction type object

```julia-repl
julia> h = TextHashFunction(10)
TextHashFunction(hash, 10)

julia> index_hash("a", h)
8

julia> index_hash("b", h)
7
```
"""
function index_hash(s::AbstractString, h::TextHashFunction)
    return Int(rem(h.hash_function(s), h.cardinality)) + 1
end


text:
module DEACfrequency
using Random;
using JLD;
const functionNames = ["exponential_frequency_bins",
                       "fixed_frequency_bins",
                       "load_frequency_bins"]

function exponential_frequency_bins(start::Float64,stop::Float64,length::Int64;epsilon::Float64=1e-5)
    frequency_bins = Array{Float64,1}(undef,length + 1);
    exponential_frequency_bins!(frequency_bins,start,stop,epsilon=epsilon);
    frequency_bins
end

function exponential_frequency_bins!(frequency_bins::Array{Float64,1},start::Float64,stop::Float64;epsilon::Float64=1e-5)
    length = size(frequency_bins,1) - 1;
    frequency_bins[1] = start;
    _logspace = range(log10(start+epsilon), stop=log10(stop), length=length)
    @inbounds for i in 2:length+1
        frequency_bins[i] = exp10(_logspace[i-1]);
    end
    nothing
end

function fixed_frequency_bins(start::Float64,stop::Float64,length::Int64;epsilon::Float64=1e-5)
    frequency_bins = Array{Float64,1}(undef,length+1);
    fixed_frequency_bins!(frequency_bins,start,stop);
    frequency_bins
end

function fixed_frequency_bins!(frequency_bins::Array{Float64,1},start::Float64,stop::Float64;epsilon::Float64=1e-5)
    length = size(frequency_bins,1);
    _linspace = range(start,stop=stop,length=length);
    for i in 1:length
        frequency_bins[i]=_linspace[i];
    end
    nothing
end

function load_frequency_bins(start::Float64,stop::Float64,length::Int64;filename::String="./frequency.jld")
    frequency_bins = load_frequency_bins(filename);
    frequency_bins
end

function load_frequency_bins(filename::String="./frequency.jld")
    frequency_bins = JLD.load(filename,"frequency_bins");
    frequency_bins
end

end


text:
tol = 1e-13

pmax = 3

for ii = 1:pmax

  (ploc,) = luteos.genlocal3D( ii )

  (ϕ,∇ϕ) = luteos.basisFuncTetLag( Val{ii}, ploc[:,2], ploc[:,3], ploc[:,4] )
  println( "      3D")

  sz  = size( ϕ )
  for jj = 1:sz[1], kk = 1:sz[2]
    if jj == kk
      @test abs(ϕ[jj,kk] - 1.0) < tol
    else
      @test abs(ϕ[jj,kk]) < tol
    end
  end

  # check evaluation for random point
  xloc   = [0.11698  0.2678  0.05263;
            0.06905  0.1376  0.08579]
  (ϕ,∇ϕ) = luteos.basisFuncTetLag( Val{ii}, xloc[:,1], xloc[:,2], xloc[:,3] )

  for jj in 1:2, dd in 1:3
    xloctest = ϕ[:,jj]' * ploc[:,dd+1]
    @test abs( xloctest - xloc[jj,dd] ) < tol

  end

end


text:
function test_format(fmt, before, after)
    sandbox_dir = joinpath(@__DIR__, "tmp")
    mkdir(sandbox_dir)
    try
        code_path = joinpath(sandbox_dir, "code.jl")
        open(io -> write(io, before), code_path, "w")

        @test fmt(code_path) == false
        @test strip(read(code_path, String)) == after
    finally
        rm(sandbox_dir; recursive = true)
    end
end

@testset "$f interface" for f in (format, format_file)
    before = "foo(; k =v)"
    @testset "style can be positional arg or keyword" begin
        pos = path -> f(path, BlueStyle())
        key = path -> f(path; style = BlueStyle())
        after = "foo(; k=v)"
        test_format(pos, before, after)
        test_format(key, before, after)
    end
    @testset "other keywords take affect" begin
        pos = path -> f(path, BlueStyle(); whitespace_in_kwargs = true)
        key = path -> f(path; style = BlueStyle(), whitespace_in_kwargs = true)
        after = "foo(; k = v)"
        test_format(pos, before, after)
        test_format(key, before, after)
    end
end


text:
using Libdl

const DEPS_FILE = joinpath(@__DIR__, "deps.jl")

if isfile(DEPS_FILE)
    rm(DEPS_FILE)
end

function write_depsfile(path)
    open(DEPS_FILE, "w") do io
        println(io, "const libgurobi = \"$(escape_string(path))\"")
    end
end

const ALIASES = [
    "gurobi90",
    "gurobi91",
]

paths_to_try = copy(ALIASES)

for a in ALIASES
    if haskey(ENV, "GUROBI_HOME")
        if Sys.isunix()
            push!(paths_to_try, joinpath(ENV["GUROBI_HOME"], "lib", string("lib", a, ".so")))
        end
        if Sys.iswindows()
            push!(paths_to_try, joinpath(ENV["GUROBI_HOME"], "bin", string(a, ".", Libdl.dlext)))
        end
        if Sys.isapple()
            push!(paths_to_try, joinpath(ENV["GUROBI_HOME"], "lib", string("lib", a, ".dylib")))
        end
    end
    # gurobi uses .so on OS X for some reason
    if Sys.isapple()
        push!(paths_to_try, string("lib$a.so"))
        push!(paths_to_try, string("lib$a.dylib"))
    end
end

found = false
for l in paths_to_try
    d = Libdl.dlopen_e(l)
    if d != C_NULL
        global found = true
        write_depsfile(l)
        break
    end
end

function _print_GUROBI_HOME_help()
    println("""
    You should set the `GUROBI_HOME` environment variable to point to the
    install location then try again. For example (updating the path to the
    correct location if needed):
    ```
    # On Windows, this might be
    ENV["GUROBI_HOME"] = "C:\\\\Program Files\\\\gurobi910\\\\win64\\\\"
    import Pkg
    Pkg.add("Gurobi")
    Pkg.build("Gurobi")

    # On OSX, this might be
    ENV["GUROBI_HOME"] = "/Library/gurobi910/mac64/"
    import Pkg
    Pkg.add("Gurobi")
    Pkg.build("Gurobi")

    # On Unix, this might be
    ENV["GUROBI_HOME"] = "/opt/gurobi910/linux64/"
    import Pkg
    Pkg.add("Gurobi")
    Pkg.build("Gurobi")
    ```
    **Note: your path may differ. Check which folder you installed the Gurobi
    binary in, and update the path accordingly.**
    """)
end

function diagnose_gurobi_install()
    println("""

    **Unable to locate Gurobi installation. Running some common diagnostics.**

    Gurobi.jl only supports the following versions:
    """)
    println.(" - ", ALIASES)
    println("""

    Did you download and install one of these versions from gurobi.com?
    Installing Gurobi.jl via the Julia package manager is _not_ sufficient!
    """)
    if haskey(ENV, "GUROBI_HOME")
        dir = joinpath(ENV["GUROBI_HOME"], Sys.isunix() ? "lib" : "bin")
        println("""
        Found GUROBI_HOME =  $(ENV["GUROBI_HOME"])

        Does this point to the correct install location?

        We're going to look for the Gurobi library in this directory:
            $(dir)

        That directory has the following files:
        """)
        try
            for file in readdir(dir)
                println(" - ", joinpath(dir, file))
            end
            println("""

            We were looking for (but could not find) a file named like
            `libgurobiXXX.so`, `libgurobiXXX.dylib`, or `gurobiXXX.dll`.\n\n""")

            _print_GUROBI_HOME_help()
        catch ex
            if typeof(ex) <: SystemError
                println("""
                Aha! We tried looking in `$(dir)`, but something went wrong. Are
                you sure that your GUROBI_HOME environment variable is correct?
                When combined with the appropriate suffix (e.g., `lib` or
                `bin`, it needs to point to a valid directory.\n\n""")
                _print_GUROBI_HOME_help()
            else
                rethrow(ex)
            end
        end
    else
        try
            # Try to call `gurobi_cl`. This should work if Gurobi is on the
            # system path. If it succeeds, it will print out the version.
            io = IOBuffer()
            run(pipeline(`gurobi_cl --version`; stdout = io))
            seekstart(io)
            println("""
            We couldn't find the `GUROBI_HOME` environment variable, but we
            found this version of Gurobi on your `PATH`.

            $(read(io, String))
            Is this version one of the supported versions listed above? If so,
            we found the executable, but not the libraries we need. Follow the
            advice below to set the `GUROBI_HOME` environment variable. If not,
            you should edit your `PATH` to point to the correct version, or set
            the `GUROBI_HOME` environment variable.\n""")
            _print_GUROBI_HOME_help()
        catch
            println("""

            We could not find a version of Gurobi in your `PATH`, and we could
            not find the environment variable `GUROBI_HOME`.\n\n""")
            _print_GUROBI_HOME_help()
        end
    end
end

if haskey(ENV, "GUROBI_JL_SKIP_LIB_CHECK")
    # Skip!
elseif get(ENV, "JULIA_REGISTRYCI_AUTOMERGE", "false") == "true"
    write_depsfile("julia_registryci_automerge")
elseif !found
    diagnose_gurobi_install()
    error("""
    Unable to locate Gurobi installation. If the advice above did not help,
    open an issue at https://github.com/jump-dev/Gurobi.jl and post the full
    print-out of this diagnostic attempt.
    """)
end


text:
import NearestNeighbors: KDTree, BallTree, BruteTree

KDTree(E::AbstractEmbedding,
        metric = Euclidean(); kwargs...) =
    NearestNeighbors.KDTree(E.points, metric; kwargs...)

BallTree(E::AbstractEmbedding,
        metric = Euclidean(); kwargs...) =
    NearestNeighbors.BallTree(E.points, metric; kwargs...)

BruteTree(E::AbstractEmbedding,
        metric = Euclidean(); kwargs...) =
    NearestNeighbors.BruteTree(E.points, metric; kwargs...)

export KDTree, BallTree, BruteTree


text:
module TimedLoop

using Dates
using Printf

export @tml, @tmlv

######## time to h:m:s
function hms(dt)
    (h, r) = divrem(dt, 60 * 60)
    (m, r) = divrem(r, 60)
    string(Int(h), ":", Int(m), ":", round(Int, r))
end

# ref: threaded for loop Threads.@threads
# https://github.com/JuliaLang/julia/blob/46ce4d79337bdd257ee2e3d2f4bb1c55ff0a5030/base/threadingconstructs.jl

function m_tml(iter, lbody, v_out)
    lidx = iter.args[1]         # index
    range = iter.args[2]

    return quote
        vout = $(esc(v_out))
        range_ = $(esc(range))
        n_loop = length(range_)
        local timing_total = 0.0
        println(">>>")
        for i in 1:n_loop
            timing_now = now()
            ## loop body
            local $(esc(lidx)) = range_[i]
            $(esc(lbody))
            ## end of loop body
            ##### timing logic
            timing_trj = (now() - timing_now).value / 1000
            timing_total += timing_trj
            timing_avg = timing_total / i

            should_print = vout || i < 10 || mod(i,10) == 0 || timing_avg > 5.0

            if should_print
                timing_tot_est = timing_avg * n_loop
                timing_rem = timing_tot_est - timing_total
                timing_eta = Dates.format(
                    now() + Dates.Millisecond(round(Int64, 1000 * timing_rem)),
                    "e, dd u yyyy HH:MM:SS",
                )
                @printf("%i/%i in %.3f sec.\n", i, n_loop, timing_trj)
                @printf("ETA: %s\ntotal run est. %s, %.3f s per loop, %s to go\n", timing_eta, hms(timing_tot_est), timing_avg, hms(timing_rem))
                println("="^7)
            end #if
        end #for
    end #quote
end

function construct_tml(args...)
    na = length(args)

    if !(na in (1,2))
        throw(ArgumentError("wrong number of arguments"))
    end

    if na == 1
        ex = args[1]
        v_out = false
    else
        ex = args[2]
        v_out = args[1]
        if !isa(v_out, Bool)
            throw(ArgumentError("verbose switch should be Boolean"))
        end
    end

    if !isa(ex, Expr)
        throw(ArgumentError("need an expression argument"))
    end

    if ex.head === :for
        iter = ex.args[1]
        lbody = ex.args[2]
        m_tml(iter, lbody, v_out)
    else
        throw(ArgumentError("unrecognized arguments"))
    end
end

macro tml(args...)
    construct_tml(args...)
end

macro tmlv(args...)
    construct_tml(true, args...)
end

end
#mod

max_stars_repo_name: daviehh/TimedLoop

max_issues_repo_name: daviehh/TimedLoop

max_forks_repo_name: daviehh/TimedLoop