text:
#########################################################################################################################################
######### knotvector.jl contains all basic functionality in order to construct and manipulate B-spline knotvectors ######################
#########################################################################################################################################
#
# Functions:
#       buildvector()           -  construct a knotvector using a vector with increasing values of type T and a vector #                                          specifying the knot multiplicity
#       uniqueknots()           -  decompose a knotvector into an increasing vector of Float64, the knotmultiplicity, and vectors specifying indexing into the knotvector
#       globalrefine()          -  put k new knots inside the non-zero knotspans of the knotvector
#       findspan()              -  determine the knot span index of a point 'u' in the knotvector 'kts'
#       findnonzerospan()       -  find span of all non-zero span elements in knotvector kts
#       getelementknotvectors() -  determine the local knotvectors of each non-zero-span element in the knovector
#
#########################################################################################################################################
#
# File part of the Julia IsoGeometric Analysis toolbox
# R.R.Hiemstra
# 14-09-2015
#
########################################################################################################################################

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

# typealiases
const Degree = Int
const KnotVector = Vector{Float64}

"""
construct a knotvector using a vector with increasing values of type T and a vector specifying the knot multiplicity

## input:
    uvals :: Vector{T}   - vector with increasing values of type T
    mult  :: Vector{Int} - vector specifying the knot multiplicity

## output:
    kts   :: Vector{T}   - knotvector with knot values of type T
"""
function buildvector(uvals::Array{T,1},mult::Array{Int64,1}) where {T<:Number}
  # allocate space
  v = zeros(T,sum(mult))

  # set values
  ind = 1
  for i in 1:length(uvals)
      v[ind:ind+mult[i]-1] = fill(uvals[i],mult[i])
      ind = ind + mult[i]
  end

  return v
end

"""
decompose a knotvector into an increasing vector of Float64, the knotmultiplicity, and vectors specifying indexing into the knotvector

## input:
    kts :: Vector{Float64}  - knotvector

## output:
    ukts   :: Vector{Float64} - unique knot values in increasing order
    mult   :: Vector{Int}     - corresponding knot multiplicity
    inda   :: Vector{Int}     - indexing such that kts==ukts[inda]
    indb   :: Vector{Int}     - indexing such that kts[indb]==ukts
"""
function uniqueknots(knots::Array{Float64,1})
    n = length(knots)
    mult = ones(Int64,n)
    uknots = zeros(Float64,n)
    inda = zeros(Int64,n)
    indc = zeros(Int64,n)

    j = 1; uknots[1] = knots[1]; inda[1] = 1; indc[1] = 1
    for i in 2:n
        if knots[i]==knots[i-1]
            mult[j] += 1
        else
            j+=1
            uknots[j] = knots[i]
            indc[j] = i
        end
        inda[i] = j
    end
    return uknots[1:j], mult[1:j], inda, indc[1:j]
end

function uniqueknots(kts)
  n = length(kts)
  temp = ntuple(i -> uniqueknots(kts[i]), n)
  return ukts, umult, inda, indc = ntuple(j -> ntuple(i -> temp[i][j], n), 4)
end

"""
global refine knotvector: put k new knots in non-zero knotspans

## input:
    kts  :: Vector{Float64}    - knotvector
    k    :: Int                - integer prescribing the number of new knots to put in current non-zero knot-spans

## output:
    newkts :: Vector{Float64}  - refined knotvector
"""
function globalrefine(kts::Vector{Float64},k::Int)
    # get unique knots and multiplicity
    ukts, umult, inda, indc = uniqueknots(kts)
    n = length(kts); m = length(ukts)

    # compute new knotvector
    newkts = zeros(n + k*(m-1)); j = 1
    for i in 1:m-1
        newkts[j:j+umult[i]-1] = ukts[i]
        j += umult[i]-1
        newkts[j:j+k+1] = linspace(ukts[i],ukts[i+1],k+2)
        j += k+1
    end
    newkts[j:j+umult[m]-1] = ukts[m]
    return newkts
end

"""
Determine the knot span index of a point 'u' in th knotvector 'kts'

## input:
    p      :: Int              - polynomial degree
    kts    :: Vector{Float64}  - knotvector
    u      :: Float64          - evaluation point

## output:
    span   :: Int              - knot span index of a point 'u' in th knotvector 'kts'
"""
function findspan(p::Integer, kts::Vector{Float64}, u::Float64)
    if u<kts[1] || u>kts[end]
        error("Boundserror findspan(p::Integer, kts::Vector{Float64}, u::Float64)")
    end

    if u==kts[end]
        span = length(kts)-1
        while kts[span]==kts[span+1]
            span+=-1
        end
    else
        low = 0
        high = length(kts)
        mid = round(Int64,(low + high)/2)

        while u < kts[mid] || u >= kts[mid + 1]
          if u < kts[mid]
            high = mid
          else
            low = mid
          end
          mid = round(Int64,(low + high)/2)
        end
        span = round(Int64,(low + high)/2)
    end
    return span
end

"""
Determine the knot span index of a vector of points 'u' in the knotvector 'kts'

## input:
    p      :: Int              - polynomial degree
    kts    :: Vector{Float64}  - knotvector
    u      :: Vector{Float64}  - evaluation point

## output:
    span   :: Vector{Int}      - vector with knot span index of each point u[i] in the knotvector 'kts'
"""
function findspan(p::Integer, kts::Vector{Float64}, u::Vector{Float64})
    # spline space dimension
    n = length(kts)-p-1

    # allocate space for span
    m = length(u)
    span = zeros(Int64,m)
    span[1] = findspan(p,kts,u[1])
    span[m] = findspan(p,kts,u[m])

    # find span at u[i]
    for i in 2:m-1
        span[i] = span[i-1]
        while u[i] > kts[span[i]+1]
            span[i] += 1
        end
    end

    return span
end

function findfunswithcompactsupp(p::Int,kts::Vector{Float64},a::Float64,b::Float64)
    n = length(kts); k = 1
    while a!=kts[k] && k<n
        k+=1
    end
    a = k
    while b>=kts[k] && k<=n
        k+=1
    end
    b = k-2-p
    return a:b
end

"""
find span of all non-zero span elements in knotvector 'kts'

## input:
    p      :: Int              - polynomial degree
    kts    :: Vector{Float64}  - knotvector

## output:
    span   :: Vector{Int}      - vector with knot span index of each non-zero-span elements in the knotvector 'kts'
"""
function findnonzerospan(p::Int,kts::Vector{Float64})
    m = length(kts)
    span = zeros(Int64,m)

    count = 1
    for j in p+1:m-p-1
        if kts[j]!=kts[j+1]
            span[count] = j
            count+=1
        end
    end

    return span[1:count-1]
end
findnonzerospan(kts::NTuple{N,Vector{Float64}}) where {N} = ntuple(length(kts),(j)->findnonzerospan(kts[j]))

"find the elements in support of basisfunctions"
function findsupp(p::Int, kts::Vector{Float64})
    span = findnonzerospan(p,kts)
    n    = dimsplinespace(p,kts)
    m    = length(span)

    graph = zeros(Bool,m,n)
    for k in 1:m
        graph[k,span[k]-p:span[k]] = true
    end

    return [[1:m][graph[:,k]] for k in 1:n]
end

"""
get element knot vectors from global knot partition

## input:
    p      :: Int              - polynomial degree
    kts    :: Vector{Float64}  - knotvector
    span   :: Vector{Int}      - span-index of each non-zero knot-span in the knotvector

## output:
    ekts   :: Vector{Vector{Float64}} - vector knotvectors local to each element in knotvector kts
"""
getelementknotvectors(p::Int, kts::Vector{Float64}, span::Vector{Int}) = [kts[span[i]-p:span[i]+1+p] for i in 1:length(span)]

"""
get element knot vectors from global knot partition

## input:
    p      :: Int              - polynomial degree
    kts    :: Vector{Float64}  - knotvector

## output:
    vkts   :: Vector{Vector{Float64}} - vector knotvectors local to each element in knotvector kts
"""
getelementknotvectors(p::Int, kts::Vector{Float64}) = getelementknotvectors(p,kts,findnonzerospan(p,kts))

max_stars_repo_name: michakraus/iga-demos

max_issues_repo_name: michakraus/iga-demos

max_forks_repo_name: michakraus/iga-demos


text:
# import CSV
import FileIO

ext: jl


text:
using Munkres
using Base.Test
using Combinatorics # Required for permutations() on 0.5

"""
Solve assignment problem via brute force exploration over all permutations
"""
function brute_force_optima(cost_matrix)
    (n,m) = size(cost_matrix)
    if n != m
        error("Non-square cost matrix, I haven't gotten around to doing that yet!")
    end
    initial_permuation = 1:n
    best_permuation = 1:n
    min_cost = typemax(eltype(cost_matrix))
    for p in permutations(initial_permuation)
        current_cost = zero(eltype(cost_matrix))
        for i = 1:n
            current_cost += cost_matrix[initial_permuation[i],p[i]]
        end
        if current_cost <= min_cost
            min_cost = current_cost
            best_permuation = p
        end
    end
    return best_permuation
end

# Small matrix designed to hit all steps
tst = [1 2 3;
       2 4 6;
       3 6 9]
@test munkres(tst) == brute_force_optima(tst)

for i=1:10
    tst = rand(8,8)
    @test munkres(tst) == brute_force_optima(tst)
end

for i=1:10
    # Exponentially distributed test for pathological floating point truncation
    tst = 10.^(50*rand(8,8))
    @test munkres(tst) == brute_force_optima(map(BigFloat, tst))
end

tst = -eye(10)
@test munkres(tst) == brute_force_optima(tst)

#test non-square behaviour on 1 x n and n x 1 arrays

for i=1:10
    tst = rand(1,10)
    @test munkres(tst)[1] == indmin(tst)
end

for i=1:10
    tst = rand(10,1)
    @test findfirst(munkres(tst).==1) == indmin(tst)
end

#test against solutions from Yi Cao's matlab code
#(see http://mathworks.com/matlabcentral/fileexchange/20328-munkres-assignment-algorithm)

tst = [ 678.0 24.0 601.0 34.0 747.0 455.0 871.0 714.0 64.0 547.0
        729.0 1.0 160.0 183.0 695.0 161.0 38.0 790.0 120.0 251.0
        983.0 472.0 960.0 63.0 994.0 761.0 375.0 287.0 959.0 698.0
        681.0 307.0 627.0 6.0 265.0 296.0 605.0 696.0 977.0 283.0
        685.0 203.0 811.0 665.0 422.0 430.0 202.0 317.0 879.0 934.0
        372.0 481.0 839.0 569.0 869.0 898.0 10.0 50.0 376.0 277.0
        395.0 634.0 85.0 156.0 219.0 590.0 924.0 194.0 434.0 294.0
        126.0 132.0 505.0 119.0 637.0 903.0 593.0 487.0 816.0 978.0
        48.0 665.0 766.0 707.0 834.0 291.0 853.0 813.0 504.0 716.0
        702.0 507.0 562.0 685.0 628.0 767.0 708.0 124.0 962.0 442.0
        827.0 71.0 946.0 672.0 117.0 47.0 592.0 876.0 614.0 115.0
        881.0 807.0 375.0 748.0 417.0 958.0 898.0 830.0 448.0 740.0
        207.0 681.0 506.0 159.0 294.0 547.0 109.0 876.0 802.0 231.0
        307.0 125.0 837.0 460.0 467.0 419.0 577.0 687.0 904.0 445.0
        757.0 190.0 26.0 601.0 49.0 71.0 814.0 476.0 268.0 998.0
        87.0 858.0 870.0 73.0 11.0 855.0 51.0 750.0 317.0 942.0
        747.0 687.0 424.0 630.0 175.0 446.0 371.0 14.0 902.0 406.0
        307.0 893.0 131.0 664.0 689.0 482.0 607.0 369.0 283.0 98.0
        442.0 59.0 979.0 463.0 161.0 240.0 328.0 949.0 409.0 299.0
        116.0 313.0 820.0 929.0 347.0 958.0 918.0 806.0 999.0 400.0]

p = [9, 2, 0, 4, 0, 7, 0, 0, 1, 0, 6, 0, 0, 0, 3, 5, 8, 10, 0, 0]

@test munkres(tst) == p

tst = [351.0 300.0 662.0 337.0 384.0 935.0 650.0 446.0 57.0 295.0 168.0 398.0 810.0 561.0 458.0 395.0 188.0 352.0 297.0 602.0
 221.0 662.0 123.0 448.0 255.0 124.0 724.0 935.0 422.0 340.0 335.0 945.0 368.0 423.0 347.0 809.0 823.0 769.0 879.0 886.0
 746.0 625.0 661.0 40.0 541.0 704.0 376.0 324.0 406.0 979.0 327.0 98.0 163.0 742.0 975.0 758.0 646.0 19.0 461.0 836.0
 791.0 962.0 62.0 123.0 813.0 394.0 364.0 297.0 266.0 192.0 863.0 298.0 569.0 462.0 698.0 876.0 317.0 227.0 400.0 869.0
 198.0 145.0 512.0 26.0 352.0 813.0 138.0 43.0 309.0 58.0 868.0 345.0 665.0 538.0 874.0 198.0 118.0 537.0 706.0 787.0
 370.0 975.0 516.0 944.0 310.0 175.0 11.0 895.0 932.0 818.0 264.0 754.0 727.0 938.0 844.0 451.0 11.0 828.0 974.0 654.0
 300.0 506.0 732.0 960.0 303.0 286.0 826.0 454.0 284.0 399.0 179.0 842.0 396.0 966.0 979.0 151.0 548.0 386.0 609.0 480.0
 479.0 209.0 945.0 947.0 493.0 906.0 97.0 488.0 204.0 817.0 681.0 201.0 583.0 740.0 218.0 237.0 662.0 417.0 814.0 672.0
 170.0 471.0 203.0 371.0 218.0 310.0 52.0 882.0 742.0 553.0 704.0 366.0 16.0 774.0 715.0 951.0 810.0 10.0 755.0 455.0
 536.0 155.0 664.0 400.0 807.0 385.0 837.0 439.0 804.0 885.0 371.0 929.0 118.0 210.0 553.0 104.0 367.0 767.0 255.0 50.0]

p = [9, 6, 18, 3, 4, 17, 16, 7, 13, 20]

@test munkres(tst) == p


text:
using Gilbert
using Test

@testset "size $m,$n" for m = 1:20, n = 1:20
    list = gilbertindices((m,n))
    @test length(list) == m*n
    @test list[1] == CartesianIndex(1,1)
    if m >= n
        if !(n == 2 && isodd(m))
            @test list[end] == CartesianIndex(m,1)
        end
    else
        if !(m == 2 && isodd(n))
            @test list[end] == CartesianIndex(1,n)
        end
    end

    ndiag = 0
    for i = 1:m*n-1
        Δ = map(abs, (list[i+1] - list[i]).I)
        @test Δ[1] <= 1
        @test Δ[2] <= 1
        if Δ[1] + Δ[2] > 1
            ndiag += 1
        end
    end
    if m > n && isodd(m) && iseven(n) || m < n && isodd(n) && iseven(m)
        @test ndiag <= 1
    else
        @test ndiag == 0
    end

    L = Gilbert.linearindices(list)
    @test !any(iszero, L)
    @test sum(L) == sum(1:m*n)
end


text:
using TFIMSampler

N = 10^4
h = 1.0
# Ls = [4, 8, 16]
# Ls = [32, 64, 128]
Ls = [512]#, 1024]
# 512 takes about 10 hours for 10^4 samples (5h warm up and 5h run)
# 1024 takes about 4 or 5 days
repeats = [i for i in 1:100]

arg = parse(Int, ARGS[1])
list = Iterators.product(Ls, repeats) |> collect
println("length of array jobs is:\t", length(list))

L, r = list[arg]
println("\nL is $L")
println("h is $h")
println("r is $r")
println("N is $N")

# pre-compile once
@time configs = TFIMSampler.single_sample(L=2, h=2.0, N=1)

# make directory if none exists
h = round(h, digits=2)
dir = "/scratch/mbeach/TFIM_samples/PBC/h-$h/L-$L/"
mkpath(dir)
file = dir * "run-$r"

@time configs = TFIMSampler.single_sample(L=L, h=h, N=N, file=file)


text:
"""
***
# Viscosity of Gases
***
Data are calculated by solving the equation:

``\\mu_{G}= A+B*T+C*T^2``

This function can be used by two methods:
* Get the gas viscosity by introducing the component's name and the temperature.
* Get the full list of components' names supported.
***
### Example:
`julia> GasViscosity("Water",298.15)`

`89.64027655549998`

`julia> GasViscosity()`

`String["Carbon monoxide", "1,3-Butadiene", "Cyclohexanol", "Styrene", "1-Nonane", "1-Butanol", "n-Octane", "Acetic acid", "Hydrogen", "Dichloromethane", "1-Octene", "Methanol", "Methyl chloride", "Ethylene glycol", "Acrylic acid", "Ethanol", "1,1,1-Trichloroethane", "Oxygen", "1,4-Dioxane", "n-Nonane", "Pyridine", "Carbon disulfide", "Hydrogen cyanide", "Ammonia", "Hydrogen chloride", "1-Decene", "1-Propanol", "Acetic anhydride", "Chlorine", "1,1-Dichloroethane", "Methyl isobutyl ketone", "Chloroform", "1,1,2-Trichloroethane", "p-Xylene", "Butyric acid", "Nitrobenzene", "i-Butane", "1-Butene", "Ethylene oxide", "Cyclohexane", "Acetaldehyde", "1,2-Dichloroethane", "Methane", "n-Heptane", "n-Hexane", "Nitrogen", "Carbon tetrachloride", "Hydrogen peroxide", "n-Decane", "Ethane", "Water", "Diethyl ketone", "Ethyl chloride", "Naphthalene", "Fluorine", "o-Xylene", "n-Butane", "Toluene", "n-Pentane", "Aniline", "Methyl bromide", "m-Xylene", "Carbon dioxide", "Benzene", "Phenol", "Propane", "Methylamine", "Dimethylamine", "Acetone", "Ethyl acetate", "Ethylbenzene", "Cyclopropane", "Glycerol", "Sulfur dioxide", "Ethylene", "Diethyl ether", "Vinyl acetate"]`
"""
function GasViscosity(substance::String,T::Float64)
    try
        param=DataGasViscosity[substance]

        if T >= param[5] || T <= param[4]
            print_with_color(:red,
            @sprintf("%4.4f K is out of range %4.4f K - %4.4f K.\n",T,
            param[4],param[5]))
            return NaN
        else
            ρG=param[1]+param[2]*T+param[3]*T^2
        end
    catch
        print_with_color(:red,@sprintf("%s not found in database.\n",substance))
        return NaN
    end
end

function GasViscosity()
    println(keys(DataGasViscosity))
end


text:
export ipqp

"""
	function stepLength(x,y,dx,dy;eta)

	compute step length such that x+alpha*dx and y+alpha*dy are feasible

	then (optional) backtrack a little by multiplying with eta
"""
function stepLength{T<:AbstractFloat}(x::Vector{T},y::Vector{T},dx::Vector{T},dy::Vector{T};eta::T=0.999)
    alphax = -1./min(minimum(dx./x),-1)
    alphax = min(1,alphax)

    alphay = -1./min(minimum(dy./y),-1)
    alphay = min(1,alphay)

    return eta*min(alphax,alphay)
end

"""
	function splitVecIPQP(v,n,me,mi)

	get parts of v=[x;y;le;li]. used in ipqp.
"""
function splitVecIPQP(v::Vector{Float64},n::Int64,me::Int64,mi::Int64)
    x   = v[1:n]
    y   = v[n+(1:mi)]
    le  = v[n+mi+(1:me)]
	li  = v[n+mi+me+1:end]
    return x,y,le,li
end

"""
	function getF(G,c,Ae,be,Ai,bi,x,y,le,li;tau::Real=0.0)

	gets right hand side of KKT system in interior point method

"""
function getF{T}(G::AbstractArray{T},c::Vector{T},Ae::AbstractArray{T},be::Vector{T},Ai::AbstractArray{T},bi::Vector{T},
	             x::Vector{T},y::Vector{T},le::Vector{T},li::Vector{T};
	             dy::Vector{T}=zeros(T,length(y)),dli::Vector{T}=zeros(T,length(li)),tau::T=0.0)
    rd    = G*x - Ae'*le - Ai'*li + c
    rpe   = Ae*x - be
    rpi   = Ai*x - y - bi
    rli   = y.*li + dy.*dli - tau
    return rd,rpe,rpi,rli,y.*li
end

"""
	function getJ(G,Ae,Ai,y,li)

	gets KKT matrix for interior point method for solving

	      | G   0   -Ae'  -Ai'  |
	      | Ae  0    0     0    |
	KKT = | Ai -I    0     0    |
	      | 0   Li   0     Y    |

	where Li=sdiag(li) and Y=sdiagg(y) are diagonal matrices

		min_x .5*x'*G*c + c'*x  subject to Ae*x = be and Ai*x >= bi

"""
function getJ{T}(G::AbstractArray{T},Ae::AbstractArray{T},Ai::AbstractArray{T},y::Vector{T},li::Vector{T})
    n  = size(G,1)
    me = size(Ae,1)
    mi = size(Ai,1)

    In    = speye(n)
    Imi   = speye(mi)
    Znmi  = spzeros(n,mi)
	Znme  = spzeros(n,me)
	Znn   = spzeros(n,n)
    Zmimi = spzeros(mi,mi)
	Zmeme = spzeros(me,me)
	Zmime = spzeros(mi,me)

    Y   = sdiag(y)
    Li  = sdiag(li)

    return [G Znmi -Ae' -Ai'; Ae Zmime' Zmeme Zmime'; Ai -Imi Zmime Zmimi; Znmi' Li Zmime Y]
end

"""
	function solveNormalEqKKT(G,c,Ae,be,Ai,bi,x,y,le,li,rd,rpe,rpi,rli;KKTinv=[],doClear=true)

	brings KKT system to \"normal equation\" form and solves it.

"""
function solveNormalEqKKT{T}(G::AbstractArray{T},Ae::AbstractArray{T},Ai::AbstractArray{T},                  y::Vector{T},li::Vector{T},rd::Vector{T},rpe::Vector{T},rpi::Vector{T},rli::Vector{T};KKTinv=cholfact(speye(1)),doClear=true)

	# build normal eq matrix
    n  = size(G,1)
    me = size(Ae,1)
    mi = size(Ai,1)

 	liy = sdiag(li./y)
	yLi = sdiag(y./li)

    rh1  = -(rd+Ai'*liy*(rpi+(rli./li)));
	rh2  = -rpe

	if size(KKTinv,1)!=size(rh1,1) || doClear
		try
			KKTinv = cholfact(G+Ai'*liy*Ai)
		catch E
			if isa(E,Base.LinAlg.PosDefException)
				KKTinv = pinv(G+Ai'*liy*Ai)
			else
				throw(E)
			end
		end
	end

	if isa(KKTinv,Array)
		if me==0
			dle = zeros(0)
			dx  = KKTinv*rh1
		else
			dle = (Ae*(KKTinv*Ae'))\(rh2-Ae*(KKTinv*rh1))
			dx  = KKTinv*(rh1+Ae'*dle)::Vector{T}
		end
	else
		if me==0
			dle =zeros(0)
			dx  = KKTinv\rh1
		else
			dle =  (Ae*(KKTinv\Ae'))\(rh2-Ae*(KKTinv\rh1))
			dx  = KKTinv\(rh1+Ae'*dle)::Vector{T}
		end
	end
	dli  = yLi\(-rpi-(rli./li)-Ai*dx)
	dy   = (-rli-y.*dli)./li
	return vec([dx;dy;dle;dli])::Vector{T}, KKTinv
end

"""
	function solveAugKKT(G,c,Ae,be,Ai,bi,x,y,le,li,rd,rpe,rpi,rli;KKTinv=[],doClear=true)

	builds and solves the augmented system. works so-so in practice
"""
function solveAugKKT{T}(G::AbstractArray{T},Ae::AbstractArray{T},Ai::AbstractArray{T},
	                  y::Vector{T},li::Vector{T},rd::Vector{T},rpe::Vector{T},rpi::Vector{T},rli::Vector{T};KKTinv=lufact(speye(1)),doClear=true)
	# build augmented matrix
    n  = size(G,1)
    me = size(Ae,1)
    mi = size(Ai,1)

    In  = speye(n)
    Imi = speye(mi)
    Znmi = spzeros(n,mi)
	Znme = spzeros(n,me)
	Znn = spzeros(n,n)
    Zmimi = spzeros(mi,mi)
	Zmeme = spzeros(me,me)
	Zmime = spzeros(mi,me)

    Y   = sdiag(y)
    Li  = sdiag(li)
	yLi = sdiag(y./li)

	rhs  = -[rd;rpe;rpi+(rli./li)]
	if size(KKTinv,1)!=size(rhs,1) || doClear
		# println("solveAugKKT: rebuild!")
    	KKT  =  [G -Ae' -Ai'; Ae  Zmeme Zmime'; Ai  Zmime yLi]
		KKTinv = lufact(KKT)
	end
	dt   = KKTinv\rhs::Vector{T}

	# solve for dy
	dx   = dt[1:n]
	dle  = dt[n+(1:me)]
	dli  = dt[n+me+1:end]
	dy   = (-rli-y.*dli)./li

	return vec([dx;dy;dle;dli])::Vector{T},KKTinv
end

"""
	function solveFullKKT(G,c,Ae,be,Ai,bi,x,y,le,li,rd,rpe,rpi,rli;KKTinv=[],doClear=true)

	builds and solves the full KKT system. Not recommended in practice
"""
function solveFullKKT{T}(G::AbstractArray{T},Ae::AbstractArray{T},Ai::AbstractArray{T},
	                  y::Vector{T},li::Vector{T},rd::Vector{T},rpe::Vector{T},rpi::Vector{T},rli::Vector{T};
	                  KKTinv=lufact(speye(1)),doClear=true)
	rhs = -[rd;rpe;rpi;rli]
	if size(KKTinv,1)!=size(rhs,1) || doClear
		# println("solveFull: rebuild!")
		KKT    = getJ(G,Ae,Ai,y,li)
		KKTinv = lufact(KKT)
	end

	return vec(KKTinv\rhs)::Vector{T}, KKTinv
end

"""
	function ipqp(G,c,Ae,be,Ai,bi,x,y,le,li;kwargs...)

	interior point solver for quadratic program based on Mehrotra's predictor-corrector approach

		min_x .5*x'*G*c + c'*x  subject to Ae*x = be and Ai*x >= bi

	Implementation is based on the presentation in Chapter 16 of

	Nocedal, J., & Wright, S. (2006). Numerical Optimization.
		Springer Science & Business Media. http://doi.org/10.1007/978-0-387-40065-5

	Optional keyword arguments:

		solveKKT          - function handle for linear solver. Options are currently
		                        solveNormalEqKKT (default)  (solves spsd linear system)
		                        solveAugKKT (symmetric indefinite solver)
		                        solveFullKKT
		maxIter::Int      - maximum number of iterations (default=10)
		ftol::Real        - stopping tolerance for norm of gradient of Lagrangian (default:1e-4)
		mutol::Real       - stopping tolerance on complementarity measure (default=1e-4)
		storeInterm::Bool - flag for storing intermediate iterates (default=false)
		out::Int          - controls verbosity
		                    (<1 --> no output, ==1 --> print final status, ==2 --> print status in each iteration)
		eta::Real         - reduction factor for steps (to remain interior) default=0.999

"""
function ipqp{T}(G::AbstractArray{T},c::Vector{T},Ae::AbstractArray{T},be::Vector{T},Ai::AbstractArray{T},bi::Vector{T},
	             x::Vector{T},y::Vector{T},le::Vector{T},li::Vector{T};maxIter::Int=10,ftol::Float64=1e-4,mutol::Float64=1e-4,
                 storeInterm::Bool=false,out::Int=1,eta::Float64=.999,solveKKT::Function=solveNormalEqKKT)

    n  = size(G,1)
    me = size(Ae,1)
    mi = size(Ai,1)

    his = zeros(maxIter+2,7)
    if storeInterm
        X  = zeros(n,maxIter)
        Y  = zeros(mi,maxIter)
        Le = zeros(me,maxIter)
        Li = zeros(mi,maxIter)
    end

	# 0) compute good starting guess (make sure y and li are sufficiently far in interior)
	mu  = dot(y,li)/mi
	rd,rpe,rpi,rli, = getF(G,c,Ae,be,Ai,bi,x,y,le,li)
	d,KKTinv  = solveKKT(G,Ae,Ai,y,li,rd,rpe,rpi,rli)
    (dx0,dy0,dle0,dli0) = splitVecIPQP(d,n,me,mi)
    y  = max.(1,abs.(y+dy0))
    li = max.(1,abs.(li+dli0))
	nrmrpe = (isempty(rpe)) ? 0.0 : norm(rpe)
   his[1,:] = [norm([rd;rpe;rpi;rli]) mu nrmrpe norm(rpi) norm(rd) norm(rli)  0]

    if out==2
        @printf "=== Interior Point for QP (n=%d,me=%d,mi=%d) ===\n" n me mi
        @printf "%4s\t%3s\t\t%2s\t\t%5s\t\t%5s\t\t%4s\t\t%4s\t\t%5s\n" "iter" "|F|" "mu" "|rpe|" "|rpi|" "|rd|" "|rl|" "alpha"
        @printf "%4d\t%1.2e\t%1.2e\t%1.2e\t%1.2e\t%1.2e\t%1.2e\n" -1 his[1,1] his[1,2] his[1,3] his[1,4] his[1,5]  his[1,6]
   end

    flag = -1; iter = 1;
    for iter=1:maxIter
		# get centrality and optimality conditions
		mu  = dot(y,li)/mi
		rd,rpe,rpi,rli, = getF(G,c,Ae,be,Ai,bi,x,y,le,li)

		# save and print history
        nrmrpe = (isempty(rpe))? 0.0 : norm(rpe)
		his[iter+1,1:6] = [norm([rd;rpe;rpi;rli]) mu nrmrpe norm(rpi) norm(rd) norm(rli)]
        if out==2
            @printf "%4d\t%1.2e\t%1.2e\t%1.2e\t%1.2e\t%1.2e\t%1.2e" iter-1 his[iter,1] his[iter,2] his[iter,3] his[iter,4] his[iter,5]  his[iter,6]
        end

		# check stopping criteria
		if his[iter,1]<ftol && (his[iter,2]<mutol)
			(out==2) && @printf "\n"
            flag = 0
			iter -=1
            break
        end

		# 1) compute an affine scaling step (sigma=0)
		d,KKTinv           = solveKKT(G,Ae,Ai,y,li,rd,rpe,rpi,rli,KKTinv=KKTinv,doClear=true)
        (dxaff,dyaff,dleaff,dliaff) = splitVecIPQP(d,n,me,mi)
        alphaaff = stepLength(y,li,dyaff,dliaff;eta=eta)

 		# 2) compute centrality parameter
        muaff = dot(y+alphaaff*dyaff, li+alphaaff*dliaff)/mi
		sigma = (muaff/mu)^3

		# 3) solve for new step
		rd,rpe,rpi,rli, = getF(G,c,Ae,be,Ai,bi,x,y,le,li,dy=dyaff,dli=dliaff,tau=sigma*mu)
		d,KKTinv        = solveKKT(G,Ae,Ai,y,li,rd,rpe,rpi,rli,KKTinv=KKTinv,doClear=false)
        (dx,dy,dle,dli) = splitVecIPQP(d,n,me,mi)
        alpha = stepLength(y,li,dy,dli;eta=0.999)

        # 4) update
        x   += alpha*dx
        y   += alpha*dy
        le  += alpha*dle
        li  += alpha*dli

		his[iter,7] = alpha
		(out==2) && @printf "\t%1.2e\n" alpha

        if storeInterm
            X[:,iter]  = x
            Y[:,iter]  = y
            Le[:,iter] = le
            Li[:,iter] = li
        end
    end

    if out>=0
        if flag==-1
            @printf "ipqp iterated maxIter(=%d) times but reached only |F(x,s,u,v)|=%1.2e (ftol=%1.2e) and mu=%1.2e (mutol=%1.2e)\n" maxIter his[iter+1,1] ftol his[iter+1,2] mutol
        elseif flag==0 && out>=1
            @printf "ipqp achieved ftol=%1.2e and mutol=%1.2e at iteration %d. Returned result satisfies |F(x,y,s)|=%1.2e and mu=%1.2e.\n" ftol  mutol iter his[iter,1] his[iter,2]
        end
    end

    (x,y,le,li) = (storeInterm) ? (X[:,1:iter],Y[:,1:iter],Le[:,1:iter],Li[:,1:iter]) : (x,y,le,li)
    return x,y,le,li,his[1:iter+1,:],KKTinv
end


text:
include("twihes.jl")

############################
#
# Parameters
#
############################

#Prime and field
const la = big(2)
const lb = big(3)
const ea = 305
const eb = 192
const f = big(1)
const p = la^ea*lb^eb*f-1
const F = Fp2(p)

#cuberoot of One
const w = F(big(1319470205536631335981810349644143385091780115593611158375203778232819918005279075081612765167581266740524628378608982325666905211062864268703698577903039608673459647224627806555078655), big(1286451788956606916776387017696978305185994757974417726547130267369253169617279161336419642919745234101683127375727357361442589772170347650505223652853823514074141695255872287918030044))
if !isone(w^3)
    error("w is not a cuberoot...")
end
#Starting curve: a = 1 for speed, solve for d to get j(E) = some supersingular j-invariant. See sage-script
#Starting j-invariant = 831837561004103883868671088207491838344992437494159978595307997269883652199194625608995144292694837047770068423678474079025058740865040338182857297729622238883560881131572628427170113 + 771813078065723508849643803754084945896394345666593331261515932318350081581366794480345255301600363646900179413481795555101011900773965712484190126442309614812527433851551995117691457*i
#See sage script for construction of d0 from j-Invariant and a0
const a0 = one(F)
const d0 = F(big(526466053434329606348792315213226961488878919273384911349255872484662011698700897679785221256842528945834974679658632918959678155926738589102758106233735036087191018029528283224394826), big(322622151547826919747264788491387568702182872571850474819101171989083545313493484238990598522402415188382526845235802789803276920938914974974558706304117646595374335297821576916453987))

#Generates torsion basises from random points
function genTorsionBasis(alicePX, alicePY, aliceQX, aliceQY, bobPX, bobPY, bobQX, bobQY)
    Z = one(F)
    aliceP = EC(alicePX, alicePY, Z)
    aliceQ = EC(aliceQX, aliceQY, Z)

    if !isValid(aliceP, a0, d0)
        error("Alice' P is not a valid point!")
    end

    if !isValid(aliceQ, a0, d0)
        error("Alice' Q is not a valid point!")
    end

    Pa = normalized(mul(f*lb^eb, aliceP, a0))
    Qa = normalized(mul(f*lb^eb, aliceQ, a0))

    Patest = doublePow(Pa, a0, ea-1)
    Qatest = doublePow(Qa, a0, ea-1)

    if !isIdentity(double(Patest, a0))
        error("Pa : Something does not look right with the order of the curve...")
    end

    if !isIdentity(double(Qatest, a0))
        error("Qa : Something does not look right with the order of the curve...")
    end

    if isIdentity(Patest)
        error("Alice' point Pa does not have the right order!")
    end
    if isIdentity(Qatest)
        error("Alice' point Qa does not have the right order!")
    end

    if Patest == Qatest
        error("Alice' points are not independent!")
    end

    bobP = EC(bobPX, bobPY, Z)
    bobQ = EC(bobQX, bobQY, Z)

    if !isValid(bobP, a0, d0)
        error("Bob's P is not a valid point!")
    end

    if !isValid(bobQ, a0, d0)
        error("Bob's Q is not a valid point!")
    end

    Pb = normalized(mul(f*la^ea, bobP, a0))
    Qb = normalized(mul(f*la^ea, bobQ, a0))

    Pbtest = triplePow(Pb, a0, d0, eb-1)
    Qbtest = triplePow(Qb, a0, d0, eb-1)

    if !isIdentity(triple(Pbtest, a0, d0))
        error("Pb : Something does not look right with the order of the curve...")
    end

    if !isIdentity(triple(Qbtest, a0, d0))
        error("Qb : Something does not look right with the order of the curve...")
    end

    if isIdentity(Pbtest)
        error("Bob's point Pb does not have the right order!")
    end

    if isIdentity(Qbtest)
        error("Bob's point Qb does not have the right order!")
    end

    if Pbtest == Qbtest
        error("Bob's points are not independent!")
    end

    return Pa, Qa, Pb, Qb
end

#See sage script for generating random points
PX1 = F(big(2492006491871983676452451839020564198147973496177388329872836720660635968389647718878527574597574695774796088332985623054551169799526739171989607694478941255607446159858803583994964802), big(1106461312022518869103531097009806634520011809887491682622823938979378314609312783672420779078369845610390072881748838641128935933666342275393037280723333398214385754330629136835945615))
PY1 = F(big(568934704890564342552989237837373499705027664549344609870746088351623756422430685800281891502113550175496052961717514133893321580563982611360102558084198226260195186543398872742166956), big(1719319066282017378159238737975467015356375587815788392734006098067712564812437626919347516402785703024563499386845184938500959018676538501237827797206911141692035339948528241047463062))
QX1 = F(big(2517628532562506074245481407647129152141059415595858026166161117305445674436190840142297568565926716458956505364881489417252719704505645622987495305020004927497493134683221021490239959), big(1433263917761027493460060659243243110082141390402415821474572364169645458452572424878937372751197553782045370168747546443025397983411476342832168979291996203472762061817060320239732514))
QY1 = F(big(1268486271669556135524139358247425138875339494127874128975087540722306664933940157281323943300180893898834305285190175720304656654248228384119897648002761039198456950277535984384435960), big(870717938688314703034698794917729951973369728070119897658763707633485814633592046219916287841159157678092192037611181964671011367024700114197655484075793009679462073178107401411271146))
PX2 = F(big(1087251779026402284780690514599816326569275126413001776550994011768620358684555955678133912577235126758668170004393171497761389593621955822052207714093613781329468475697631867387474483), big(2469211885106793446942611087302463642283587716018469215650568305955990628425421892841320028890624561652810865126897501340849958642894467779823256283556264901094677686850284342042454424))
PY2 = F(big(2042684845558020702068009872943604896265659382991530782022341450605024204133331153581602581183717652980627685160557204787602548705930777904538999384048875962096503211454749544788971900), big(1952407962161458071632111111545938309344853088503297399264610051714768899866468987194505212765542390314727724800985868656215669828095236952855844761115022653991396268869521633770184697))
QX2 = F(big(1948557883143239148216542337109344092480664013194035763418796831802277689220591717171651875237624876590110990487339142795957213635061021311223506194896335710514651491187248696544839890), big(1561987540323249359134988383755115761443665949729569112934702895896188229271478800280961840021334352088226396987001996010282799910199376720872821614759744240212908203945883767928868921))
QY2 = F(big(2434803430031357671431182080606412625586548768754456399219302240585675185680813127143693009941726516214200502993859529022071607523524471311775260681161873786254554974167028608047308730), big(1246986093659258949315414845812705415673341880493879348339825574550001003997591027583897079561051808540633480867564839740453013245250501290325921658949100205817371785062094851400235056))
const Pa, Qa, Pb, Qb = genTorsionBasis(PX1, PY1, QX1, QY1, PX2, PY2, QX2, QY2)

# Compute the optimal strategy of size n, given weights p and q
function computeStrat(n, p, q)
    S = Dict{Int64, Array{Int64, 1}}()
    C = Dict{Int64, Int64}()
    S[1] = []
    C[1] = 0
    for i in 2:n
        b = argmin([C[i-b] + C[b] + b*p + (i-b)*q for b in 1:i-1])
        S[i] = append!(append!([b], S[i-b]), S[b])
        C[i] = C[i-b] + C[b] + b*p + (i-b)*q
    end
    return S[n]
end

#The weights are estimated using @benchmark
const strat_A = Vector(computeStrat(ea, 23358, 32132))

const strat_B = Vector(computeStrat(eb, 45620, 34278))


text:
# Script which loads the DtN nep.
import Base.size
using SpecialFunctions,Printf, SparseArrays;
import NonlinearEigenproblems.compute_Mder;
import NonlinearEigenproblems.compute_Mlincomb;

# A DtN-nep
struct BesselNEP <: NEP
    Q::Matrix{ComplexF64};
    P::Vector{SparseMatrixCSC{ComplexF64}} # Full matrix coefficents
    ind2::Vector;
    n::Int;
end

size(nep::BesselNEP) = (nep.n,nep.n)
size(nep::BesselNEP,dim::Int) = nep.n

function besselh_quotient(nu,S)
# Compute: besselh'(nu,s)/besselh(nu,s)
    Fder=0.5*(besselh(nu-1,S)-besselh(nu+1,S));
    return besselh(nu,S)\Fder;
end

function besselh_quotient_der(nu,S)
    # Compute: derivative of besselh'(nu,s)/besselh(nu,s)
    Fderder=0.25*(besselh(nu-2,S)-2*besselh(nu,S)+besselh(nu+2,S));
    Fder=0.5*(besselh(nu-1,S)-besselh(nu+1,S))
    F=besselh(nu,S);
    return (F^2)\(Fderder*F-Fder*Fder);
end

function compute_Mder(nep::BesselNEP, λ::Number, der::Int=0)
    if (der>0)
        error("Not implemented");
    end
    A=zero(nep.P[1]);
    a=1.0;
    for i=1:length(nep.ind2)
        # add:  - s B'_m(s)/B_m(s)
        m=nep.ind2[i];
        f= S ->  -S*besselh_quotient(m,a*S);
        A += f(λ)*nep.P[i]
    end
    return A;
end

function compute_Mlincomb(nep::BesselNEP, λ::Number, V::AbstractVecOrMat)

    if (size(V,2)>2)
        error("Higher derivatives not implemented")
    end
    TT=promote_type(eltype(V),typeof(λ),eltype(nep.P[1]));
    n=size(nep,1);

    v = zeros(TT, size(nep,1));
    a=1.0;
    if (norm(V[:,1])>0)
        # Zero'th derivative:
        W0::Vector{TT}=nep.Q[:,1:size(nep.ind2,1)]'*V[:,1];
        z=Vector{TT}(undef,size(nep.ind2,1));

        @inbounds for i=1:size(nep.ind2,1)
            # add:  - s B'_m(s)/B_m(s)
            m=nep.ind2[i];
            f = S ->  -S*besselh_quotient(m,a*S);
            z[i] = W0[i]*f(λ);
        end
        v[:] .= nep.Q[:,1:size(nep.ind2,1)]*z;
    end
    if (size(V,2)>1)
        # First derivative:
        W1::Vector{TT}=nep.Q[:,1:size(nep.ind2,1)]'*V[:,2];
        @inbounds for i=1:size(nep.ind2,1)
          # add derivative of: - s B'_m(s)/B_m(s)
          m=nep.ind2[i];
          f= S ->  -besselh_quotient(m,a*S)-a*S*besselh_quotient_der(m,a*S)
          v[:] .+= nep.Q[:,i]*(W1[i]*f(λ))
       end
    end
    return v;
end

## Load the matrices

include("petsc_naive_bin_read.jl")

# l=40; # Number of terms

"""
    function load_dtn_dimer(data_dir::String,l::Int)

 Loads the DtN example in the NEP described in J. Araujo-Cabarcas, C. Engström and E. Jarlebring, Efficient resonance computations for Helmholtz problems based on a Dirichlet-to-Neumann map, J. Comput. Appl. Math., 330:177-192, 2018  (http://arxiv.org/pdf/1606.09547)
"""
function load_dtn_dimer(data_dir::String,l::Int)

    A = naive_petsc_read(joinpath(data_dir,"K.bin"))
    M = naive_petsc_read(joinpath(data_dir,"M.bin"))
    n = size(A,1);

    # config_fem.m contains dofs and dofs_boundaries.
    # This should be the same as dof-dofs_boundaries
    q1 = naive_petsc_read(joinpath(data_dir,"q1.bin"))
    start_dtn_dofs=findfirst(abs.(q1).>0)-1;

    dtn_a = 1;

    # Determine the number of saved terms.
    searchdir(path,key) = filter(x->startswith(x,key), readdir(path))
    files=searchdir(data_dir,"q") # All the Q-vectors
    # Mid of the number of saved terms.
    mid=round(Int,(size(files,1)-1)/2+1);

    l=min((mid-1),l); # make sure we don't exceed the dimension

    ind=Vector((mid-l):(mid+l));
    ind2=ind .- mid;

    Q=Matrix{ComplexF64}(undef,n,size(ind,1)); # Factorized matrices
    P=Vector{SparseMatrixCSC}(undef,size(ind,1)) # A vector of rank one matrices
    for i=1:size(ind,1)
        # This corresponds to the term q*q'*bessel_quotient(ind2[i],s)
        # Load the q-vector.
        fname=@sprintf("q%d.bin",ind[i]);
        q = naive_petsc_read(joinpath(data_dir,fname))/(sqrt(2*pi)*sqrt(dtn_a))
        Q[:,i]= q;
        # Multiply it out (a bit complicated since since we want to save memory)
        # NNZ part of the q-vector
        qnz = q[(start_dtn_dofs+1):end]
        Qnz=sparse(qnz*qnz')
        (I,J,V)=findnz(Qnz)
        P[i]=sparse(I.+start_dtn_dofs,J.+start_dtn_dofs,V);
    end

    nep1=SPMF_NEP([A, M], [S-> one(S), S->-S^2]);
    nep2=BesselNEP(Q,P,ind2,n);
    nep=SumNEP(nep1,nep2);
    return nep

end


text:
# name = "GAMS"
# GAMS = "1ca51c6a-1b4d-4546-9ae1-53e0a243ab12"

using GAMS

GAMS_OPT = GAMS.Optimizer()
MOI.set(GAMS_OPT, MOI.Silent(), true)
MOI.set(GAMS_OPT, MOI.TimeLimitSec(), 5)
# GAMS_OPT = MOI.Bridges.Constraint.SOCtoNonConvexQuad{Float64}(GAMS_OPT)

push!(solvers_complements, (opt = GAMS_OPT, mode = BilevelJuMP.ComplementMode()))


text:
export get_index, ensure_index

typealias BamIndex Dict{String, IntRangeDict{i32, i32}}

function get_index(bam::AbstractBam)::BamIndex
    index = load_index(bam)
    index.isnull ? make_index(bam) : index.value
end

function make_index(bam::AbstractBam)::BamIndex
    index = BamIndex()
    chr = -2
    local dict::IntRangeDict{i32, i32}
    for (idx, read) in enumerate(bam) @when read.refID >= 0
        if read.refID != chr
            chr = read.refID
            index[bam.refs[chr+1] |> car] = dict = IntRangeDict{i32, i32}()
        end

        start = read.pos |> i32
        stop = read.pos + calc_distance(read) - 1 |> i32

        push!(dict[start:stop], i32(idx))
    end
    index
end

function save_index(bam::AbstractBam, index::BamIndex)
    if isempty(bam.file)
        error("cannot save index of stream bam")
    end

    bamtime = mtime(bam.file)
    f = open(bam.file * ".fbi", "w")

    f << b"FBIv1"
    write(f, bamtime)

    f = ZlibDeflateOutputStream(f)

    for (key, value) in index
        f << key << '\0'
        save(f, value)
    end

    close(f)
end

function load_index(bam::AbstractBam)::Nullable{BamIndex}
    if isempty(bam.file) || !isfile(bam.file * ".fbi")
        nothing
    end

    bamtime = mtime(bam.file)
    f = open(bam.file * ".fbi")

    try
        index = BamIndex()
        @assert f >> 5 == b"FBIv1"
        @assert f >> f64 == bamtime

        f = ZlibInflateInputStream(f)

        while !eof(f)
            name = readuntil(f, '\0') |> del_end!
            dict = IntRangeDict{i32, i32}(f)
            index[name] = dict
        end
        index
    catch
        nothing
    finally
        close(f)
    end
end

function ensure_index(bam::AbstractBam)
    if isempty(bam.file)
        error("cannot save index of stream bam")
    end

    bamtime = mtime(bam.file)

    if isfile(bam.file * ".fbi")
        try
            open(bam.file * ".fbi") do f
                f >> 5 == b"FBIv1" && f >> f64 == bamtime
            end && return
        end
    end

    save_index(bam, make_index(bam))
end


text:
# common facilities

function augment_rows{T}(X::StridedMatrix{T}, v::T)
    m, n = size(X)
    Xa = Array(T, m, n+1)
    copy!(view(Xa, :, 1:n), X)
    fill!(view(Xa, :, n+1), v)
    Xa
end

function augment_cols{T}(X::StridedMatrix{T}, v::T)
    m, n = size(X)
    Xa = Array(T, m+1, n)
    copy!(view(Xa, 1:m, :), X)
    fill!(rowvec_view(Xa, m+1), v)
    Xa
end

## verbosity indicator

const VERBOSE_NONE  = 0
const VERBOSE_FINAL = 1
const VERBOSE_ITER  = 2

verbosity_level(s::Symbol) = s == :none  ? VERBOSE_NONE :
                             s == :final ? VERBOSE_FINAL :
                             s == :iter  ? VERBOSE_ITER :
                             error("Invalid verbosity symbol.")


text:
maxiter = 1000

using Stopping

stp = NLPStopping(nlp, NLPAtX(nlp.meta.x0)  )

Lp = Inf

my_unconstrained_check(nlp, st; kwargs...) = unconstrained_check(nlp, st, pnorm = Lp; kwargs...)
stp.meta.optimality_check = my_unconstrained_check

#stp.meta.optimality_check = unconstrained_check
stp.meta.max_iter = maxiter
stp.meta.rtol = 0

using OneDmin

function test_Stp(algo::Function, nlp; kwargs...) where T

    stats = @timed bidon=0
    logger = Logging.ConsoleLogger(stderr,Logging.Warn)
    Logging.with_logger(logger) do
        stats = @timed stp = algo(nlp; kwargs...)
    end

    return stats, stp
end

function test_noStp(algo::Function, nlp; kwargs...) where T

    stats = @timed bidon=0
    iter=0
    f=0.0
    normg=0.0
    logger = Logging.ConsoleLogger(stderr,Logging.Warn)
    Logging.with_logger(logger) do
        stats = @timed iter, f, normg, B = algo(nlp; kwargs...)
    end

    return stats, iter, f, normg
end

using Test
@info log_header([:name, :time, :iter, :f, :dual], [String, Float64, Int, Float64, Float64],
                 hdr_override=Dict(:name => "Solver", :f=>"f(x)", :dual=>"‖∇f‖"))

println("\n CG   ")

stats, stp = test_Stp(CG_generic, nlp, stp=stp, scaling = false, strongWolfe = true, LS_algo=bracket_N, CG_formula=formula_HZ)

@info log_row(Any["CG_HZ", stats.time,  stp.meta.nb_of_stop, stp.current_state.fx, stp.current_state.current_score])
@test stp.current_state.current_score < 1e-6

reset!(nlp)
reinit!(stp)

stats, stp = test_Stp(CG_HZ, nlp, stp=stp, scaling = false, strongWolfe = true, LS_algo=bracket_N)

@info log_row(Any["CG_HZ wrapper", stats.time,  stp.meta.nb_of_stop, stp.current_state.fx, stp.current_state.current_score])
@test stp.current_state.current_score < 1e-6

reset!(nlp)
reinit!(stp)

stats, stp = test_Stp(CG_HS, nlp, stp=stp, scaling = false, strongWolfe = true, LS_algo=bracket_N)

@info log_row(Any["CG_HS", stats.time,  stp.meta.nb_of_stop, stp.current_state.fx, stp.current_state.current_score])
@test stp.current_state.current_score < 1e-6

reset!(nlp)
reinit!(stp)

stats, stp = test_Stp(CG_PR, nlp, stp=stp, scaling = false, strongWolfe = true, LS_algo=bracket_N)

@info log_row(Any["CG_PR", stats.time,  stp.meta.nb_of_stop, stp.current_state.fx, stp.current_state.current_score])
@test stp.current_state.current_score < 1e-6

reset!(nlp)
reinit!(stp)

stats, stp = test_Stp(CG_generic, nlp, stp=stp, scaling = true, strongWolfe = true, LS_algo=bracket_N, CG_formula=formula_FR)

@info log_row(Any["CG_FR", stats.time,  stp.meta.nb_of_stop, stp.current_state.fx, stp.current_state.current_score])
@test stp.current_state.current_score < 1e-6

println("\n bfgs   ")

reset!(nlp)
reinit!(stp)
Lp = 2.0
my_unconstrained_check(nlp, st; kwargs...) = unconstrained_check(nlp, st, pnorm = Lp; kwargs...)
stp.meta.optimality_check = my_unconstrained_check

stats, stp = test_Stp(bfgs_StopLS, nlp, stp=stp, LS_algo=bracket_B)

@info log_row(Any["bfgsSLS-L2", stats.time,  stp.meta.nb_of_stop, stp.current_state.fx, stp.current_state.current_score])
@test norm(stp.current_state.gx, Lp) < 1e-6

reset!(nlp)
reinit!(stp)
Lp = Inf
my_unconstrained_check(nlp, st; kwargs...) = unconstrained_check(nlp, st, pnorm = Lp; kwargs...)
stp.meta.optimality_check = my_unconstrained_check

stats, stp = test_Stp(bfgs_StopLS, nlp, stp=stp, LS_algo=bracket_B)

@info log_row(Any["bfgsSLS-L∞", stats.time,  stp.meta.nb_of_stop, stp.current_state.fx, stp.current_state.current_score])
@test norm(stp.current_state.gx, Lp) < 1e-6

reset!(nlp)
reinit!(stp)

stats, stp = test_Stp(bfgs_Stop, nlp, stp=stp)

@info log_row(Any["bfgsS", stats.time,  stp.meta.nb_of_stop, stp.current_state.fx, stp.current_state.current_score])
@test stp.current_state.current_score < 1e-6

reset!(nlp)

stats, iter, f, g = test_noStp(bfgs, nlp, scaling = true, maxiter = maxiter, Lp = Inf)

@info log_row(Any["bfgs", stats.time,  iter, f, g])
@test g < 1e-6

println("\n L-bfgs   ")

mem = 70
reset!(nlp)
reinit!(stp)

stats, stp = test_Stp(L_bfgs_StopLS, nlp, stp=stp, LS_algo=bracket_B, mem = mem)

@info log_row(Any["L-bfgsSLS", stats.time,  stp.meta.nb_of_stop, stp.current_state.fx, stp.current_state.current_score])
@test stp.current_state.current_score < 1e-6

reset!(nlp)
reinit!(stp)

stats, stp = test_Stp(L_bfgs_Stop, nlp, stp=stp, mem = mem)

@info log_row(Any["L-bfgsS", stats.time,  stp.meta.nb_of_stop, stp.current_state.fx, stp.current_state.current_score])
@test stp.current_state.current_score < 1e-6

reset!(nlp)

stats, iter, f, g = test_noStp(L_bfgs, nlp, scaling = true, maxiter = maxiter, Lp = Inf, mem = mem)

@info log_row(Any["L-bfgs", stats.time,  iter, f, g])
@test g < 1e-6

println("\n M-bfgs   ")

reset!(nlp)
reinit!(stp)

stats, stp = test_Stp(M_bfgs_StopLS, nlp, stp=stp, LS_algo=bracket_B)

@info log_row(Any["M-bfgsSLS", stats.time,  stp.meta.nb_of_stop, stp.current_state.fx, stp.current_state.current_score])
@test stp.current_state.current_score < 1e-6

reset!(nlp)
reinit!(stp)

stats, stp = test_Stp(M_bfgs_Stop, nlp, stp=stp)

@info log_row(Any["M-bfgsS", stats.time,  stp.meta.nb_of_stop, stp.current_state.fx, stp.current_state.current_score])
@test stp.current_state.current_score < 1e-6

reset!(nlp)

stats, iter, f, g = test_noStp(M_bfgs, nlp, scaling = true, maxiter = maxiter, Lp = Inf)

@info log_row(Any["M-bfgs", stats.time,  iter, f, g])
@test g < 1e-6

println("\n scaling = false")

reset!(nlp)
reinit!(stp)

stats, stp = test_Stp(M_bfgs_StopLS, nlp, stp=stp, LS_algo=bracket_B, scaling = false)

@info log_row(Any["M-bfgsSLS", stats.time,  stp.meta.nb_of_stop, stp.current_state.fx, stp.current_state.current_score])
@test stp.current_state.current_score < 1e-6

mem = 70
reset!(nlp)
reinit!(stp)

stats, stp = test_Stp(L_bfgs_StopLS, nlp, stp=stp, LS_algo=bracket_B, mem = mem, scaling = false)

@info log_row(Any["L-bfgsSLS", stats.time,  stp.meta.nb_of_stop, stp.current_state.fx, stp.current_state.current_score])
@test stp.current_state.current_score < 1e-6

println("\n Newton   ")

reset!(nlp)
reinit!(stp)

stats, stp = test_Stp(Newton_StopLS, nlp, stp=stp)

@info log_row(Any["NwtSLS", stats.time,  stp.meta.nb_of_stop, stp.current_state.fx, stp.current_state.current_score])
@test stp.current_state.current_score < 1e-6

reset!(nlp)
reinit!(stp)

stats, stp = test_Stp(Newton_Stop, nlp, stp=stp)

@info log_row(Any["NwtS", stats.time,  stp.meta.nb_of_stop, stp.current_state.fx, stp.current_state.current_score])
@test stp.current_state.current_score < 1e-6

reset!(nlp)

stats, iter, f, g = test_noStp(Newton_Spectral, nlp,  maxiter = maxiter, Lp = Inf, ϵ = 1e-6)

@info log_row(Any["Nwt", stats.time,  iter, f, g])
@test g < 1e-6

reset!(nlp)
reinit!(stp)

stats, stp = test_Stp(NewtonSpectralAbs, nlp, stp=stp)

@info log_row(Any["Nwt_spec_abs", stats.time,  stp.meta.nb_of_stop, stp.current_state.fx, stp.current_state.current_score])
@test stp.current_state.current_score < 1e-6

reset!(nlp)
reinit!(stp)

stats, stp = test_Stp(NewtonSpectralAbs, nlp, stp=stp, scale_abs = false)

@info log_row(Any["Nwt_spec", stats.time,  stp.meta.nb_of_stop, stp.current_state.fx, stp.current_state.current_score])
@test stp.current_state.current_score < 1e-6

reset!(nlp)
reinit!(stp)

stats, stp = test_Stp(NewtonLDLtAbsLS, nlp, stp=stp)

@info log_row(Any["Nwt_LDLt_abs", stats.time,  stp.meta.nb_of_stop, stp.current_state.fx, stp.current_state.current_score])
@test stp.current_state.current_score < 1e-6

reset!(nlp)
reinit!(stp)

stats, stp = test_Stp(NewtonLDLtAbsLS, nlp, stp=stp, scale_abs = false)

@info log_row(Any["Nwt_LDLt", stats.time,  stp.meta.nb_of_stop, stp.current_state.fx, stp.current_state.current_score])
@test stp.current_state.current_score < 1e-6

reset!(nlp)
reinit!(stp)

stats, stp = test_Stp(NewtonCG, nlp, stp=stp)

@info log_row(Any["NwtCG", stats.time,  stp.meta.nb_of_stop, stp.current_state.fx, stp.current_state.current_score])
@test stp.current_state.current_score < 1e-6

max_stars_repo_name: vepiteski/LSDescent

max_issues_repo_name: vepiteski/LSDescent

max_forks_repo_name: vepiteski/LSDescent


text:
using ChebyshevQuantum
using LinearAlgebra

function go(N)

    D = ChebyshevQuantum.Interp.getD(N);
    pts = ChebyshevQuantum.Interp.getCpts(N);
    v = eigvals( (D*D )[2:N,2:N]);
    v / (pi^2/4)
end


text:
module DefaultGates

using QXZoo.GateOps

export x, y, z, h, c_x, c_y, c_z, u
export r_x, r_y, r_z, r_phase, c_r_x, c_r_y, c_r_z, c_u, c_r_phase

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% #
#                            Default gate symbols
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% #

module GateSymbols
    using QXZoo.GateOps

    """Projection operators"""
    p00 = GateOps.GateSymbol(:p00) # |0><0|
    p10 = GateOps.GateSymbol(:p10) # |1><0|
    p01 = GateOps.GateSymbol(:p01) # |0><1|
    p11 = GateOps.GateSymbol(:p11) # |1><1|

    """Single qubit"""
    x = GateOps.GateSymbol(:x)
    y = GateOps.GateSymbol(:y)
    z = GateOps.GateSymbol(:z)
    h = GateOps.GateSymbol(:h)
    s = GateOps.GateSymbol(:z, 1, false)
    t = GateOps.GateSymbol(:z, 2, false)
    I = GateOps.GateSymbol(:I)

    r_x(θ) = GateOps.GateSymbolP(:r_x, 0, false, θ)
    r_y(θ) = GateOps.GateSymbolP(:r_y, 0, false, θ)
    r_z(θ) = GateOps.GateSymbolP(:r_z, 0, false, θ)
    r_phase(θ) = GateOps.GateSymbolP(:r_phase, 0, false, θ)

    """Two qubit"""
    c_x = GateOps.GateSymbol(:c_x)
    c_y = GateOps.GateSymbol(:c_y)
    c_z = GateOps.GateSymbol(:c_z)
    c_H = GateOps.GateSymbol(:c_h)

    c_r_x(θ) = GateOps.GateSymbolP(:c_r_x, 0, false, θ)
    c_r_y(θ) = GateOps.GateSymbolP(:c_r_y, 0, false, θ)
    c_r_z(θ) = GateOps.GateSymbolP(:c_r_z, 0, false, θ)
    c_r_phase(θ) = GateOps.GateSymbolP(:c_r_phase, 0, false, θ)
end

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% #
#                            1Q Gate calls ops
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% #

"""
    x(q_target::Int)

Generate a single qubit Pauli-x GateCall (GateCall1) applied to the target qubit

# Examples
```julia-repl
julia> QXZoo.DefaultGates.x(0)
QXZoo.GateOps.GateCall1(QXZoo.GateOps.GateSymbol(:x, 0, false), 0)
```
"""
function x(q_target::Integer)
    return GateOps.GateCall1(GateSymbols.x, q_target)
end

"""
    y(q_target::Int)

Generate a single qubit Pauli-y GateCall (GateCall1) applied to the
target qubit

# Examples
```julia-repl
julia> QXZoo.DefaultGates.y(0)
QXZoo.GateOps.GateCall1(QXZoo.GateOps.GateSymbol(:y, 0, false), 0)
```
"""
function y(q_target::Int)
    return GateOps.GateCall1(GateSymbols.y, q_target)
end

"""
    z(q_target::Int)

Generate a single qubit Pauli-z GateCall (GateCall1) applied to the
target qubit

# Examples
```julia-repl
julia> QXZoo.DefaultGates.z(0)
QXZoo.GateOps.GateCall1(QXZoo.GateOps.GateSymbol(:z, 0, false), 0)
```
"""
function z(q_target::Int)
    return GateOps.GateCall1(GateSymbols.z, q_target)
end

"""
    h(q_target::Int)

Generate a single qubit Hadamard GateCall (GateCall1) applied to the
target qubit

# Examples
```julia-repl
julia> QXZoo.DefaultGates.h(0)
QXZoo.GateOps.GateCall1(QXZoo.GateOps.GateSymbol(:h, 0, false), 0)
```
"""
function h(q_target::Int)
    return GateOps.GateCall1(GateSymbols.h, q_target)
end

"""
    s(q_target::Int)

Generate a single qubit S GateCall (GateCall1) applied to the
target qubit (s == sqrt(z))

# Examples
```julia-repl
julia> QXZoo.DefaultGates.s(0)
QXZoo.GateOps.GateCall1(QXZoo.GateOps.GateSymbol(:z, 1, false), 0)
```
"""
function s(q_target::Int)
    return GateOps.GateCall1(GateSymbols.s, q_target)
end

"""
    s(q_target::Int)

Generate a single qubit T GateCall (GateCall1) applied to the
target qubit (t == sqrt(s))

# Examples
```julia-repl
julia> QXZoo.DefaultGates.t(0)
QXZoo.GateOps.GateCall1(QXZoo.GateOps.GateSymbol(:z, 2, false), 0)
```
"""
function t(q_target::Int)
    return GateOps.GateCall1(GateSymbols.t, q_target)
end

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% #

"""
    u(label::GateSymbol, q_target::Int)

Generate a single qubit arbitrary unitary GateCall (GateCall1) applied to the
target qubit

# Examples
```julia-repl
julia> QXZoo.DefaultGates.u(QXZoo.GateOps.GateSymbol(:mygate), 1)
QXZoo.GateOps.GateCall1(QXZoo.GateOps.GateSymbol(:mygate, 0, false), 1)
```
"""
function u(label::GateSymbol, q_target::Int)
    return GateOps.GateCall1( label, q_target)
end
"""
    u(label::String, q_target::Int)

Generate a single qubit arbitrary unitary GateCall (GateCall1) applied to the
target qubit

# Examples
```julia-repl
julia> QXZoo.DefaultGates.u("mygate", 1)
QXZoo.GateOps.GateCall1(QXZoo.GateOps.GateSymbol(:mygate, 0, false), 0)
```
"""
function u(label::String, q_target::Int)
    return GateOps.GateCall1( GateOps.GateSymbol(Symbol(label)), q_target)
end

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% #

"""
    r_x(q_target::Int, theta::Number)

Generate a single qubit R_x(θ) GateCall (GateCall1) applied to the target qubit

# Examples
```julia-repl
julia> QXZoo.DefaultGates.r_x(3,pi/2)
QXZoo.GateOps.GateCall1(QXZoo.GateOps.GateSymbolP(:r_x, 0, false, 1.5707963267948966), 3)
```
"""
function r_x(q_target::Int, theta::Number)
    return GateOps.GateCall1(GateSymbols.r_x(theta), q_target )
end

"""
    r_y(q_target::Int, theta::Number)

Generate a single qubit R_y(θ) GateCall (GateCall1) applied to the target qubit.

# Examples
```julia-repl
julia> QXZoo.DefaultGates.r_y(3,pi/2)
QXZoo.GateOps.GateCall1(QXZoo.GateOps.GateSymbolP(:r_y, 0, false, 1.5707963267948966), 3)
```
"""
function r_y(q_target::Int, theta::Number)
    return GateOps.GateCall1(GateSymbols.r_y(theta), q_target )
end

"""
    r_z(q_target::Int, theta::Number)

Generate a single qubit R_z(θ) GateCall (GateCall1) applied to the target qubit

# Examples
```julia-repl
julia> QXZoo.DefaultGates.r_z(3,pi/2)
QXZoo.GateOps.GateCall1(QXZoo.GateOps.GateSymbolP(:r_z, 0, false, 1.5707963267948966), 3)
```
"""
function r_z(q_target::Int, theta::Number)
    return GateOps.GateCall1(GateSymbols.r_z(theta), q_target )
end

"""
    r_phase(q_target::Int, theta::Number)

Generate a single qubit phase shift GateCall ( diag(1, exp(1im*theta)) , GateCall1) applied to the target qubit

# Examples
```julia-repl
julia> QXZoo.DefaultGates.r_phase(3,pi/2)
QXZoo.GateOps.GateCall1(QXZoo.GateOps.GateSymbolP(:r_phase, 0, false, 1.5707963267948966), 3)
```
"""
function r_phase(q_target::Int, theta::Number)
    return GateOps.GateCall1(GateSymbols.r_phase(theta), q_target )
end

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% #
#                            2Q Gate calls ops
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% #

"""
    swap(q_target::Int, q_ctrl::Int)

Apply swap gate between given qubits

# Examples
```julia-repl
julia> QXZoo.DefaultGates.swap(0,1)
QXZoo.GateOps.GateCall2(QXZoo.GateOps.GateSymbol(:swap, 0, false), 0, 1, nothing)
```
"""
function swap(q_target::Int, q_ctrl::Int)
    return GateOps.GateCall2( GateOps.GateSymbol(:swap), q_target, q_ctrl )
end

"""
    c_x(q_target::Int, q_ctrl::Int)

Generate a controlled Pauli-x (two qubit) GateCall (GateCall2), controlled on the index `q_ctrl` applied to the target `q_target`.

# Examples
```julia-repl
julia> QXZoo.DefaultGates.c_x(0, 1)
QXZoo.GateOps.GateCall2(QXZoo.GateOps.GateSymbol(:c_x, 0, false), 0, 1, QXZoo.GateOps.GateSymbol(:x, 0, false))
```
"""
function c_x(q_target::Int, q_ctrl::Int)
    return GateCall2( GateSymbols.c_x, q_target, q_ctrl, GateSymbols.x )
end

"""
    c_y(q_target::Int, q_ctrl::Int)

Generate a controlled Pauli-y (two qubit) GateCall (GateCall2), controlled on the index `q_ctrl` applied to the target `q_target`.

# Examples
```julia-repl
julia> QXZoo.DefaultGates.c_y(0, 1)
QXZoo.GateOps.GateCall2(QXZoo.GateOps.GateSymbol(:c_y, 0, false), 0, 1, QXZoo.GateOps.GateSymbol(:y, 0, false))
```
"""
function c_y(q_target::Int, q_ctrl::Int)
    return GateCall2( GateSymbols.c_y, q_target, q_ctrl, GateSymbols.y )
end

"""
    c_z(q_target::Int, q_ctrl::Int)

Generate a controlled Pauli-x (two qubit) GateCall (GateCall2), controlled on the index `q_ctrl` applied to the target `q_target`.

# Examples
```julia-repl
julia> QXZoo.DefaultGates.c_z(0, 1)
QXZoo.GateOps.GateCall2(QXZoo.GateOps.GateSymbol(:c_z, 0, false), 0, 1, QXZoo.GateOps.GateSymbol(:z, 0, false))
```
"""
function c_z(q_target::Int, q_ctrl::Int)
    return GateCall2( GateSymbols.c_z, q_target, q_ctrl, GateSymbols.z )
end

"""
    c_u(label::GateSymbol, q_target::Int, q_ctrl::Int)

Generate a controlled unitary (two qubit) GateCall (GateCall2), controlled on the index `q_ctrl` applied to the target `q_target`.

# Examples
```julia-repl
julia> QXZoo.DefaultGates.c_u(QXZoo.GateOps.GateSymbol(:myCU), 0, 1)
QXZoo.GateOps.GateCall2(QXZoo.GateOps.GateSymbol(:myCU, 0, false), 0, 1, nothing)
```
"""
function c_u(label::GateOps.GateSymbol, q_target::Int, q_ctrl::Int)
    return GateCall2( label, q_target, q_ctrl)
end

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% #

"""
    c_r_x(q_target::Int, q_ctrl::Int, theta::Number)

Generate a controlled rotation about x (exp(-iθx/2)) GateCall (GateCall2), controlled on the index `q_ctrl` applied to the target `q_target`.

# Examples
```julia-repl
julia> QXZoo.DefaultGates.c_r_x( 0, 1, pi/2)
QXZoo.GateOps.GateCall2(QXZoo.GateOps.GateSymbolP(:c_r_x, 0, false, 1.5707963267948966), 0, 1, nothing)
```
"""
function c_r_x(q_target::Int, q_ctrl::Int, theta::Number)
    return GateCall2( GateSymbols.c_r_x(theta), q_target, q_ctrl )
end

"""
    c_r_y(q_target::Int, q_ctrl::Int, theta::Number)

Generate a controlled rotation about y (exp(iθy/2)) GateCall (GateCall2), controlled on the index `q_ctrl` applied to the target `q_target`.

# Examples
```julia-repl
julia> QXZoo.DefaultGates.c_r_y( 0, 1, pi/3)
QXZoo.GateOps.GateCall2(QXZoo.GateOps.GateSymbolP(:c_r_y, 0, false, 1.0471975511965976), 0, 1, nothing)
```
"""
function c_r_y(q_target::Int, q_ctrl::Int, theta::Number)
    return GateCall2( GateSymbols.c_r_y(theta), q_target, q_ctrl )
end

"""
    c_r_z(q_target::Int, q_ctrl::Int, theta::Number)

Generate a controlled rotation about z (exp(iθz/2)) GateCall (GateCall2), controlled on the index `q_ctrl` applied to the target `q_target`

# Examples
```julia-repl
julia> QXZoo.DefaultGates.c_r_z( 0, 1, pi/4)
QXZoo.GateOps.GateCall2(QXZoo.GateOps.GateSymbolP(:c_r_z, 0, false, 0.7853981633974483), 0, 1, nothing)
```
"""
function c_r_z(q_target::Int, q_ctrl::Int, theta::Number)
    return GateCall2( GateSymbols.c_r_z(theta), q_target, q_ctrl )
end

"""
    c_r_phase(q_target::Int, q_ctrl::Int, theta::Real)

Generate a controlled phase rotation about (controlled [1 0; 0 exp(iθ)] GateCall (GateCall2), controlled on the index `q_ctrl` applied to the target `q_target`.

# Examples
```julia-repl
julia> QXZoo.DefaultGates.c_r_phase( 0, 1, pi/7 )
QXZoo.GateOps.GateCall2(QXZoo.GateOps.GateSymbolP(:c_r_phase, 0, false, 0.4487989505128276), 0, 1, nothing)
```
"""
function c_r_phase(q_target::Int, q_ctrl::Int, theta::Number)
    return GateCall2( GateSymbols.c_r_phase(theta), q_target, q_ctrl )
end

end


text:
using ProgressMeter, CSV
function write_results(algo, prob, optval = NaN, elapsed_seconds = NaN; warnings = false, errors = false)
    row = (;algo.algo, algo.settings, prob.numeric_type, prob.problem, optval, elapsed_seconds, warnings, errors)
    CSV.write(joinpath(@__DIR__, "results.csv"), [row], append = isfile(joinpath(@__DIR__, "results.csv")))
end

isfile(joinpath(@__DIR__, "results.csv")) && error("`results.csv` already exists")

include("common.jl")

# The idea is these can be run in parallel, and the thunks that write to the CSV returned,
# then those could be executed serially.
# It's probably more robust to have each process write a separate one-line CSV
# and then collect them at the end, though.
function run_problem(algo_idx, problem_idx; verbose = true)
    algo = algos[algo_idx]
    prob = problems[problem_idx]
    out = Pipe()
    err = Pipe()
    proc = run(pipeline(`julia --project=$(@__DIR__) do_problem.jl $(algo_idx) $(problem_idx)`, stdout=out, stderr=err), wait=false)
    time_elapsed = 0.0
    if verbose
        meter = Progress(timeout*2, desc = "Problem $problem_idx/$n_problems, algorithm $algo_idx/$n_algos ")
    end
    while time_elapsed < timeout && Base.process_running(proc)
        sleep(0.5)
        verbose && next!(meter)
        time_elapsed += 0.5
    end
    verbose && finish!(meter)
    if Base.process_running(proc)
        verbose && println("Timed out!")
        Base.kill(proc)

        f = () -> write_results(algo, prob)
    else
        verbose && println("Finished!")
        close(err.in)
        close(out.in)
        errors_text = String(read(err))
        if !isempty(errors_text) && verbose
            @error errors_text
        end
        results = split(chomp(String(read(out))), '\n')
        warnings = false
        if length(results) > 2
            verbose && @warn join(results[1:end-2], "\n")
            warnings = true
        end
        verbose && @info results
        if length(results) < 2
            verbose && @error "Not enough returns, something went wrong"
            errors = true
            optval, elapsed = NaN, NaN
        else
            optval, elapsed = parse.(Float64, results[end-1:end])
        end
        f = () -> write_results(algo, prob, optval, elapsed; errors = !isempty(errors_text), warnings)
    end
    return f
end

n_problems = length(problems)
n_algos = length(algos)
timeout = 60*5
for problem_idx in 1:n_problems, algo_idx in 1:n_algos
    # To run just a subset of the problems, a filter can be put such as
    # the following, which were used when I added algorithms after the first run.
    ## algos[algo_idx].algo == "MISDP_dB" || continue
    ## algos[algo_idx].algo == "MISDP (dB^2)" || continue
    f = run_problem(algo_idx, problem_idx; verbose = true)
    f()
end

df = CSV.read(joinpath(@__DIR__, "results.csv"))


text:
export structural_weight, structural_mean

# Projection on vector field

struct ProjVectorField{T,N}<:AbstractLinearOperator{T,N,N}
    ξ::AbstractArray{T,N}
end

AbstractLinearOperators.domain_size(P::ProjVectorField) = size(P.ξ)
AbstractLinearOperators.range_size(P::ProjVectorField) = size(P.ξ)
AbstractLinearOperators.matvecprod(P::ProjVectorField{T,N}, u::AbstractArray{T,N}) where {T,N} = u-P.ξ.*ptdot(u,P.ξ)
AbstractLinearOperators.matvecprod_adj(P::ProjVectorField{T,N}, u::AbstractArray{T,N}) where {T,N} = P*u

function structural_weight(u::AbstractArray{T,N}; ∇::Union{Nothing,AbstractLinearOperator}=nothing, h::NTuple{N,S}=tuple(ones(Float32,N)...), η::U=0f0) where {T,N,S<:Real,U<:Real}
    ∇ === nothing && (∇ = gradient_operator(size(u), h; T=T); u isa CuArray && (∇ = gpu(∇)))
    ∇u = ∇*u
    return ProjVectorField{T,N+1}(∇u./ptnorm2(∇u; η=real(T)(η)))
end

Flux.gpu(P::ProjVectorField{T,N}) where {T,N} = ProjVectorField{T,N}(gpu(P.ξ))
Flux.cpu(P::ProjVectorField{T,N}) where {T,N} = ProjVectorField{T,N}(cpu(P.ξ))

# Utils

function structural_mean(u::AbstractArray{T,N}; ∇::Union{Nothing,AbstractLinearOperator}=nothing, h::NTuple{N,S}=tuple(ones(Float32,N)...)) where {T,N,S<:Real}
    ∇ === nothing && (∇ = gradient_operator(size(u), h; T=T); u isa CuArray && (∇ = gpu(∇)))
    return sum(ptnorm2(∇*u))[1]/prod(size(u))
end

max_stars_repo_name: grizzuti/FastSolversForWeightedTV

max_issues_repo_name: grizzuti/FastSolversForWeightedTV

max_forks_repo_name: grizzuti/FastSolversForWeightedTV


text:
completion_test(line, char) = LanguageServer.textDocument_completion_request(LanguageServer.CompletionParams(LanguageServer.TextDocumentIdentifier("testdoc"), LanguageServer.Position(line, char), missing), server, server.jr_endpoint)

@testset "latex completions" begin
    settestdoc("""
    \\therefor
    .\\therefor
    #\\therefor
    "\\therefor"
    \"\"\"\\therefor\"\"\"
    ^\\therefor
    \\:water_buffal
    """)
    @test completion_test(0, 9).items[1].textEdit.newText == "∴"
    @test completion_test(0, 9).items[1].textEdit.range == LanguageServer.Range(0, 0, 0, 9)

    @test completion_test(1, 10).items[1].textEdit.newText == "∴"
    @test completion_test(1, 10).items[1].textEdit.range == LanguageServer.Range(1, 1, 1, 10)

    @test completion_test(2, 10).items[1].textEdit.newText == "∴"
    @test completion_test(2, 10).items[1].textEdit.range == LanguageServer.Range(2, 1, 2, 10)

    @test completion_test(3, 10).items[1].textEdit.newText == "∴"
    @test completion_test(3, 10).items[1].textEdit.range == LanguageServer.Range(3, 1, 3, 10)

    @test completion_test(4, 12).items[1].textEdit.newText == "∴"
    @test completion_test(4, 12).items[1].textEdit.range == LanguageServer.Range(4, 3, 4, 12)

    @test completion_test(5, 10).items[1].textEdit.newText == "∴"
    @test completion_test(5, 10).items[1].textEdit.range == LanguageServer.Range(5, 1, 5, 10)

    @test completion_test(6, 14).items[1].textEdit.newText == "🐃"
    @test completion_test(6, 14).items[1].textEdit.range == LanguageServer.Range(6, 0, 6, 14)
end

@testset "path completions" begin
end

@testset "import completions" begin
    settestdoc("import Base: r")
    @test any(item.label == "rand" for item in completion_test(0, 14).items)

    settestdoc("import ")
    @test all(item.label in ("Main", "Base", "Core") for item in completion_test(0, 7).items)

    settestdoc("""module M end
    import .""")
    @test_broken completion_test(1, 8).items[1].label == "M"

    settestdoc("import Base.")
    @test any(item.label == "Meta" for item in completion_test(0, 12).items)

    settestdoc("import Base.M")
    @test any(item.label == "Meta" for item in completion_test(0, 13).items)

    settestdoc("import Bas")
    @test any(item.label == "Base" for item in completion_test(0, 10).items)
end

@testset "getfield completions" begin
    settestdoc("Base.")
    @test any(item.label == "Base" for item in completion_test(0, 5).items)

    settestdoc("Base.r")
    @test any(item.label == "rand" for item in completion_test(0, 6).items)

    settestdoc("""
    using Base.Meta
    Base.Meta.
    """)
    @test any(item.label == "quot" for item in completion_test(1, 10).items)

    settestdoc("""
    module M
    inner = 1
    end
    M.
    """)
    @test any(item.label == "inner" for item in completion_test(3, 2).items)

    settestdoc("""
    x = Expr()
    x.
    """)
    @test all(item.label in ("head", "args") for item in completion_test(1, 2).items)

    settestdoc("""
    struct T
        f1
        f2
    end
    x = T()
    x.
    """)
    @test all(item.label in ("f1", "f2") for item in completion_test(1, 2).items)
end

@testset "token completions" begin
    settestdoc("B")
    @test any(item.label == "Base" for item in completion_test(0, 1).items)

    settestdoc("r")
    @test any(item.label == "rand" for item in completion_test(0, 1).items)

    settestdoc("@t")
    @test any(item.label == "@time" for item in completion_test(0, 2).items)

    settestdoc("i")
    @test any(item.label == "if" for item in completion_test(0, 1).items)

    settestdoc("i")
    @test any(item.label == "in" for item in completion_test(0, 1).items)

    settestdoc("for")
    @test any(item.label == "for" for item in completion_test(0, 3).items)

    settestdoc("in")
    @test any(item.label == "in" for item in completion_test(0, 2).items)

    settestdoc("isa")
    @test any(item.label == "isa" for item in completion_test(0, 3).items)
end

@testset "scope var completions" begin
    settestdoc("""
    myvar = 1
    βbb = 2
    bβb = 3
    myv
    βb
    bβ
    """)
    @test any(item.label == "myvar" for item in completion_test(3, 3).items)
    @test any(item.label == "βbb" for item in completion_test(4, 2).items)
    @test any(item.label == "bβb" for item in completion_test(5, 2).items)
end


text:
#This file contains a list of the currently implemented likelihood function

abstract type Likelihood end

include("bernoulli.jl")
include("exponential.jl")
include("gaussian.jl")
include("studentT.jl")
include("poisson.jl")
include("binomial.jl")

#————————————————————————————————————————————
#Predict observations at test locations

""" Computes the predictive mean and variance given a Gaussian distribution for f using quadrature"""
function predict_obs(lik::Likelihood, fmean::AbstractVector, fvar::AbstractVector)
    n_gaussHermite = 20
    nodes, weights = gausshermite(n_gaussHermite)
    weights /= sqrtπ
    f = fmean .+ sqrt.(2*fvar)*nodes'

    mLik = Array{Float64}(undef, size(f)); vLik = Array{Float64}(undef, size(f));
    @inbounds for i in 1:n_gaussHermite
        fi = view(f, :, i)
        mLik[:, i] = mean_lik(lik, fi)
        vLik[:, i] = var_lik(lik, fi)
    end
    μ = mLik*weights
    σ² = (vLik + mLik.^2)*weights - μ.^2
    return μ, σ²
end


text:
isinteractive() && include("common.jl")

####
#### utility functions
####

"Test kinetic energy gradient by automatic differentiation."
function test_KE_gradient(κ::DynamicHMC.EuclideanKineticEnergy, p)
    ∇ = ∇kinetic_energy(κ, p)
    ∇_AD = ForwardDiff.gradient(p -> kinetic_energy(κ, p), p)
    @test ∇ ≈ ∇_AD
end

####
#### testsets
####

@testset "Gaussian KE full" begin
    for _ in 1:100
        K = rand(2:10)
        Σ = rand_Σ(Symmetric, K)
        κ = GaussianKineticEnergy(inv(Σ))
        @unpack M⁻¹, W = κ
        @test W isa LowerTriangular
        @test M⁻¹ * W * W' ≈ Diagonal(ones(K))
        m, C = simulated_meancov(()->rand_p(RNG, κ), 10000)
        @test Matrix(Σ) ≈ C rtol = 0.1
        test_KE_gradient(κ, randn(K))
    end
end

@testset "Gaussian KE diagonal" begin
    for _ in 1:100
        K = rand(2:10)
        Σ = rand_Σ(Diagonal, K)
        κ = GaussianKineticEnergy(inv(Σ))
        @unpack M⁻¹, W = κ
        @test W isa Diagonal
        # FIXME workaround for https://github.com/JuliaLang/julia/issues/28869
        @test M⁻¹ * (W * W') ≈ Diagonal(ones(K))
        m, C = simulated_meancov(()->rand_p(RNG, κ), 10000)
        @test Matrix(Σ) ≈ C rtol = 0.1
        test_KE_gradient(κ, randn(K))
    end
end

@testset "phasepoint internal consistency" begin
    # when this breaks, interface was modified, rewrite tests
    @test fieldnames(PhasePoint) == (:Q, :p)
    "Test the consistency of cached values."
    function test_consistency(H, z)
        @unpack q, ℓq, ∇ℓq = z.Q
        @unpack ℓ = H
        ℓ2, ∇ℓ2 = logdensity_and_gradient(ℓ, q)
        @test ℓ2 == ℓq
        @test ∇ℓ2 == ∇ℓq
    end
    @unpack H, z, Σ = rand_Hz(rand(3:10))
    test_consistency(H, z)
    ϵ = find_stable_ϵ(H.κ, Σ)
    for _ in 1:10
        z = leapfrog(H, z, ϵ)
        test_consistency(H, z)
    end
end

@testset "leapfrog calculation" begin
    # Simple leapfrog implementation. `q`: position, `p`: momentum, `ℓ`: neg_energy, `ϵ`:
    # stepsize. `m` is the diagonal of the kinetic energy ``K(p)=p'M⁻¹p``, defaults to `I`.
    function leapfrog_Gaussian(q, p, ℓ, ϵ, m = ones(length(p)))
        u = .√(1 ./ m)
        pₕ = p .+ ϵ/2 .* last(logdensity_and_gradient(ℓ, q))
        q′ = q .+ ϵ * u .* (u .* pₕ) # mimic numerical calculation leapfrog performs
        p′ = pₕ .+ ϵ/2 .* last(logdensity_and_gradient(ℓ, q′))
        q′, p′
    end

    n = 3
    M = rand_Σ(Diagonal, n)
    m = diag(M)
    κ = GaussianKineticEnergy(inv(M))
    q = randn(n)
    p = randn(n)
    Σ = rand_Σ(n)
    ℓ = multivariate_normal(randn(n), cholesky(Σ).L)
    H = Hamiltonian(κ, ℓ)
    ϵ = find_stable_ϵ(H.κ, Σ)
    z = PhasePoint(evaluate_ℓ(ℓ, q), p)

    @testset "arguments not modified" begin
        q₂, p₂ = copy(q), copy(p)
        q′, p′ = leapfrog_Gaussian(q, p, ℓ, ϵ, m)
        z′ = leapfrog(H, z, ϵ)
        @test p == p₂               # arguments not modified
        @test q == q₂
        @test z′.Q.q ≈ q′
        @test z′.p ≈ p′
    end

    @testset "leapfrog steps" begin
        for i in 1:100
            q, p = leapfrog_Gaussian(q, p, ℓ, ϵ, m)
            z = leapfrog(H, z, ϵ)
            @test z.Q.q ≈ q
            @test z.p ≈ p
        end
    end

    @testset "leapfrog on invalid values" begin
        p = randn(n)
        q = fill(NaN, n)
        Q = evaluate_ℓ(ℓ, q)
        err = ArgumentError("Internal error: leapfrog called from non-finite log density")
        @test_throws err leapfrog(H, PhasePoint(Q, p), 0.01)
    end
end

@testset "leapfrog Hamiltonian invariance" begin
    "Test that the Hamiltonian is invariant using the leapfrog integrator."
    function test_hamiltonian_invariance(H, z, L, ϵ; atol)
        π₀ = logdensity(H, z)
        warned = false
        for i in 1:L
            z = leapfrog(H, z, ϵ)
            Δ = logdensity(H, z) - π₀
            if abs(Δ) ≥ atol && !warned
                @warn "Hamiltonian invariance violated: step $(i) of $(L), Δ = $(Δ)"
                show(H)
                show(z)
                warned = true
            end
            @test Δ ≈ 0 atol = atol
        end
    end

    for _ in 1:100
        @unpack H, z = rand_Hz(rand(2:5))
        ϵ = find_initial_stepsize(InitialStepsizeSearch(), H, z)
        test_hamiltonian_invariance(H, z, 100, ϵ/100; atol = 0.5)
    end
end

@testset "leapfrog back and forth" begin
    for _ in 1:1000
        @unpack H, z = rand_Hz(5)
        z1 = z
        N = 5
        ϵ = 0.1
        z1 = leapfrog(H, z1, ϵ)
        z1 = leapfrog(H, z1, -ϵ)
        @test z.p ≈ z1.p norm = x -> norm(x, Inf)
        @test z.Q.q ≈ z1.Q.q norm = x -> norm(x, Inf)
    end

    for _ in 1:100
        @unpack H, z = rand_Hz(2)
        z1 = z
        N = 3
        ϵ = 0.1

        # forward
        for _ in 1:N
            z1 = leapfrog(H, z1, ϵ)
        end

        # backward
        for _ in 1:N
            z1 = leapfrog(H, z1, -ϵ)
        end

        @test z.p ≈ z1.p norm = x -> norm(x, Inf) rtol = 0.001
        @test z.Q.q ≈ z1.Q.q norm = x -> norm(x, Inf) rtol = 0.001
    end
end

@testset "PhasePoint building blocks and infinite values" begin
    # wrong gradient length
    @test_throws ArgumentError EvaluatedLogDensity([1.0, 2.0], 1.0, [1.0])

    # wrong p length
    Q = EvaluatedLogDensity([1.0, 2.0], 1.0, [1.0, 2.0])
    @test_throws ArgumentError PhasePoint(Q, [1.0])
    @test PhasePoint(Q, [1.0, 2.0]) isa PhasePoint

    # fallback constructors
    Q1 = EvaluatedLogDensity([1.0, 2.0], -2.0, [3.0, 3.0])       # standard
    Q2 = EvaluatedLogDensity([1, 2], -2.0, [3.0, 3.0])           # promote
    Q3 = EvaluatedLogDensity((i for i in 1:2), -2.0, [3.0, 3.0]) # generator
    @test Q1.q == Q2.q == Q3.q
    @test Q1.ℓq == Q2.ℓq == Q3.ℓq
    @test Q1.∇ℓq == Q2.∇ℓq == Q3.∇ℓq

    # infinity fallbacks
    h = Hamiltonian(GaussianKineticEnergy(1), multivariate_normal(zeros(1)))
    @test logdensity(h, PhasePoint(EvaluatedLogDensity([1.0], -Inf, [1.0]), [1.0])) == -Inf
    @test logdensity(h, PhasePoint(EvaluatedLogDensity([1.0], NaN, [1.0]), [1.0])) == -Inf
    @test logdensity(h, PhasePoint(EvaluatedLogDensity([1.0], 9.0, [1.0]), [NaN])) == -Inf

end

@testset "Hamiltonian and KE printing" begin
    κ = GaussianKineticEnergy(Diagonal([1.0, 4.0]))
    @test repr(κ) == "Gaussian kinetic energy (Diagonal), √diag(M⁻¹): [1.0, 2.0]"
    H = Hamiltonian(κ, multivariate_normal(zeros(2)))
    @test repr(H) ==
        "Hamiltonian with Gaussian kinetic energy (Diagonal), √diag(M⁻¹): [1.0, 2.0]"
    @test_throws ArgumentError Hamiltonian(κ, multivariate_normal(zeros(1)))
end

####
#### test Hamiltonian/leapfrog using HMC
####

"""
$(SIGNATURES)

Simple Hamiltonian Monte Carlo transition, for testing.
"""
function HMC_transition(H, z::PhasePoint, ϵ, L)
    π₀ = logdensity(H, z)
    z′ = z
    for _ in 1:L
        z′ = leapfrog(H, z′, ϵ)
    end
    Δ = logdensity(H, z′) - π₀
    accept = Δ > 0 || (rand() < exp(Δ))
    accept ? z′ : z
end

"""
$(SIGNATURES)

Simple Hamiltonian Monte Carlo sample, for testing.
"""
function HMC_sample(H, q, N, ϵ; L = 10)
    qs = similar(q, N, length(q))
    for i in 1:N
        z = PhasePoint(evaluate_ℓ(H.ℓ, q), rand_p(RNG, H.κ))
        q = HMC_transition( H, z, ϵ, L).Q.q
        qs[i, :] = q
    end
    qs
end

@testset "unit normal simple HMC" begin
    # Tests the leapfrog and Hamiltonian code with HMC.
    K = 2
    ℓ = multivariate_normal(zeros(K), Diagonal(ones(K)))
    q = randn(K)
    H = Hamiltonian(GaussianKineticEnergy(Diagonal(ones(K))), ℓ)
    qs = HMC_sample(H, q, 10000, find_stable_ϵ(H.κ, Diagonal(ones(K))) / 5)
    m, C = mean_and_cov(qs, 1)
    @test vec(m) ≈ zeros(K) atol = 0.1
    @test C ≈ Matrix(Diagonal(ones(K))) atol = 0.1
end

@testset "normal NUTS HMC transition mean and cov" begin
    # A test for sample_tree with a fixed ϵ and κ, which is perfectly adapted and should
    # provide excellent mixing
    for _ in 1:10
        K = rand(2:8)
        N = 10000
        μ = randn(K)
        Σ = rand_Σ(K)
        L = cholesky(Σ).L
        ℓ = multivariate_normal(μ, L)
        Q = evaluate_ℓ(ℓ, randn(K))
        H = Hamiltonian(GaussianKineticEnergy(Σ), ℓ)
        qs = Array{Float64}(undef, N, K)
        ϵ = 0.5
        algorithm = NUTS()
        for i in 1:N
            Q = first(sample_tree(RNG, algorithm, H, Q, ϵ))
            qs[i, :] = Q.q
        end
        m, C = mean_and_cov(qs, 1)
        @test vec(m) ≈ μ atol = 0.1 rtol = maximum(diag(C))*0.02 norm = x -> norm(x,1)
        @test cov(qs, dims = 1) ≈ L*L' atol = 0.1 rtol = 0.1
    end
end


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

name = "libbluray"
version = v"1.3.0"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://code.videolan.org/videolan/libbluray/-/archive/1.3.0/libbluray-1.3.0.tar.gz", "b305e453f8c6c409465a2d970857ff47bddaa0906486d66b4b8c725a404bbc69")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/libbluray-*
apk add apache-ant
./bootstrap
./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target} --disable-bdjava-jar
make -j${nproc}
make install
install_license COPYING
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms =  filter!(!Sys.iswindows, supported_platforms())

# The products that we will ensure are always built
products = [
    LibraryProduct("libbluray", :libbluray),
    ExecutableProduct("bd_list_titles", :bd_list_titles),
    ExecutableProduct("bd_splice", :bd_splice),
    ExecutableProduct("bd_info", :bd_info)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="XML2_jll", uuid="02c8fc9c-b97f-50b9-bbe4-9be30ff0a78a"))
    Dependency(PackageSpec(name="FreeType2_jll", uuid="d7e528f0-a631-5988-bf34-fe36492bcfd7"))
    Dependency(PackageSpec(name="Fontconfig_jll", uuid="a3f928ae-7b40-5064-980b-68af3947d34b"))
    Dependency(PackageSpec(name="libudfread_jll", uuid="037e6697-03b9-52b7-b841-7aee0d773eb5"))
]

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

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
# Tests in-place products with the implicit Given's rotations we have by forming
# the Given's rotation as a matrix explicitly.

using Test, LinearAlgebra
using ArnoldiMethod: Hessenberg, Rotation2, Rotation3

@testset "Givens rotation" begin
    @testset "Single rotation" begin
        @testset "lmul!" for T in (Float64, ComplexF64)
            A = rand(T, 6, 5)

            G = Rotation2(rand(real(T)), rand(T), 2)
            G_mat = Matrix(G, 6)

            @test [A[:,1:1] G_mat * A[:,2:4] A[:,5:5]] ≈ lmul!(G, copy(A), 2, 4)
            @test G_mat * A ≈ lmul!(G, copy(A))
        end

        @testset "rmul!" for T in (Float64, ComplexF64)
            A = rand(T, 10, 5)

            G = Rotation2(rand(real(T)), rand(T), 2)
            G_mat = Matrix(G, 5)

            @test [A[1:1,:]; A[2:4,:]*G_mat'; A[5:10,:]] ≈ rmul!(copy(A), G, 2, 4)
            @test A * G_mat' ≈ rmul!(copy(A), G)
        end
    end

    @testset "Double rotation" begin
        @testset "lmul!" for T in (Float64, ComplexF64)
            A = rand(T, 6, 5)

            G = Rotation3(rand(real(T)), rand(T), rand(real(T)), rand(T), 2)
            G_mat = Matrix(G, 6)

            @test [A[:,1:1] G_mat * A[:,2:4] A[:,5:5]] ≈ lmul!(G, copy(A), 2, 4)
            @test G_mat * A ≈ lmul!(G, copy(A))
        end

        @testset "rmul!" for T in (Float64, ComplexF64)
            A = rand(T, 10, 5)

            G = Rotation3(rand(real(T)), rand(T), rand(real(T)), rand(T), 2)
            G_mat = Matrix(G, 5)

            @test [A[1:1,:]; A[2:4,:]*G_mat'; A[5:10,:]] ≈ rmul!(copy(A), G, 2, 4)
            @test A * G_mat' ≈ rmul!(copy(A), G)
        end
    end
end


text:
module Throttles

abstract type Throttle end

"""
    IntervalThrottle(interval)

Instruct [`with_preview`](@ref) etc. to call `on_preview` only if the
previous invocation of `on_preview` is at least `interval` seconds ago.
"""
struct IntervalThrottle <: Throttle
    interval::Float64
end

init(::IntervalThrottle) = time_ns()
step(th::IntervalThrottle, prev) =
    if (time_ns() - prev) / 1e9 > th.interval
        (true, time_ns())
    else
        (false, prev)
    end

end # module Throttles


text:
using Pkg

export haspkg

"""
    haspkg(pkgname::String)::Bool

Determines if the package `pkgname` available in the current environment.
"""
haspkg(pkgname::String) = haskey(Pkg.installed(), pkgname)

max_stars_repo_name: akshaysridhar/CLIMA

max_issues_repo_name: akshaysridhar/CLIMA

max_forks_repo_name: akshaysridhar/CLIMA


text:
using HTTP, Sockets

const HOST = ip"0.0.0.0"
const PORT = 9999

router = HTTP.Router()
server = HTTP.Server(router)

HTTP.register!(router, "/", HTTP.HandlerFunction(req -> HTTP.Messages.Response(200, "Hello World")))
HTTP.register!(router, "/bye", HTTP.HandlerFunction(req -> HTTP.Messages.Response(200, "Bye")))
HTTP.register!(router, "*", HTTP.HandlerFunction(req -> HTTP.Messages.Response(404, "Not found")))

HTTP.serve(server, HOST, PORT)

max_stars_repo_name: PacktPublishing/Julia-Projects

max_issues_repo_name: PacktPublishing/Julia-Projects

max_forks_repo_name: PacktPublishing/Julia-Projects


text:
import accrue.cryptoerase.runtime.Condition;

class C {
    static final Condition{L /c T} c = new Condition();
    int{L /c H} x = 0;

    void put(int v) {
	this.x = v;
    }

    int get() {
	return this.x;
    }

    public static void main(String[] args) {
	C o1 = new C();
	C o2 = new C();
	o1.put(42);
	o2.put(o1.get());
    }
}

max_stars_repo_name: HarvardPL/cryptoerase

max_issues_repo_name: HarvardPL/cryptoerase

max_forks_repo_name: HarvardPL/cryptoerase


text:
using LinearAlgebra
using StaticArrays
using PyPlot
const Vec2f = SVector{2, Float64}
const Vec4f = SVector{4, Float64}

@inline function bisection_newton(f, df, left::Float64, right::Float64, eps=0.05, itr_max=20)::Float64
    x_est = (left + right)*0.5
    for itr = 1:itr_max
        f_df = f(x_est)/df(x_est)
        if (right-(x_est-f_df))*((x_est-f_df)-left)<0.0 || abs(f_df)<(right-left)/4.0
            f(x_est)>0 ? right = x_est : left = x_est
            f_df = x_est - (right + left)*0.5
        end
        x_est -= f_df
        abs(f_df)<eps && break
    end
    return x_est
end

@inline function cost_optimal(s0::Vec4f, s1::Vec4f)
    @views a = Vec2f(s0[1:2])
    @views b = Vec2f(s0[3:4])
    @views c = Vec2f(s1[1:2])
    @views d =Vec2f(s1[3:4])
    return cost_optimal(a, b, c, d)
end

@inline function cost_optimal(x0::Vec2f, v0::Vec2f, x1::Vec2f, v1::Vec2f)
    x01 = x1 - x0
    v01 = v1 - v0
    p = -4*(dot(v0, v0)+dot(v1, v1)+dot(v0, v1))
    q = 24*dot(v0+v1, x01)
    r = -36*dot(x01, x01)
    cost(t) = t + dot(v01, 4.0*v01/t-6(-v0*t+x01)/t^2)+dot(-6*v01/t^2+12(x01-v0*t)/t^3, -v0*t+x01)
    @fastmath d_cost(t) = t^4+p*t*t+q*t+r # df(t)/dt
    @fastmath dd_cost(t) = 4.0*t*t*t+2*p*t+q # ddf(t)/dt^2

    t_min = 0.0
    t_max = 10.0
    t_star = bisection_newton(d_cost, dd_cost, t_min, t_max)
    return cost(t_star), t_star
end

# see ICRA paper: D.J.Webb et al Kinodynamic RRT* (2013)
function forward_reachable_box(x0::Vec2f, v0::Vec2f, r::Float64)
    @fastmath tau_x_plus = 2/3*(-(v0.^2).+r + v0.*sqrt.((v0.^2).+r))
    @fastmath tau_x_minus = 2/3*(-(v0.^2).+r - v0.*sqrt.((v0.^2).+r))
    @fastmath xmax = v0.*tau_x_plus + x0 + sqrt.(1/3*(tau_x_plus.^2).*(-tau_x_plus.+r))
    @fastmath xmin = v0.*tau_x_minus + x0 - sqrt.(1/3*(tau_x_minus.^2).*(-tau_x_minus.+r))

    @fastmath tau_v_plus = 0.5*r
    @fastmath vmax = v0 .+ sqrt.(tau_v_plus.*(-tau_v_plus.+r))
    @fastmath vmin = v0 .- sqrt.(tau_v_plus.*(-tau_v_plus.+r))
    return xmin, xmax, vmin, vmax
end

function backward_reachable_box(x0::Vec2f, v0::Vec2f, r::Float64)
    @fastmath tau_x_plus = 2/3*((v0.^2).-r + v0.*sqrt.((v0.^2).+r))
    @fastmath tau_x_minus = 2/3*((v0.^2).-r - v0.*sqrt.((v0.^2).+r))
    @fastmath xmax = v0.*tau_x_plus + x0 + sqrt.(1/3*(tau_x_plus.^2).*(tau_x_plus.+r))
    @fastmath xmin = v0.*tau_x_minus + x0 - sqrt.(1/3*(tau_x_minus.^2).*(tau_x_minus.+r))

    @fastmath tau_v_plus = 0.5*r
    @fastmath vmax = v0 .+ sqrt.(tau_v_plus.*(tau_v_plus.+r))
    @fastmath vmin = v0 .- sqrt.(tau_v_plus.*(tau_v_plus.+r))
    return xmin, xmax, vmin, vmax
end

"""
function filter_freachable_exact(s_set, s_c, r)
    s_set_filtered = Vec4f[]
    for s in s_set
        @views xq = Vec2f(s[1:2])
        @views vq = Vec2f(s[3:4])
        @views xc = Vec2f(s_c[1:2])
        @views vc = Vec2f(s_c[3:4])
        cost_optimal(xc, vc, xq, vq)<r && push!(s_set_filtered, s)
    end
    return s_set_filtered
end
"""

function filter_reachable(Sset::Vector{Vec4f}, idxset::Vector{Int64},
                          s_c::Vec4f, r::Float64, ForR::Symbol)
    s_set_filtered = Vec4f[]
    @views x_c = Vec2f(s_c[1:2])
    @views v_c = Vec2f(s_c[3:4])
    xmin, xmax, vmin, vmax = ((ForR==:F) ?
                              forward_reachable_box(x_c, v_c, r) :
                              backward_reachable_box(x_c, v_c, r))

    function isinside(s::Vec4f)
        @inbounds !(xmin[1]<s[1]<xmax[1]) && return false
        @inbounds !(xmin[2]<s[2]<xmax[2]) && return false
        @inbounds !(vmin[1]<s[3]<vmax[1]) && return false
        @inbounds !(vmin[2]<s[4]<vmax[2]) && return false
        return true
    end
    idx_filter = Int64[]
    dist_filter = Float64[]
    time_filter = Float64[]
    for idx in idxset
        @inbounds s = Sset[idx]
        if isinside(s)
            cost, tau = (ForR==:B) ? cost_optimal(s, s_c) : cost_optimal(s_c, s)
            if cost < r
                push!(idx_filter, idx)
                push!(dist_filter, cost)
                push!(time_filter, tau)
            end
        end
    end
    return idx_filter, dist_filter, time_filter
end

# see ICRA paper: D.J.Webb et al Kinodynamic RRT* (2013), Eq.(20)
function gen_trajectory(s0::Vec4f, s1::Vec4f, tau, N_split = 10)
    @views x0 = Vec2f(s0[1:2])
    @views v0 = Vec2f(s0[3:4])
    @views x1 = Vec2f(s1[1:2])
    @views v1 = Vec2f(s1[3:4])
    x01 = x1 - x0
    v01 = v1 - v0
    d = Vec4f([-6*v01/tau^2 + 12*(-tau*v0+x01)/tau^3;
               4*v01/tau - 6*(-tau*v0+x01)/tau^2])
    @inline function f(t)
        s = t - tau
        eye = Matrix{Int}(I, 2, 2)
        @fastmath M_left = vcat(hcat(eye, eye*s), hcat(eye*0, eye))
        @fastmath M_right = vcat(hcat(eye*(-s^3)*(1.0/6.0), eye*s^2*0.5),
                       hcat(eye*(-s^2*0.5), eye*s))
        return  Vec4f(M_left*s1 + M_right*d)
    end
    waypoints = Vec4f[]
    for n in 1:N_split+1
        t = (n-1)*tau/N_split
        push!(waypoints, f(t))
    end
    return waypoints
end

function show_trajectory(s0, s1, tau, N_split = 20, c_ = :gray, linewidth_ = 0.5)
    waypoints =  gen_trajectory(s0, s1, tau, N_split)
    M = zeros(4, N_split+1)
    for i in 1:N_split+1
        @inbounds M[:, i] = waypoints[i]
    end
    @views plot(M[1, :], M[2, :], c=c_, linewidth = linewidth_)
end

function test()
    N = 10^4
    s_set = Vec4f[]
    for i = 1:N
        push!(s_set, Vec4f(rand()-0.5, rand()-0.5, 3*rand()-0.5, 3*rand()-0.5))
    end
    s_c = Vec4f(+0, +0, 0.2, 0.2)
    idxset = [x for x in 1:N]
    idx_for, trash, _ = filter_reachable(s_set, idxset, s_c, 1.2, :F)
    idx_back, trash = filter_reachable(s_set, idxset, s_c, 1.2, :B)
    s_for = s_set[idx_for]
    s_back = s_set[idx_back]

    setA = zeros(2, length(idx_for))
    setB = zeros(2, length(idx_back))
    for i = 1:length(idx_for)
        setA[:, i] = s_for[i][1:2]
    end
    for i = 1:length(idx_back)
        setB[:, i] = s_back[i][1:2]
    end
    scatter(setA[1, :], setA[2, :], s=3)
    scatter(setB[1, :], setB[2, :], s=3)

    xlim(-0.6, 0.6)
    ylim(-0.6, 0.6)
    return setA
end

function test_connection()
    """
    s0 = Vec4f(0.17, 0.08, 0.19, -0.06)
    s1 = Vec4f(0.22, 0.03, 0.139, -0.2719)
    """
    s0 = Vec4f(0.0, 0.0, 0.3, 0.0)
    s1 = Vec4f(0.1, 0.1, -0.1, -0.01)
    show_trajectory(s0, s1, 0.5)
end
#test_connection()


text:
module Histograms

export Bin, Histogram

using Statistics: quantile

"""
    Bin(lb, ub; lbf = (>=), ubf = (<))

# Description
A bin specifying a lowerbound, upperbound, and optional comparison operators.

Passing a value to an instantiated `Bin` returns a logical value indicating whether the
value falls within the interval.

# Examples
```jldoctest
julia> b = Bin(1, 20.5)
Bin(1, 20.5, >=, <)

julia> b(2)
true

julia> b.([0, 1, 20.5, 30])
4-element BitArray{1}:
 0
 1
 0
 0

julia> Bin(0.5, 0.8, lbf=(>), ubf=(<=))
Bin(0.5, 0.8, >, <=)
```
"""
struct Bin
    lb
    ub
    lbf::Function
    ubf::Function
end

Bin(lb, ub; lbf=(>=), ubf=(<)) = Bin(lb, ub, lbf, ubf)
(b::Bin)(v) = b.lbf(v, b.lb) & b.ubf(v, b.ub)

"""
    Histogram([Bin(1, 2), Bin(3, 4)])
    Histogram(randn(1000), 10, lower_tail=true, upper_tail=true; nargs...)

A set of contiguous `Bin` that describes a distribution.

Can be instantiated manually by passing a `Vector{Bin}`, or by passing a data vector
and the number of desired bins. Cut points are determined via `Statistics.quantile` and
by default extends the tails out to infinity. Default behavior will not include positive Inf.

Passing a vector to an instantiated `Histogram` returns a frequency vector indicating the percent
of values that fall into each positional `Bin`.

# Examples

```jldoctest
julia> h = Histogram([Bin(1, 5), Bin(5, 10)])
Histogram(Bin[Bin(1, 5, >=, <), Bin(5, 10, >=, <)])

julia> h[1]
Bin(1, 5, >=, <)

julia> h([2, 2, 6, 9])
2-element Array{Float64,1}:
 0.5
 0.5

julia> h = Histogram(randn(500), 4)
Histogram(Bin[Bin(-Inf, -0.6970630349428795, >=, <), Bin(-0.6970630349428795, 0.059893308055948793, >=, <), Bin(0.059893308055948793, 0.612623377088124, >=, <), Bin(0.612623377088124, 4.053317940567178, >=, <), Bin(4.053317940567178, Inf, >=, <)])

julia> h(randn(500))
5-element Array{Float64,1}:
 0.252
 0.282
 0.21
 0.256
 0.0
```
"""
struct Histogram
    bins::Vector{Bin}
end

function Histogram(v, n, lower_tail=true, upper_tail=true; nargs...)
    @assert n > 1

    q = quantile(v, collect(1:n) ./ n)
    if lower_tail q=vcat(-Inf, q) end
    if upper_tail q=vcat(q, Inf) end

    bins = Vector{Bin}(undef, length(q)-1)
    @inbounds for i in 2:length(q)
        bins[i-1] = Bin(q[i-1], q[i]; nargs...)
    end

    Histogram(bins)
end

Base.length(h::Histogram) = length(h.bins)
Base.getindex(h::Histogram, i) = h.bins[i]

function (h::Histogram)(v)
    k = length(h)
    n = length(v)
    freq = Vector{Float64}(undef, k)
    @inbounds for i in 1:k
        freq[i] = sum(h[i].(v))/n
    end
    freq
end

end # module


text:
function is_pull_required(gl::GitLink)
    gl_repo = repo_dir(gl)
    url = remote_url(gl)
    rhash = _remote_HEAD_hash(url)
    chash = _HEAD_hash(gl_repo)
    return rhash != chash
end


text:
using Plots
using DifferentialEquations

@userplot LorenzPlot
@recipe function f(lp::LorenzPlot)
    x, y, z, i = lp.args
    n = length(x)
    inds = circshift(1:n, 1-i)
    seriesalpha --> range(0, 1, length=n)
    aspect_ratio --> 1
    label --> false
    x[inds], y[inds], z[inds]
end

function anim()

    function lorenz(du, u, p, t)

        du[1] = p[1]*(u[2] - u[1])
        du[2] = u[1]*(p[2] - u[3]) - u[2]
        du[3] = u[1]*u[2] - p[3]*u[3]

    end

    u0 = [1, 5, 10]
    tspan = (0, 100)
    p = (10, 28, 8/3)

    prob = ODEProblem(lorenz, u0, tspan, p, reltol=1e-10, abstol=1e-10)
    sol = solve(prob)
    n = length(sol.t)

    plt = plot3d(1, xaxis=("x", (-30, 30)), yaxis=("y", (-30, 30)), zaxis=("z", (0, 60)),
                 title="Lorenz Attractor", marker=2, markercolor="black", dpi=300)

    # anim = @animate for i in 1:n
    #     push!(plt, sol[1,i], sol[2,i], sol[3,i])
    # end
    # gif(anim, "de_anim.gif", fps=30)

    anim = @animate for i in 1:10:n
        lorenzplot(sol[1,:], sol[2,:], sol[3,:], i, c = :thermal)
    end
    gif(anim, "de_anim_60.gif", fps=60)

end

anim()

    # anim = @animate for i in 1:10:n
    #     lorenzplot(sol[1,:], sol[2,:], sol[3,:], i)
    # end
    # gif(anim, "de_anim60.gif", fps=60)

max_stars_repo_name: nunesmelo/djs-office-hours

max_issues_repo_name: nunesmelo/djs-office-hours

max_forks_repo_name: nunesmelo/djs-office-hours


text:
#!julia --color=yes
push!(LOAD_PATH,"../src/")

using Documenter, Pipelines

makedocs(
    sitename="Pipelines.jl",
    pages = [
        "Home" => "index.md",
        "Manual" => [
            "command_program.md",
            "julia_program.md",
            "command_dependency.md",
        ],
        "API.md"
    ]
)

deploydocs(
    repo = "github.com/cihga39871/Pipelines.jl.git",
    devbranch = "main"
)


text:
using NMRPulses
using Test

id = IdlePulse(1, 1e-3)
Ω = 2*pi
pi_2 = PiHalfPulse(Ω, 0, 1e-3)
pi_pulse = PiPulse(Ω, 0, 1e-3)

s = Sequence([pi_2, id, pi_pulse, id, pi_2])

a, b, c = prepare_sequence(s)

@testset "NMRPulses.jl" begin
    # Write your tests here.
    @test s(0.0)[1] == a[1]
    @test pi_2(0.0)[1] == a[1]

end