text:
using NBodySimulator
using StochasticDiffEq

const T = 370 # °K
const T0 = 300 # °K
const kb = 8.3144598e-3 # kJ/(K*mol)
const ϵOO = 0.1554253*4.184 # kJ
const σOO = 0.3165492 # nm
const ρ = 997/1.6747# Da/nm^3
const mO = 15.999 # Da
const mH = 1.00794 # Da
const mH2O = mO+2*mH
const N = 216#floor(Int, ρ * L^3 / m)
const L = (mH2O*N/ρ)^(1/3)
const R = 0.9 # ~3*σOO
const Rel = 0.49*L
const v_dev = sqrt(kb * T /mH2O)
const τ = 0.5e-4 # ps
const t1 = 0τ
const t2 = 40000τ # ps
const k_bond = 1059.162*4.184*1e2 # kJ/(mol*nm^2)
const k_angle = 75.90*4.184 # kJ/(mol*rad^2)
const rOH = 0.1012 # nm
const ∠HOH = 113.24*pi/180 # rad
const qH = 0.41
const qO = -0.82
const k = 138.935458 #
#bodies = load_water_molecules_from_pdb("C:/Users/Michael/Desktop/GSoC/pdbs/output_4.pdb")
bodies = generate_bodies_in_cell_nodes(N, mH2O, v_dev, L)
jl_parameters = LennardJonesParameters(ϵOO, σOO, R)
e_parameters = ElectrostaticParameters(k, Rel)
spc_paramters = SPCFwParameters(rOH, ∠HOH, k_bond, k_angle)
pbc = CubicPeriodicBoundaryConditions(L)
water = WaterSPCFw(bodies, mH, mO, qH, qO,  jl_parameters, e_parameters, spc_paramters);
#thermostat = BerendsenThermostat(T0, 200τ)
thermostat = NoseHooverThermostat(T0, 200τ)
#thermostat = AndersenThermostat(90, 0.01/τ)
#thermostat = LangevinThermostat(T0, 75)
simulation = NBodySimulation(water, (t1, t2), pbc, thermostat, kb);
result = @time run_simulation(simulation, VelocityVerlet(), dt=τ)
#result = @time run_simulation(simulation, EM(),  dt=τ)

#time_now = Dates.format(now(), "yyyy_mm_dd_HH_MM_SS")
#Nactual = length(bodies)
#timesteps = round(length(result.solution.t))

(rs, grf) = @time rdf(result)
(ts, dr2) = @time msd(result)

#using JLD
#save("D:/water $Nactual molecules $timesteps steps $time_now.jld", "rs", rs, "grf", grf, "ts", ts, "dr2", dr2)
#save("D:/water $Nactual molecules $timesteps steps $time_now.jld", "rs", rs, "grf", grf, "ts", ts, "dr2", dr2, "e_tot", e_tot, "e_kin", e_kin, "e_pot", e_pot)
#save("D:/!_good_water $Nactual molecules $timesteps steps $time_now.jld", "rs", rs, "grf", grf, "ts", ts, "dr2", dr2, "τ", τ, "t1", t1, "t2", t2, "N", N, "L", L)
#@time save_to_pdb(result, "D:/water simulation $Nactual molecules and $timesteps steps $time_now.pdb" )

#using Plots
#import GR
#plot(rs, grf, xlim=[0, 0.4999L], label=["Radial distribution function"],ylabel="g(r)", xlabel="r, nm")

#@time animate(result, "D:/$Nactual H20 particles with $timesteps timesteps $time_now.gif")

#=
t = t1:τ:result.solution.t[end-1]
temper = temperature.(result, t)

using JLD
save("d:/$(typeof(thermostat)) thermostat for water $T0 _$time_now.jld", "t",t, "temper", temper, "T0", T0, "thermostat", thermostat)

using Plots
pl=plot(t, temper, ylim=[0,500], xlabel="t, ps", ylabel = "T, °K", label="Temperature, °K", linewidth=2)
plot!(pl, t, T0*ones(length(t)), label = "$T0 °K", linewidth=2, legend = :bottomright)
#plot!(pl, title="Berendsen thermostat at tau=$(thermostat.τ) ps")
plot!(pl, title="Andersen thermostat at dt*v=$(thermostat.ν*τ) ")
#plot!(pl, title="Nose-Hoover thermostat at tau=$(thermostat.τ) ps")
=#

#(rs, grf) = @time rdf(result)
#(ts, dr2) = @time msd(result)

#=
t = t1:τ:result.solution.t[end-1]
temper = @time temperature.(result, t)

time_now = Dates.format(now(), "yyyy_mm_dd_HH_MM_SS")
Nactual = length(bodies)
timesteps = round(length(result.solution.t))

using MAT
file = matopen("d:/water $(length(bodies)) omm units at temperature $T0 and gamma $(thermostat.γ) and $timesteps timesteps with tau $τ _$time_now.mat", "w")
write(file, "t", collect(t))
write(file, "temper", temper)
if thermostat isa LangevinThermostat
    write(file, "gamma", thermostat.γ)
end
#write(file, "rs", rs)
#write(file, "ts", ts)
#write(file, "grf", grf)
#write(file, "dr2", dr2)
close(file)

@time NBodySimulator.save_to_pdb(result, "d:/water from PDB $timesteps timesteps $time_now.pdb" )
=#


text:
using CUDArt, Knet
import Base: isapprox

function isapprox(x, y;
                  maxeps::Real = max(eps(eltype(x)), eps(eltype(y))),
                  rtol::Real=maxeps^(1/4), atol::Real=maxeps^(1/2))
    size(x) == size(y) || (warn("isapprox: $(size(x))!=$(size(y))"); return false)
    x = convert(Array, x)
    y = convert(Array, y)
    d = abs(x-y)
    s = abs(x)+abs(y)
    all(d .< (atol + rtol * s))
end

ext: jl


text:
@testset "add two numbers" begin
    a = ListNode{Int}([2, 4, 3])
    b = ListNode{Int}([5, 6, 4])
    expected = ListNode{Int}([7, 0, 8])
    @test expected == add_two_numbers(a, b)
    @test add_two_numbers(ListNode{Int}([0]), ListNode{Int}([0])) == ListNode{Int}([0])
    @test add_two_numbers(
        ListNode{Int}([9, 9, 9, 9, 9, 9, 9]),
        ListNode{Int}([9, 9, 9, 9]),
    ) == ListNode{Int}([8, 9, 9, 9, 0, 0, 0, 1])
end


text:
function get_additional_state_dimension(obj::AbstractJacobiansMapper)
    return jcall(obj, "getAdditionalStateDimension", jint, ())
end

function get_name(obj::AbstractJacobiansMapper)
    return jcall(obj, "getName", JString, ())
end

function get_parameters(obj::AbstractJacobiansMapper)
    return jcall(obj, "getParameters", jint, ())
end

function get_parameters_jacobian(obj::AbstractJacobiansMapper, arg0::SpacecraftState, arg1::Vector{Vector{jdouble}})
    return jcall(obj, "getParametersJacobian", void, (SpacecraftState, Vector{Vector{jdouble}}), arg0, arg1)
end

function get_state_jacobian(obj::AbstractJacobiansMapper, arg0::SpacecraftState, arg1::Vector{Vector{jdouble}})
    return jcall(obj, "getStateJacobian", void, (SpacecraftState, Vector{Vector{jdouble}}), arg0, arg1)
end

function set_initial_jacobians(obj::AbstractJacobiansMapper, arg0::SpacecraftState, arg1::Vector{Vector{jdouble}}, arg2::Vector{Vector{jdouble}}, arg3::Vector{jdouble})
    return jcall(obj, "setInitialJacobians", void, (SpacecraftState, Vector{Vector{jdouble}}, Vector{Vector{jdouble}}, Vector{jdouble}), arg0, arg1, arg2, arg3)
end


text:
#=
The type of a patch difference.

Is the sum over the patch of point differences.
=#

PatchDifferenceType = Float32

max_stars_repo_name: bootchk/resynthesizerJulia

max_issues_repo_name: bootchk/resynthesizerJulia

max_forks_repo_name: bootchk/resynthesizerJulia


text:
# test all solvers with the well known Woods test function
nlp = woods(n = 10)

solver = ALL_solvers[25]
(x, f, gNorm, iter, optimal, tired, status) = solver(nlp, verbose = false);
@printf(
    "%-15s  %8d  %9.2e  %7.1e  %5d  %5d  %6d  %s\n",
    nlp.meta.name,
    nlp.meta.nvar,
    f,
    gNorm,
    nlp.counters.neval_obj,
    nlp.counters.neval_grad,
    nlp.counters.neval_hprod,
    status
)    #stats = run_solver(solver, model, verbose=false)
#@test (all([stats...] .>= 0))
reset!(nlp)
solver = ALL_solvers[26]
(x, f, gNorm, iter, optimal, tired, status) = solver(nlp, verbose = false);
@printf(
    "%-15s  %8d  %9.2e  %7.1e  %5d  %5d  %6d  %s\n",
    nlp.meta.name,
    nlp.meta.nvar,
    f,
    gNorm,
    nlp.counters.neval_obj,
    nlp.counters.neval_grad,
    nlp.counters.neval_hprod,
    status
)    #stats = run_solver(solver, model, verbose=false)

#reset!(nlp);


text:
function procrustes(X::Shape, Y::Shape)
    @assert size(X) == size(Y)
    n, m = size(X)
    mu_x = mean(X, 1)
    mu_y = mean(Y, 1)

    X0 = X .- mu_x
    Y0 = Y .- mu_y

    ss_x = sum(X0.^2)
    ss_y = sum(Y0.^2)

    norm_x = sqrt(ss_x)
    norm_y = sqrt(ss_y)

    X0 ./= norm_x
    Y0 ./= norm_y

    A = X0'Y0
    U, s, Vt = svd(A)
    V = Vt'
    T = V*U'

    traceTA = sum(s)

    b = traceTA * norm_x / norm_y
    d = 1 - traceTA ^ 2
    Z = norm_x * traceTA * (Y0 * T) .+ mu_x

    c = mu_x - b * (mu_y * T)

    tform = @compat(Dict("rotation" => T, "scale" => b, "translation" => c))
    return d, Z, tform

end

function align_shapes(shapes::Vector{Shape}; n_iter=100)
    shapes = copy(shapes)
    # mean_shape = mean(shapes)
    mean_shape = shapes[1]
    ref_shape = mean_shape
    for it=1:n_iter
        for i=1:length(shapes)
            _, shapes[i], _ = procrustes(mean_shape, shapes[i])
        end
        new_mean_shape = mean(shapes)
        _, Y, _ = procrustes(ref_shape, new_mean_shape)
        new_mean_shape = Y
        mean_shape = new_mean_shape
    end
    mean_shape = mean(shapes)
    return mean_shape, shapes
end


text:
@inline function add_ddfp_dd(x::Tuple{T,T}, y::T) where {T<:IEEEFloat}
    hi, lo = x
    hi, lo = two_sumof3(hi, y, lo)
    return hi, lo
end

@inline function sub_ddfp_dd(x::Tuple{T,T}, y::T) where {T<:IEEEFloat}
    hi, lo = x
    hi, lo = two_sumof3(hi, lo, -y)
    return hi, lo
end

@inline function mul_ddfp_dd(x::Tuple{T,T}, y::T) where {T<:IEEEFloat}
    hi, lo = x
    hihi, hilo = two_prod(y, hi)
    lohi, lolo = two_prod(y, lo)
    hi, lo = two_sumof4(hihi, hilo, lohi, lolo)
    return hi, lo
end

@inline function dvi_ddfp_dd(x::Tuple{T,T}, y::T) where {T<:IEEEFloat}
    yinv = inv_fp_dd(y)
    zhi, zlo = mul_dddd_dd(x, yinv)
    return zhi, zlo
end


text:
function event_occurred(obj::FieldEventDetector, arg0::FieldSpacecraftState, arg1::jboolean)
    return jcall(obj, "eventOccurred", Action, (FieldSpacecraftState, jboolean), arg0, arg1)
end

function g(obj::FieldEventDetector, arg0::FieldSpacecraftState)
    return jcall(obj, "g", RealFieldElement, (FieldSpacecraftState,), arg0)
end

function get_max_check_interval(obj::FieldEventDetector)
    return jcall(obj, "getMaxCheckInterval", RealFieldElement, ())
end

function get_max_iteration_count(obj::FieldEventDetector)
    return jcall(obj, "getMaxIterationCount", jint, ())
end

function get_threshold(obj::FieldEventDetector)
    return jcall(obj, "getThreshold", RealFieldElement, ())
end

function init(obj::FieldEventDetector, arg0::FieldSpacecraftState, arg1::FieldAbsoluteDate)
    return jcall(obj, "init", void, (FieldSpacecraftState, FieldAbsoluteDate), arg0, arg1)
end

function reset_state(obj::FieldEventDetector, arg0::FieldSpacecraftState)
    return jcall(obj, "resetState", FieldSpacecraftState, (FieldSpacecraftState,), arg0)
end


text:
include("$(@__DIR__)/../src/FiniteGroups.jl")
using .FiniteGroups, LinearAlgebra

py"""
from mpl_toolkits import mplot3d
import matplotlib as mpl
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
import matplotlib.pyplot as plt
import time
from sympy.combinatorics import Permutation
from sympy.interactive import init_printing
import copy
import scipy
from scipy.stats import ortho_group

##——————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————##
## def fundemental var and fun
sigma_0 = np.identity(2, dtype=complex);
sigma_x = np.array([[0, 1], [1, 0]], dtype=complex);
sigma_y = np.array([[0, -1j], [1j, 0]], dtype=complex);
sigma_z = np.array([[1, 0], [0, -1]], dtype=complex);

sigma_1 = sigma_x
sigma_2 = sigma_y
sigma_3 = sigma_z

def direct_sum(rep1, rep2):
    _new_matrix = np.zeros(np.add(rep1.shape,rep2.shape), dtype = complex)
    _new_matrix[:rep1.shape[0],:rep1.shape[1]] = rep1
    _new_matrix[rep1.shape[0]:,rep1.shape[1]:] = rep2
    return _new_matrix

Sigma = np.array([sigma_0, sigma_1, sigma_2, sigma_3])

#因为一些特殊的原因，这个地方的ij指标反了，由这样的函数生成出来的mass是和S = sigma_z sigma_0 sigma_0反对易的
def paulimatrix3(j,i,k):
    _matrix1 = np.kron(Sigma[i], Sigma[j])
    return np.kron(_matrix1, Sigma[k])

def paulimatrix2(i,j):
    return np.kron(Sigma[i], Sigma[j])

def chop(expr, delta=10**-5):
    return np.ma.masked_inside(expr, -delta, delta).filled(0)

##——————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————##
Mass_term = np.array([paulimatrix3(0,1,0),
                      paulimatrix3(0,1,1),
                      paulimatrix3(0,1,3),
                      paulimatrix3(0,2,2),

                      paulimatrix3(1,1,0),
                      paulimatrix3(1,1,1),
                      paulimatrix3(1,1,3),
                      paulimatrix3(1,2,2),

                      paulimatrix3(2,1,2),
                      paulimatrix3(2,2,0),
                      paulimatrix3(2,2,1),
                      paulimatrix3(2,2,3),

                      paulimatrix3(3,1,0),
                      paulimatrix3(3,1,1),
                      paulimatrix3(3,1,3),
                      paulimatrix3(3,2,2),
                      ])

def four_dim_rep_pointgroup(any_point_group_rep):
    if any_point_group_rep.shape[0] >=4:
        print('矩阵维度大于四，点群的线性表示的维度最高是3，所以出错')
        final_matrix = any_point_group_rep
    elif any_point_group_rep.shape[0] ==3:
        final_matrix = direct_sum(any_point_group_rep, np.array([[1]]))
    elif any_point_group_rep.shape[0] ==2:
        final_matrix = direct_sum(any_point_group_rep, sigma_0)
    elif any_point_group_rep.shape[0] ==1:
         final_matrix = direct_sum(any_point_group_rep, sigma_0)
         final_matrix = direct_sum(final_matrix, np.array([[1]]))
    else:
        final_matrix = np.kron(sigma_0, sigma_0)
    return final_matrix

def flavor_real_copro_rep(rep, dia_or_offdia):
    if dia_or_offdia == 1:
        _fla_rep = np.kron(sigma_0, rep)
    elif dia_or_offdia == -1:
        _fla_rep = np.kron(sigma_x, rep)
    else:
        print('错误的和S的对易关系，只能是1或者-1')
        _fla_rep = np.zeros((8,8))
    return _fla_rep

def mass_rep_fun(_PG_rep, dia_or_offdia):
    _four_dim_rep = four_dim_rep_pointgroup(_PG_rep)
    _eight_dim_copro_rep = flavor_real_copro_rep(_four_dim_rep, dia_or_offdia)
    _mass_rep = np.zeros((Mass_term.shape[0],Mass_term.shape[0]), dtype = complex)
    for j in range(Mass_term.shape[0]):
        _mass_after_transfor = np.dot(_eight_dim_copro_rep,np.dot(Mass_term[j],_eight_dim_copro_rep.conjugate().T))
        for i in range(Mass_term.shape[0]):
            _mass_rep[i,j]= np.trace(np.dot(Mass_term[i],_mass_after_transfor))

    return _mass_rep/8

def judge_comm_anticomm(m1,m2):
    if np.dot(m1,m2)-np.dot(m2,m1) == np.zeros(m1.shape):
        print
##——————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————##
##transform the matrix form to "the array form" of permutation

def tran_Mat_to_Permu(matrix):
    _permu_array_form = list(range(matrix.shape[0]))
    for i in range(matrix.shape[0]):
        vec = np.zeros(matrix.shape[0])
        vec[i] = 1
        _permu_array_form[i]=np.nonzero(np.dot(matrix, vec))[0][0]
    return _permu_array_form

def index_to_conv(list_p):
    _list = copy.deepcopy(list_p)
    for i in range(len(_list)):
        for j in range(len(_list[i])):
            _list[i][j] +=1
    return _list

##——————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————##

T_three_dim_rep = np.array([[0.5,-np.sqrt(2)/2,0.5],
                            [-np.sqrt(2)/2,0,np.sqrt(2)/2],
                            [0.5,np.sqrt(2)/2,0.5]])
# T_three_dim_rep = np.array([[0,0,1],
#                             [-1,0,0],
#                             [0,1,0]])
rep = T_three_dim_rep

##——————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————##

# opera = np.load("Oh_7.npy")
# pairingsymm = np.load("Oh_pairing_4.npy")
# mass_rep_char =np.zeros(len(opera),dtype=complex)
# mass_rep_eigval=np.zeros(len(opera),dtype=object)
# for i in range(len(opera)):
#     massrep = mass_rep_fun(opera[i],int(pairingsymm[i][0][0]))
#     mass_rep_char[i] = np.trace(massrep)
#     mass_rep_eigval[i],eigvec = np.linalg.eig(massrep)
#     print(chop(np.prod(mass_rep_eigval[i])))
#     print(mass_rep_eigval[i])
# np.save("Oh_7irrep_4Pairingsymm.npy",chop(mass_rep_char))

##——————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————

##当mass在操作下会变成另一个mass而非一组mass的线性组合时，可以化成置换操作向轮换操作的分解问题。

# def mass_invar_space1():
#     _mass_rep = chop(np.real(mass_rep_fun(rep,-1)))

#     init_printing(perm_cyclic=True,pretty_print=False)
#     p = Permutation(tran_Mat_to_Permu(_mass_rep)).full_cyclic_form
#     p_conv = index_to_conv(p)
#     print('mass的不变子空间', p_conv)
#     print('———————————————————————————————————————————————————————————')

#     for i in range(len(p)):
#         if len(p[i])>4:
#             print('mass 的最小子空间维数大于4,')
#         for k in range(len(p[i])):
#             p[i][(k+1) %len(p[i])]
#             print('第%d个不变子空间的基的变换,第%d个mass-->第%d个mass的系数是%d+%dj'%(i+1, p_conv[i][k], p_conv[i][(k+1) %len(p[i])], _mass_rep[p[i][(k+1) %len(p[i])],p[i][k]].real, _mass_rep[p[i][(k+1) %len(p[i])],p[i][k]].imag))
#         for l in range(len(p[i])-1):
#             for j in range(len(p[i])-l-1):
#                 if (np.dot(Mass_term[p[i][j]],Mass_term[p[i][j+l+1]])-np.dot(Mass_term[p[i][j+l+1]],Mass_term[p[i][j]])==np.zeros(Mass_term[p[i][j]].shape)).all():
#                     print('第%d个不变子空间里面的第%d个基和第%d个基是对易的'%(i+1,j+1,j+l+2))
#                 elif (np.dot(Mass_term[p[i][j]],Mass_term[p[i][j+l+1]])+np.dot(Mass_term[p[i][j+l+1]],Mass_term[p[i][j]])==np.zeros(Mass_term[p[i][j]].shape)).all():
#                     print('第%d个不变子空间里面的第%d个基和第%d个基是反对易的'%(i+1,j+1,j+l+2))
#                 else:
#                     print('第%d个不变子空间里面的第%d个基和第%d个基即不对易也不反对易'%(i+1,j+1,j+l+2))
#         print('———————————————————————————————————————————————————————————')

##——————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————
##当mass在操作下会变成另一个mass而非一组mass的线性组合时，可以化成置换操作向轮换操作的分解问题。

mass_rep = chop((mass_rep_fun(rep,-1)))
inv_space=[[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[]]
for i in range(mass_rep.shape[0]):
    inv_space[i].insert(len(inv_space[i]), i)
    for j in range(mass_rep.shape[1]):
        if mass_rep[i][j] != 0:
            inv_space[i].insert(len(inv_space[i]), j)

invar_space = [[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[]]
invar_space_nonoverlap=[]
invar_space_nonoverlap_before = []
invar_space_nonoverlap_after =inv_space
for overindex in range(100):
    if invar_space_nonoverlap_before == invar_space_nonoverlap_after:
        break
    else:
        invar_space_nonoverlap_before = invar_space_nonoverlap_after
        for i2 in range(len(invar_space_nonoverlap_before)):
            set_union = set(invar_space_nonoverlap_before[i2])
            for j2 in range(len(invar_space_nonoverlap_before)):
                if set_union & set(invar_space_nonoverlap_before[j2])!=set():
                    set_union = set_union | set(invar_space_nonoverlap_before[j2])

            invar_space[i2]=list(np.sort(list(set_union)))
            invar_space_nonoverlap_after = invar_space

invar_space_nonoverlap = []
for l in range(len(invar_space_nonoverlap_after)):
    if invar_space_nonoverlap_after[l] not in invar_space_nonoverlap:
        invar_space_nonoverlap.append(invar_space_nonoverlap_after[l])

p = invar_space_nonoverlap
p_conv = index_to_conv(p)
print('mass的不变子空间', p_conv)
print('———————————————————————————————————————————————————————————')
for i in range(len(p)):
    mass_invar_space_rep=np.zeros((len(p[i]),len(p[i])),dtype=float)
    for k in range(len(p[i])):
        for kk in range(len(p[i])):
            mass_invar_space_rep[k,kk] = np.real(mass_rep[p[i][k],p[i][kk]])
    print('第%d个不变子空间的基矢是mass'%i,p_conv[i],'mass不变子空间的变换矩阵是\n',mass_invar_space_rep)

    for l in range(len(p[i])-1):
        for j in range(len(p[i])-l-1):
            if (np.dot(Mass_term[p[i][j]],Mass_term[p[i][j+l+1]])-np.dot(Mass_term[p[i][j+l+1]],Mass_term[p[i][j]])==np.zeros(Mass_term[p[i][j]].shape)).all():
                print('第%d个不变子空间里面的第%d个基和第%d个基是对易的'%(i+1,j+1,j+l+2))
            elif (np.dot(Mass_term[p[i][j]],Mass_term[p[i][j+l+1]])+np.dot(Mass_term[p[i][j+l+1]],Mass_term[p[i][j]])==np.zeros(Mass_term[p[i][j]].shape)).all():
                print('第%d个不变子空间里面的第%d个基和第%d个基是反对易的'%(i+1,j+1,j+l+2))
            else:
                print('第%d个不变子空间里面的第%d个基和第%d个基即不对易也不反对易'%(i+1,j+1,j+l+2))
        print('-----------------------------------------------------')
    print('———————————————————————————————————————————————————————————')

"""


text:
using Distributed

# We have to ask for other processes in our local cluster to start
addprocs(exeflags="--project")
# and we have to load the core module--or we deadlock
@everywhere using RxJulia

struct CustomTestSet <: Test.AbstractTestSet
    description::String
    delegate::Test.DefaultTestSet
end

function indent()
    last_test_depth = get(task_local_storage(), :last_test_depth, 0)
    depth = Test.get_testset_depth()
    task_local_storage(:last_test_depth, depth)
    nl = (depth != last_test_depth) ?  "\n" : ""
    repeat("$nl  ", depth)
end

function CustomTestSet(description)
    printstyled(stdout, "$(indent())Testing $description", color=:light_blue)
    CustomTestSet(description, Test.DefaultTestSet(description))
end

function Test.record(ts::CustomTestSet, set::Test.AbstractTestSet)
    Test.record(ts.delegate, set)
end

function Test.record(ts::CustomTestSet, res::Test.Fail)
    printstyled(stdout, "!"; bold = true, color=Base.error_color())
    Test.record(ts.delegate, res)
end

function Test.record(ts::CustomTestSet, err::Test.Error)
    printstyled(stdout, "E"; bold = true, color=Base.error_color())
    Test.record(ts.delegate, err)
end

function Test.record(ts::CustomTestSet, res::Test.Pass)
    printstyled(stdout, "."; bold = true, color=:green)
    Test.record(ts.delegate, res)
end

function Test.finish(ts::CustomTestSet)
    printstyled(stdout, "DONE\n", color=:light_blue)
    Test.finish(ts.delegate)
end


text:
using Pkg; Pkg.activate("../")
using Morse
import Morse: CharDict, StrDict

function CounterVocabulary(sets::Vector)
    char2ix, tag2ix, word2ix, ctag2ix = CharDict(), StrDict(), StrDict(), StrDict()

    for (i,T) in enumerate(specialTokens)
        get!(tag2ix,T,1); get!(word2ix,T,1); get!(word2ix,T,1)
        get!(char2ix,T[1],1);
    end
    specialIndicies = (unk=1, mask=1, eow=1, bow=1)

    for data in sets, sentence in data, par::ParsedIO in sentence
        for c::Char in par.chars
            char2ix[c] = get(char2ix, c, 0) + 1
        end
        for a::Analysis in par.analyses
            for c::Char in a.lemma
                tag2ix[string(c)] = get(tag2ix, string(c), 0) + 1
            end
            ctag2ix[join(a.tags,'+')] = get(ctag2ix,join(a.tags,'+'),0) + 1
        end
        word2ix[join(par.chars)] =  get(word2ix,join(par.chars),0) + 1
    end

    for data in sets, sentence in data, par::ParsedIO in sentence
        for a::Analysis in par.analyses
            for t::String in a.tags
                tag2ix[t] = get(tag2ix,t,0) + 1
            end
        end
    end
    return (chars=char2ix, tags=tag2ix, ctags=ctag2ix, words=word2ix)
end

config = Morse.intro(ARGS)

p_cls = string("morse.Classifier_lemma_",config[:lemma],"_lang_UD-",config[:langcode],"_size_full.gen")
p_seq = string("morse.S2SContext_lemma_",config[:lemma],"_lang_UD-",config[:langcode],"_size_full.gen")

DataSet  = eval(Meta.parse(config[:dataSet]))
ModelType = eval(Meta.parse(config[:modelType]))

code = config[:langcode]
lang = CODE_TO_LANG[code]
files = Morse.dir("data","ud-treebanks-v2.1/UD_$(lang)") * "/$(code)-ud-" .* ["train","dev","test"] .* ".conllu"

parser     = Parser{DataSet}(Val(config[:version]))
trn        = parseFile(first(files); p=parser, withLemma=config[:lemma], parseAll=false)
tst        = parseFile(last(files); p=parser, withLemma=config[:lemma], parseAll=false)

trncounter = CounterVocabulary([trn])

tst_pred_seq = parseFile(p_seq; p=parser, withLemma=config[:lemma], parseAll=false)
tst_pred_cls = parseFile(p_cls; p=parser, withLemma=config[:lemma], parseAll=false)

@assert length(tst_pred_seq) == length(tst) "seq gens broken"
@assert length(tst_pred_cls) == length(tst) "classification gens broken $(length(tst_pred_cls)) != $(length(tst))"
counter = trncounter.ctags

function stats(tst,seq,cls,counter)
    totaltokens = [0,0,0]
    seqtrues = [0,0,0]
    clstrues = [0,0,0]

    for (s,sent) in enumerate(tst)
        for (i,wordio) in enumerate(sent)
            if wordio.analyses[1].isValid
                correct_analysis = join(wordio.analyses[1].tags,'+')
                cnt = get(counter,correct_analysis,0)
                seq_cnt   = join(seq[s][i].analyses[1].tags,'+')  == correct_analysis
                cls_cnt   = join(cls[s][i].analyses[1].tags,'+')  == correct_analysis

                if cnt == 0
                    totaltokens[1] += 1
                    seqtrues[1]    += seq_cnt
                    clstrues[1]    += cls_cnt
                elseif cnt < 5
                    totaltokens[2] += 1
                    seqtrues[2]    += seq_cnt
                    clstrues[2]    += cls_cnt
                else cnt >= 5
                    totaltokens[3] += 1
                    seqtrues[3]    += seq_cnt
                    clstrues[3]    += cls_cnt
                end
            end
        end
    end
    return totaltokens, seqtrues, clstrues
end

total, seqt, clst =stats(tst,tst_pred_seq,tst_pred_cls,counter)

println(uppercase(code)," & ",total[1],"/",total[2],"/",total[3]," &\t ",100clst[1]/total[1]," & ",100seqt[1]/total[1],
        " &\t ",100clst[2]/total[2]," & ",100seqt[2]/total[2],
        " &\t ",100clst[3]/total[3]," & ",100seqt[3]/total[3])

# function pred_upos_trim!(data)
#     for sent in data
#         for wordio in sent
#             analysis = wordio.analyses[1]
#             if startswith(analysis.tags[1],"Upos=Upos=")
#                 analysis.tags[1] = String(SubString(analysis.tags[1],6:length(analysis.tags[1])))
#             end
#         end
#     end
# end

#pred_upos_trim!(tst_pred_seq)
#pred_upos_trim!(tst_pred_cls)

# @show length(tst)
# @show length(tst_pred_seq)
# @show length(tst_pred_cls)
# @show last(tst_pred_seq)[1]
# @show last(tst)[1]


text:
module test_julia_rounding

using Test

@test round(pi, digits=2) == 3.14
@test round(pi, sigdigits=3) == 3.14

end # module test_julia_rounding


text:
include("Util_single_RJ.jl")
include("Util_multi.jl")

@everywhere module ReversibleJumpMultipleParticle
using StatsBase
using KernelDensity
using ReversibleJumpSingleParticle
using MultipleParticle

export resample, get_values, get_weights, get_obj_func_vals, get_stepsizes,
       normalize_weights, reset_weights, drop_zero_weights, update_weights, calc_ESS, calcPA, update_stepsize_inRKHS_all, MCstep_all_parallel
export generate_init_particles, get_model_indecies, get_subset, get_proportion_of_models, MCstep_all_rj,
       MCstep_all_parallel_rj

function generate_init_particles(N, lbs_set, ubs_set, obj_funcs, model_indecies)
    """
    N: number of samplers
    lbs_set, ubs_set:: sets of bounds. Array{Float64,1} or Array{Array{Float64, 1}, 1}
    """
    particles = Array{ReversibleParticle}(length(N))
    if typeof(lbs_set) <: Array{Array{Float64, 1},1} || typeof(ubs_set) <: Array{Array{Float64, 1},1}
        particles = [init_particle(lbs_set[model_indecies[i]], ubs_set[model_indecies[i]], obj_funcs[model_indecies[i]], 1 / N, model_indecies[i]) for i in 1:N]
    elseif typeof(lbs_set) <: Array{Float64, 1} || typeof(ubs_set) <: Array{Float64, 1}
        particles = [init_particle(lbs_set, ubs_set, obj_funcs[model_indecies[i]], 1 / N, model_indecies[i]) for i in 1:N]
    else
        ErrorException(@sprintf("TypeError lbs: %s ≠ ubs: %s",typeof(lbs_set), typeof(ubs_set)))
    end
    return particles
end

function get_model_indecies(particles)
    return [particle.model_id for particle in particles]
end

function get_subset(particles, model_id)
    model_indecies = get_model_indecies(particles)
    return particles[model_indecies .== model_id]
end

function get_proportion_of_models(particles, N)
    total = length(particles)
    model_indecies = get_model_indecies(particles)
    proportion = zeros(N)
    for id in unique(model_indecies)
        particles_subset = get_subset(particles, id)
        proportion[id] = length(particles_subset) / total
    end
    return proportion
end

function MCstep_all_rj(particles, ϵ)
    model_indecies = get_model_indecies(particles)
    flags = zeros(length(particles))
    for id in unique(model_indecies)
        particles_subset = get_subset(particles, id)
        flags_subset = zeros(length(particles_subset))

        # kernel density estimation of prior
        params = hcat(get_values(particles_subset)...)'
        prior = [kde(params[:, i], npoints=16) for i in 1:size(params)[2]]

        N = length(particles_subset)
        for i in 1:N
            particles_subset[i], flags_subset[i] = MCstep(particles_subset[i], prior, ϵ)
        end
        particles[model_indecies .== id] = particles_subset
        flags[model_indecies .== id] = flags_subset
    end
    return particles, flags
end

function MCstep_all_parallel_rj(particles, ϵ)
    """
    to fetch results of @parallel for macro from several worker process,
    SharedArray must be required.
    """
    model_indecies = get_model_indecies(particles)
    flags = zeros(length(particles))
    for id in unique(model_indecies)
        particles_subset = get_subset(particles, id)
        particles_subset, flags_subset = MCstep_all_parallel(particles_subset, ϵ)
        particles[model_indecies .== id] = particles_subset
        flags[model_indecies .== id] = flags_subset
    end

    return particles, flags
end

end


text:
# Benchmark

using Plots;
using Random;
using Statistics;
using BenchmarkTools;

Random.seed!(1234);

include("GenerateQuadraticProgram.jl");
include("LinearSystemSolversFun.jl");

NS_TO_SEC_FCTR          = 1e-9;
BYTE_TO_MEGA_BYTE_FCTR  = 2 ^ -20;

## Parameters

# Simulaion
vNumElements    = [250; 500];
vNumConstraints = [125; 250];

# problemClass    = rand(instances(ProblemClass));
problemClass    = randomQp;

## Generating Model
ρ = 1e6;
σ = 1e-6;

# Solver Parameters
numFactor       = 3;
numIterations   = 50;
ϵSolver         = 1e-6;
numItrSolver    = 1e-6;

# Benchmark
benchMarkSamples    = 100;
benchMarkEvals      = 1;
benchMarkSeconds    = 2;

# Solution Threshold
ϵSol = 1e-5;

numDims = size(vNumElements, 1);

vF      = Vector{Any}(undef, 0);
vLabel  = Vector{Any}(undef, 0);

# LaLdlt!(vX, mP, mA, vV; numFactor = 3, numIterations = 50, ϵSolver = 1e-6, numItrSolver = 1e-6, ρ = 1e6, σ = 1e-6);
# ItrSolCg!(vX, mP, mA, vV; numFactor = 3, numIterations = 50, ϵSolver = 1e-6, numItrSolver = 500, ρ = 1e6, σ = 1e-6);
# KrylovCg!(vX, mP, mA, vV; numFactor = 3, numIterations = 50, ϵSolver = 1e-6, numItrSolver = 500, ρ = 1e6, σ = 1e-6);
# KrylovCr!(vX, mP, mA, vV; numFactor = 3, numIterations = 50, ϵSolver = 1e-6, numItrSolver = 500, ρ = 1e6, σ = 1e-6);
# KrylovCgLanczos!(vX, mP, mA, vV; numFactor = 3, numIterations = 50, ϵSolver = 1e-6, numItrSolver = 500, ρ = 1e6, σ = 1e-6);
# LinOpCg!(vX, mP, mA, vV; numFactor = 3, numIterations = 50, ϵSolver = 1e-6, numItrSolver = 500, ρ = 1e6, σ = 1e-6);
# KrylovTriCg!(vX, mP, mA, vV; numFactor = 3, numIterations = 50, ϵSolver = 1e-6, numItrSolver = 500, ρ = 1e6, σ = 1e-6);
# KrylovTriMr!(vX, mP, mA, vV; numFactor = 3, numIterations = 50, ϵSolver = 1e-6, numItrSolver = 500, ρ = 1e6, σ = 1e-6);

push!(vF, (vX, mP, mA, vV) -> LaLdlt!(vX, mP, mA, vV, ρ, σ; numFactor = numFactor, numIterations = numIterations, ϵSolver = ϵSolver, numItrSolver = numItrSolver));
push!(vLabel, "LDL Factorization");
push!(vF, (vX, mP, mA, vV) -> ItrSolCg!(vX, mP, mA, vV, ρ, σ; numFactor = numFactor, numIterations = numIterations, ϵSolver = ϵSolver, numItrSolver = numItrSolver));
push!(vLabel, "IterativeSolvers CG");
push!(vF, (vX, mP, mA, vV) -> KrylovCg!(vX, mP, mA, vV, ρ, σ; numFactor = numFactor, numIterations = numIterations, ϵSolver = ϵSolver, numItrSolver = numItrSolver));
push!(vLabel, "Krylov CG");
push!(vF, (vX, mP, mA, vV) -> KrylovCr!(vX, mP, mA, vV, ρ, σ; numFactor = numFactor, numIterations = numIterations, ϵSolver = ϵSolver, numItrSolver = numItrSolver));
push!(vLabel, "Krylov CR");
push!(vF, (vX, mP, mA, vV) -> KrylovCgLanczos!(vX, mP, mA, vV, ρ, σ; numFactor = numFactor, numIterations = numIterations, ϵSolver = ϵSolver, numItrSolver = numItrSolver));
push!(vLabel, "Kryloc Lanczos");
push!(vF, (vX, mP, mA, vV) -> LinOpCg!(vX, mP, mA, vV, ρ, σ; numFactor = numFactor, numIterations = numIterations, ϵSolver = ϵSolver, numItrSolver = numItrSolver));
push!(vLabel, "LinearOperators CG");
push!(vF, (vX, mP, mA, vV) -> KrylovTriCg!(vX, mP, mA, vV, ρ, σ; numFactor = numFactor, numIterations = numIterations, ϵSolver = ϵSolver, numItrSolver = numItrSolver));
push!(vLabel, "Kryloc TriCG");
push!(vF, (vX, mP, mA, vV) -> KrylovTriMr!(vX, mP, mA, vV, ρ, σ; numFactor = numFactor, numIterations = numIterations, ϵSolver = ϵSolver, numItrSolver = numItrSolver));
push!(vLabel, "Krylov TriMr");

numSolvers  = size(vF, 1);
cBenchMark  = Vector{Any}(undef, numSolvers);
vCorrSol    = falses(numSolvers);

tR = zeros(numDims, numSolvers, 5); #<! Min Time [Nano Sec], Median Time [Nano Sec], Max Time [Nano Sec], Number of Allocations, Allocations Size
for ii = 1:numDims
    mP, vQ, mA, vL, vU = GenerateRandomQP(problemClass, vNumElements[ii], numConstraints = vNumConstraints[ii]);

    # Some problems might change the actual data
    numElements     = size(mP, 1);
    numConstraints  = size(mA, 1);

    vX = randn(numElements);
    vV = randn(numElements + numConstraints);

    currTime    = time();
    vXRef       = (mP + (σ * sparse(I, numElements, numElements)) + (ρ * (mA' * mA))) \ vX;
    runTime     = time() - currTime;

    println("\nThe Reference Run Time: $runTime [Sec]\n");
    for jj = 1:numSolvers
        sBenchMark      = @benchmarkable vF[ii](vXX, $mP, $vQ, $mA, vV) setup = (vXX = copy($vX)); #<! We must interpolate the function as it is not in the global scope
        cBenchMark[ii]  = run(sBenchMark, samples = benchMarkSamples, evals = benchMarkEvals, seconds = benchMarkSeconds);
        vCorrSol[ii]    = norm(vXX - vXRef, 2) <= (numElements * ϵThr);

        mR[ii, 1] = min(cBenchMark[ii].times...);
        mR[ii, 2] = max(cBenchMark[ii].times...);
        mR[ii, 3] = median(cBenchMark[ii].times);
        mR[ii, 4] = cBenchMark[ii].allocs;
        mR[ii, 5] = cBenchMark[ii].memory;
    end
end

dataTitle = ["Solvers Min Run Time" "Solvers Max Run Time" "Solvers Median Run Time" "Solvers Allocations" "Solvers Allocations" "Solvers Success"]

for ii = 1:3
    display(scatter(1:numSolvers, mR[:, ii] * NS_TO_SEC_FCTR, title = dataTitle[ii], label = solversLabels, xlabel = "Solver", ylabel = "Time [Sec]"));
end

ii = 4;
display(scatter(1:numSolvers, mR[:, ii], title = dataTitle[ii], label = solversLabels, xlabel = "Solver", ylabel = "Number of Allocations"));

ii = 5;
display(scatter(1:numSolvers, mR[:, ii] * NS_TO_SEC_FCTR, title = dataTitle[ii], label = solversLabels, xlabel = "Solver", ylabel = "Size [MB]"));

ii = 6;
display(scatter(1:numSolvers, vCorrSol, title = dataTitle[ii], label = solversLabels, xlabel = "Solver", ylabel = "Solver Success Flag"));

# @benchmark LaLdlt!($vX, $mP, $mA, $vV; numFactor = 3, numIterations = 50, ϵSolver = 1e-6, numItrSolver = 1e-6, ρ = 1e6, σ = 1e-6);
# @benchmark ItrSolCg!($vX, $mP, $mA, $vV; numFactor = 3, numIterations = 50, ϵSolver = 1e-6, numItrSolver = 500, ρ = 1e6, σ = 1e-6);
# @benchmark KrylovCg!($vX, $mP, $mA, $vV; numFactor = 3, numIterations = 50, ϵSolver = 1e-6, numItrSolver = 500, ρ = 1e6, σ = 1e-6);
# @benchmark KrylovCr!($vX, $mP, $mA, $vV; numFactor = 3, numIterations = 50, ϵSolver = 1e-6, numItrSolver = 500, ρ = 1e6, σ = 1e-6);
# @benchmark KrylovCgLanczos!($vX, $mP, $mA, $vV; numFactor = 3, numIterations = 50, ϵSolver = 1e-6, numItrSolver = 500, ρ = 1e6, σ = 1e-6);
# @benchmark LinOpCg!($vX, $mP, $mA, $vV; numFactor = 3, numIterations = 50, ϵSolver = 1e-6, numItrSolver = 500, ρ = 1e6, σ = 1e-6);
# @benchmark KrylovTriCg!($vX, $mP, $mA, $vV; numFactor = 3, numIterations = 50, ϵSolver = 1e-6, numItrSolver = 500, ρ = 1e6, σ = 1e-6);
# @benchmark KrylovTriMr!($vX, $mP, $mA, $vV; numFactor = 3, numIterations = 50, ϵSolver = 1e-6, numItrSolver = 500, ρ = 1e6, σ = 1e-6);

# ρ = 1e6;
# σ = 1e-6;

# vT = zeros(numConstraints); #<! Buffer
# opL = LinearOperator(Float64, numElements, numElements, true, true, (vU, vW, α, β) -> begin
#     mul!(vT, mA, vW);
#     mul!(vU, transpose(mA), vT);
#     mul!(vU, mP, vW, one(Float64), ρ);
#     vU .= vU .+ (σ .* vW);
# end);

# mL = mP + (σ * sparse(I, numElements, numElements)) + (ρ * (transpose(mA) * mA));

# vX = rand(numElements);

# println(norm((opL * vX) - (mL * vX)))

# vT = mA * vX;
# vU = transpose(mA) * vT;
# vU .= mP * vX .+ ρ .* vU;
# vU .= vU .+ (σ * vX);

# println(norm(vU - mL * vX))

max_stars_repo_name: RoyiAvital/QuadraticProgramSolver

max_issues_repo_name: RoyiAvital/QuadraticProgramSolver

max_forks_repo_name: RoyiAvital/QuadraticProgramSolver


text:
function is_func_call(x::EXPR)
    if typof(x) === Call
        return true
    elseif iswherecall(x)
        return is_func_call(x.args[1])
    elseif isbracketed(x)
        return is_func_call(x.args[2])
    elseif isunarycall(x)
        return !(isoperator(x.args[1]) && (kindof(x.args[1]) === Tokens.EX_OR || kindof(x.args[1]) === Tokens.DECLARATION))
    elseif isbinarycall(x)
        if issyntaxcall(x.args[2])
            if is_decl(x.args[2])
                return is_func_call(x.args[1])
            else
                return false
            end
        else
            true
        end
    else
        return false
    end
end

is_assignment(x) = isbinarycall(x) && kindof(x.args[2]) === Tokens.EQ

# OPERATOR
is_exor(x) = isoperator(x) && kindof(x) === Tokens.EX_OR && x.dot == false
is_decl(x) = isoperator(x) && kindof(x) === Tokens.DECLARATION
is_issubt(x) = isoperator(x) && kindof(x) === Tokens.ISSUBTYPE
is_issupt(x) = isoperator(x) && kindof(x) === Tokens.ISSUPERTYPE
is_and(x) = isoperator(x) && kindof(x) === Tokens.AND && x.dot == false
is_not(x) = isoperator(x) && kindof(x) === Tokens.NOT && x.dot == false
is_plus(x) = isoperator(x) && kindof(x) === Tokens.PLUS && x.dot == false
is_minus(x) = isoperator(x) && kindof(x) === Tokens.MINUS && x.dot == false
is_star(x) = isoperator(x) && kindof(x) === Tokens.STAR && x.dot == false
is_eq(x) = isoperator(x) && kindof(x) === Tokens.EQ && x.dot == false
is_dot(x) = isoperator(x) && kindof(x) === Tokens.DOT
is_ddot(x) = isoperator(x) && kindof(x) === Tokens.DDOT
is_dddot(x) = isoperator(x) && kindof(x) === Tokens.DDDOT
is_pairarrow(x) = isoperator(x) && kindof(x) === Tokens.PAIR_ARROW && x.dot == false
is_in(x) = isoperator(x) && kindof(x) === Tokens.IN && x.dot == false
is_elof(x) = isoperator(x) && kindof(x) === Tokens.ELEMENT_OF && x.dot == false
is_colon(x) = isoperator(x) && kindof(x) === Tokens.COLON
is_prime(x) = isoperator(x) && kindof(x) === Tokens.PRIME
is_cond(x) = isoperator(x) && kindof(x) === Tokens.CONDITIONAL
is_where(x) = isoperator(x) && kindof(x) === Tokens.WHERE
is_anon_func(x) = isoperator(x) && kindof(x) === Tokens.ANON_FUNC

is_comma(x) = ispunctuation(x) && kindof(x) === Tokens.COMMA
is_lparen(x) = ispunctuation(x) && kindof(x) === Tokens.LPAREN
is_rparen(x) = ispunctuation(x) && kindof(x) === Tokens.RPAREN
is_lbrace(x) = ispunctuation(x) && kindof(x) === Tokens.LBRACE
is_rbrace(x) = ispunctuation(x) && kindof(x) === Tokens.RBRACE
is_lsquare(x) = ispunctuation(x) && kindof(x) === Tokens.LSQUARE
is_rsquare(x) = ispunctuation(x) && kindof(x) === Tokens.RSQUARE

# KEYWORD
is_if(x) = iskw(x) && kindof(x) === Tokens.IF
is_import(x) = iskw(x) && kindof(x) === Tokens.IMPORT

# Literals
is_lit_string(x) = kindof(x) === Tokens.STRING || kindof(x) === Tokens.TRIPLE_STRING

issubtypedecl(x::EXPR) = isbinarycall(x) && is_issubt(x.args[2])

rem_subtype(x::EXPR) = issubtypedecl(x) ? x[1] : x
rem_decl(x::EXPR) = isdeclaration(x) ? x[1] : x
rem_curly(x::EXPR) = typof(x) === Curly ? x.args[1] : x
rem_call(x::EXPR) = typof(x) === Call ? x[1] : x
rem_where(x::EXPR) = iswherecall(x) ? x[1] : x
rem_wheres(x::EXPR) = iswherecall(x) ? rem_wheres(x[1]) : x
rem_where_subtype(x::EXPR) = (iswherecall(x) || issubtypedecl(x)) ? x[1] : x
rem_where_decl(x::EXPR) = (iswherecall(x) || isdeclaration(x)) ? x[1] : x
rem_invis(x::EXPR) = isbracketed(x) ? rem_invis(x[2]) : x
rem_dddot(x::EXPR) = is_splat(x) ? x[1] : x
const rem_splat = rem_dddot
rem_kw(x::EXPR) = typof(x) === Kw ? x[1] : x

is_some_call(x) = typof(x) === Call || isunarycall(x) || (isbinarycall(x) && !(kindof(x.args[2]) === Tokens.DOT || issyntaxcall(x.args[2])))
is_eventually_some_call(x) = is_some_call(x) || ((isdeclaration(x) || iswherecall(x)) && is_eventually_some_call(x[1]))

defines_function(x::EXPR) = typof(x) === FunctionDef || (is_assignment(x) && is_eventually_some_call(x[1]))
defines_macro(x) = typof(x) == Macro
defines_datatype(x) = defines_struct(x) || defines_abstract(x) || defines_primitive(x)
defines_struct(x) = typof(x) === Struct || defines_mutable(x)
defines_mutable(x) = typof(x) === Mutable
defines_abstract(x) = typof(x) === Abstract
defines_primitive(x) = typof(x) === Primitive
defines_module(x) = typof(x) === ModuleH || typof(x) === BareModule
defines_anon_function(x) = isbinarycall(x) && is_anon_func(x.args[2])

has_sig(x::EXPR) = defines_datatype(x) || defines_function(x) || defines_macro(x) || defines_anon_function(x)

"""
    get_sig(x)

Returns the full signature of function, macro and datatype definitions.
Should only be called when has_sig(x) == true.
"""
function get_sig(x::EXPR)
    if isbinarycall(x)
        return x.args[1]
    elseif typof(x) === Struct || typof(x) === FunctionDef || typof(x) === Macro
        return x.args[2]
    elseif typof(x) === Mutable || typof(x) === Abstract || typof(x) === Primitive
        return x.args[3]
    end
end

function get_name(x::EXPR)
    if typof(x) === Struct || typof(x) === Mutable || typof(x) === Abstract || typof(x) === Primitive
        sig = get_sig(x)
        sig = rem_subtype(sig)
        sig = rem_wheres(sig)
        sig = rem_subtype(sig)
        sig = rem_curly(sig)
    elseif typof(x) === ModuleH || typof(x) === BareModule
        sig = x.args[2]
    elseif typof(x) === FunctionDef || typof(x) === Macro
        sig = get_sig(x)
        sig = rem_wheres(sig)
        sig = rem_decl(sig)
        sig = rem_call(sig)
        sig = rem_curly(sig)
        sig = rem_invis(sig)
        if isbinarycall(sig) && kindof(sig.args[2]) === Tokens.DOT
            if length(sig.args) > 2 && sig.args[3].args isa Vector{EXPR} && length(sig.args[3].args) > 0
                sig = sig.args[3].args[1]
            end
        end
        return sig
    elseif isbinarycall(x)
        length(x.args) < 2 && return x
        if kindof(x.args[2]) === Tokens.DOT
            if length(x.args) > 2 && typof(x.args[3]) === Quotenode && x.args[3].args isa Vector{EXPR} && length(x.args[3].args) > 0
                return get_name(x.args[3].args[1])
            else
                return x
            end
        end
        sig = x.args[1]
        if isunarycall(sig)
            return get_name(sig.args[1])
        end
        sig = rem_wheres(sig)
        sig = rem_decl(sig)
        sig = rem_call(sig)
        sig = rem_curly(sig)
        sig = rem_invis(sig)
        return get_name(sig)
    else
        sig = x
        if isunarycall(sig)
            sig = sig.args[1]
        end
        sig = rem_wheres(sig)
        sig = rem_decl(sig)
        sig = rem_call(sig)
        sig = rem_curly(sig)
        sig = rem_invis(sig)
    end
end

function get_arg_name(arg::EXPR)
    arg = rem_kw(arg)
    arg = rem_dddot(arg)
    arg = rem_where(arg)
    arg = rem_decl(arg)
    arg = rem_subtype(arg)
    arg = rem_curly(arg)
    arg = rem_invis(arg)
end


text:
##### Interfaces #####

# We must convert ydot into an array and assign element-by-element
# This is because we are modifying by-reference an NVector in Sundials...
function interface_ode(t::Sundials.realtype, y::Sundials.N_Vector, ydot::Sundials.N_Vector, user_data::Vector{Any})
  y = Sundials.asarray(y)
  ydot = Sundials.asarray(ydot)
  user_data[1](t, y, user_data[2], user_data[3], ydot)
  return int32(0)
end

##### Simulation functions #####

# Function to simulate the model
@inbounds function simulate(times::Vector{Float64},
                            states::Vector{Float64},
                            parameters::Vector{Float64},
                            inputs::Vector{Function},
                            settings::Dict{ASCIIString,Any},
                            model::Function,
                            observer::Function)
  neq = length(states)

  if settings["method"] == "bdf"
    cvode_mem = Sundials.CVodeCreate(Sundials.CV_BDF, Sundials.CV_NEWTON)
  elseif settings["method"] == "adams"
    cvode_mem = Sundials.CVodeCreate(Sundials.CV_ADAMS, Sundials.CV_NEWTON)
  else
    error("Please choose bdf or adams as method")
  end
# Shut up Sundials
  flag = Sundials.CVodeSetErrFile(cvode_mem, C_NULL )
  if flag < zero(flag)
    error("Error in the CVodeSetErrFile function")
  end

  # Initialize the Sundials solver by passing the wrapper function
  flag = Sundials.CVodeInit(cvode_mem, cfunction(interface_ode, Int32,
    (Sundials.realtype, Sundials.N_Vector, Sundials.N_Vector,  Vector{Any})),
             times[1], Sundials.nvector(states))
  if flag < zero(flag)
    error("Error in the CVodeInit function")
  end

  # Tell Sundials the tolerance settings for error control
  flag = Sundials.CVodeSStolerances(cvode_mem, settings["rtol"], settings["atol"])
  if flag < zero(flag)
    error("Error in the CVodeSStolerances function")
  end

  # Tell Sundials the number of state variables, so that I can allocate memory for the Jacobian
  flag = Sundials.CVDense(cvode_mem, neq)
  if flag < zero(flag)
    error("Error in the CVDense function")
  end

  # Pass information for the callback
  flag = Sundials.CVodeSetUserData(cvode_mem,{model, parameters, inputs, observed})
  if flag < zero(flag)
    error("Error in the CVodeSetUserData function")
  end

  # Set maximum number of steps
  flag = Sundials.CVodeSetMaxNumSteps(cvode_mem, settings["maxsteps"])
  if flag < zero(flag)
    error("Error in the CVodeSetMaxNumSteps function")
  end

  # Set maximum order of the integration
  flag = Sundials.CVodeSetMaxOrd(cvode_mem, settings["maxord"])
  if flag < zero(flag)
    error("Error in the CVodeSetMaxOrd function")
  end

  # Set the initial step size
  flag = Sundials.CVodeSetInitStep(cvode_mem, settings["hini"])
  if flag < zero(flag)
    error("Error in the CVodeSetInitStep function")
  end

  # Set the minimum step size
  flag = Sundials.CVodeSetMinStep(cvode_mem, settings["hmin"])
  if flag < zero(flag)
    error("Error in the CVodeSetMinStep function")
  end

  # Set the maximum step size
  flag = Sundials.CVodeSetMaxStep(cvode_mem, settings["hmax"])
  if flag < zero(flag)
    error("Error in the CVodeSetMaxStep function")
  end

  # Set the maximum number of error test fails
  flag = Sundials.CVodeSetMaxErrTestFails(cvode_mem, settings["maxerr"])
  if flag < zero(flag)
    error("Error in the CVodeSetMaxErrTestFails function")
  end

  # Set the maximum number of nonlinear iterations per step
  flag = Sundials.CVodeSetMaxNonlinIters(cvode_mem, settings["maxnonlin"])
  if flag < zero(flag)
    error("Error in the CVodeSetMaxNonlinIters function")
  end

  # Set the maximum number of convergence failures
  flag = Sundials.CVodeSetMaxConvFails(cvode_mem, settings["maxconvfail"])
  if flag < zero(flag)
    error("Error in the CVodeSetMaxConvFails function")
  end

  # Make a first call to the model to check that everything is ok and retrieve the number of observed variables
  observed = zeros(Float64, settings["nobs"])
  observer(times[1], states, parameters, inputs, observed)

  # Use time as column -> Because of column-ordered
  output = zeros(Float64, (neq + settings["nobs"] + 1, length(times)))
  output[:,1] = [times[1], states, observed]

  # Main time loop. Each timestep call cvode. Handle exceptions and fill up output
  t = times[1]
  y = copy(states)
  for i in 2:length(times)
    # Note that the second t must be an array to build the pointer...
    flag = Sundials.CVode(cvode_mem, times[i], y, [t], Sundials.CV_NORMAL)
    if flag < zero(flag)
       if flag == oftype(flag, -1)
          error("The solver took mxstep internal steps but could not reach tout.")
       elseif flag == oftype(flag, -2)
          error("The solver could not satisfy the accuracy demanded by the user for some internal step")
       elseif flag == oftype(flag, -3)
          error("Error test failures occured too many times during one internal time step or minimum step size was reached")
       elseif flag == oftype(flag, -4)
          error("Convergence test failures occurred too many times during one internal time step or minimum step size was reached.")
       elseif flag == oftype(flag, -5)
          error("The linear solver’s initialization function failed.")
       elseif flag == oftype(flag, -6)
          error("The linear solver’s setup function failed in an unrecoverable manner")
       elseif flag == oftype(flag, -7)
          error("The linear solver’s solve function failed in an unrecoverable manner")
       elseif flag == oftype(flag, -8)
          error("The right hand side function failed in an unrecoverable manner")
       elseif flag == oftype(flag, -9)
          error("The right-hand side function failed at the first call.")
       elseif flag == oftype(flag, -10)
          error("The right-hand side function had repeated recoverable errors.")
       elseif flag == oftype(flag, -11)
          error("The right-hand side function had a recoverable errors but no recovery is possible.")
       elseif flag == oftype(flag, -25)
          error("The time t is outside the last step taken.")
       elseif flag == oftype(flag, -26)
          error("The output derivative vector is NULL.")
       elseif flag == oftype(flag, -27)
          error("The output and initial times are too close to each other.")
       end
     end
    output[1,i] = times[i]
    for h in 1:neq
      output[h+1,i] = y[h]
    end
  end
  # If we have observed variables we call the model function again
  if settings["nobs"] > 0
    for i in 1:length(times)
      observer(times[i], vec(output[2:(neq + 1), i]), parameters, inputs, observed)
      for h in 1:settings["nobs"]
        output[neq + 1 + h, i] = observed[h]
      end
    end
  end

  return transpose(output)
end

## Return the solution of the model at the specific timepoints
function simulate(model::OdeModel, time::Vector{Float64}, settings::Dict{ASCIIString, Any})
    simulation = simulate(time, convert(Array, model.States), convert(Array, model.Parameters),
                            convert(Array, model.Inputs), settings, model.Model, model.Observer)
end

# Generate Dict with default settings.

function CvodeSettings()
  settings = Dict{ASCIIString, Any}()
  settings["method"] = "bdf"
  settings["rtol"] = 1e-6
  settings["atol"] = 1e-6
  settings["maxsteps"] = 1e5
  settings["maxord"] = 5
  settings["hini"] = 0
  settings["hmin"] = 0
  settings["hmax"] = 0
  settings["maxerr"] = 12
  settings["maxnonlin"] = 12
  settings["maxconvfail"] = 12
  settings["jacobian"] = false
  return settings
end


text:
function optimize(d::DifferentiableFunction,
                  initial_x::Vector;
                  method::Symbol = :l_bfgs,
                  tolerance::Real = 1e-8,
                  iterations::Integer = 1_000,
                  store_trace::Bool = false,
                  show_trace::Bool = false)
    if method == :gradient_descent
        gradient_descent(d,
                         initial_x,
                         tolerance = tolerance,
                         iterations = iterations,
                         store_trace = store_trace,
                         show_trace = show_trace)
    elseif method == :bfgs
        bfgs(d,
             initial_x,
             tolerance = tolerance,
             iterations = iterations,
             store_trace = store_trace,
             show_trace = show_trace)
    elseif method == :l_bfgs
        l_bfgs(d,
               initial_x,
               tolerance = tolerance,
               iterations = iterations,
               store_trace = store_trace,
               show_trace = show_trace)
    else
        throw(ArgumentError("Unknown method $method"))
    end
end

function optimize(f::Function,
                  g!::Function,
                  h!::Function,
                  initial_x::Vector;
                  method::Symbol = :newton,
                  tolerance::Real = 1e-8,
                  iterations::Integer = 1_000,
                  store_trace::Bool = false,
                  show_trace::Bool = false)
    if method == :nelder_mead
        nelder_mead(f,
                    initial_x,
                    tolerance = tolerance,
                    iterations = iterations,
                    store_trace = store_trace,
                    show_trace = show_trace)
    elseif method == :simulated_annealing
        simulated_annealing(f,
                            initial_x,
                            tolerance = tolerance,
                            iterations = iterations,
                            store_trace = store_trace,
                            show_trace = show_trace)
    elseif method == :gradient_descent
        d = DifferentiableFunction(f, g!)
        gradient_descent(d,
                         initial_x,
                         tolerance = tolerance,
                         iterations = iterations,
                         store_trace = store_trace,
                         show_trace = show_trace)
    elseif method == :newton
        d = TwiceDifferentiableFunction(f, g!, h!)
        newton(d,
               initial_x,
               tolerance = tolerance,
               iterations = iterations,
               store_trace = store_trace,
               show_trace = show_trace)
    elseif method == :bfgs
        d = DifferentiableFunction(f, g!)
        bfgs(d,
             initial_x,
             tolerance = tolerance,
             iterations = iterations,
             store_trace = store_trace,
             show_trace = show_trace)
    elseif method == :l_bfgs
        d = DifferentiableFunction(f, g!)
        l_bfgs(d,
               initial_x,
               tolerance = tolerance,
               iterations = iterations,
               store_trace = store_trace,
               show_trace = show_trace)
    else
        throw(ArgumentError("Unknown method $method"))
    end
end

function optimize(f::Function,
                  g!::Function,
                  initial_x::Vector;
                  method::Symbol = :l_bfgs,
                  tolerance::Real = 1e-8,
                  iterations::Integer = 1_000,
                  store_trace::Bool = false,
                  show_trace::Bool = false)
    if method == :nelder_mead
        nelder_mead(f,
                    initial_x,
                    tolerance = tolerance,
                    iterations = iterations,
                    store_trace = store_trace,
                    show_trace = show_trace)
    elseif method == :simulated_annealing
        simulated_annealing(f,
                            initial_x,
                            tolerance = tolerance,
                            iterations = iterations,
                            store_trace = store_trace,
                            show_trace = show_trace)
    elseif method == :gradient_descent
        d = DifferentiableFunction(f, g!)
        gradient_descent(d,
                         initial_x,
                         tolerance = tolerance,
                         iterations = iterations,
                         store_trace = store_trace,
                         show_trace = show_trace)
    elseif method == :bfgs
        d = DifferentiableFunction(f, g!)
        bfgs(d,
             initial_x,
             tolerance = tolerance,
             iterations = iterations,
             store_trace = store_trace,
             show_trace = show_trace)
    elseif method == :l_bfgs
        d = DifferentiableFunction(f, g!)
        l_bfgs(d,
               initial_x,
               tolerance = tolerance,
               iterations = iterations,
               store_trace = store_trace,
               show_trace = show_trace)
    else
        throw(ArgumentError("Unknown method $method"))
    end
end

function optimize(f::Function,
                  initial_x::Vector;
                  method::Symbol = :nelder_mead,
                  tolerance::Real = 1e-8,
                  iterations::Integer = 1_000,
                  store_trace::Bool = false,
                  show_trace::Bool = false)
    if method == :nelder_mead
        nelder_mead(f,
                    initial_x,
                    tolerance = tolerance,
                    iterations = iterations,
                    store_trace = store_trace,
                    show_trace = show_trace)
    elseif method == :simulated_annealing
        simulated_annealing(f,
                            initial_x,
                            tolerance = tolerance,
                            iterations = iterations,
                            store_trace = store_trace,
                            show_trace = show_trace)
    elseif method == :gradient_descent
        d = DifferentiableFunction(f)
        gradient_descent(d,
                         initial_x,
                         tolerance = tolerance,
                         iterations = iterations,
                         store_trace = store_trace,
                         show_trace = show_trace)
    elseif method == :bfgs
        d = DifferentiableFunction(f)
        bfgs(d,
             initial_x,
             tolerance = tolerance,
             iterations = iterations,
             store_trace = store_trace,
             show_trace = show_trace)
    elseif method == :l_bfgs
        d = DifferentiableFunction(f)
        l_bfgs(d,
               initial_x,
               tolerance = tolerance,
               iterations = iterations,
               store_trace = store_trace,
               show_trace = show_trace)
    else
        throw(ArgumentError("Unknown method $method"))
    end
end


text:
function cooler!(Hs::AdiaComputer,gamma::Real,t::Real)
    # boundscheck(Hs,gamma,t)
    energy = eig(full(Hamiltonian(Hs)))[1]
    @assert -pi/2<minimum(energy)*t-gamma<pi/2 "bad cooling parameters"
    @assert -pi/2<maximum(energy)*t-gamma<pi/2 "bad cooling parameters"

    Hs.state = normalize!((Hs.state-im*exp(im*gamma)*
            trotter(-im*t*(1-Hs.location)*Hs.HB,
                -im*t*Hs.location*Hs.HP,3)*Hs.state)/2)
end

function heater!(Hs::AdiaComputer,gamma::Real,t::Real)
    # boundscheck(Hs,gamma,t)
    energy = eig(full(Hamiltonian(Hs)))[1]
    @assert -pi/2<minimum(energy)*t-gamma<pi/2 "bad cooling parameters"
    @assert -pi/2<maximum(energy)*t-gamma<pi/2 "bad cooling parameters"

    Hs.state = normalize!((Hs.state+im*exp(im*gamma)*
                trotter(-im*t*(1-Hs.location)*Hs.HB,
                    -im*t*Hs.location*Hs.HP,3)*Hs.state)/2)
end

function daemon!(
    Hs::AdiaComputer,
    gamma::Real,
    t::Real
    )
    dice = rand()
    if dice <= 0.5*(1+sin(gamma))
        #get |1> (higher energy)
        heater!(Hs,gamma,t)
        return 0.5*(1-sin(gamma))
    else
        #get |0> (lower energy)
        cooler!(Hs,gamma,t)
        return 0.5*(1+sin(gamma))
    end
end

function cooling!(
    Hs::AdiaComputer;
    n=5)

    count = 0
    gamma,t = CoolingPara(Hs)

    while count<n
        Hs.prob *= daemon!(Hs,gamma,t)
        count += 1
    end
    return Hs
end

function CoolingPara(Hs::AdiaComputer)
    Eigen = eig(full(Hamiltonian(Hs)))[1]
    maxEigen = maximum(Eigen)
    minEigen = minimum(Eigen)

    gamma = (maxEigen+minEigen)/(maxEigen-minEigen) * pi/2 * 0.1

    if (gamma-pi/2)>0
        t = 0.5*((gamma-pi/2)/minEigen + (gamma+pi/2)/maxEigen)
    else
        t = 0.5*((gamma-pi/2)/maxEigen + (gamma+pi/2)/maxEigen)
    end
    return gamma,t
end

export cooling!

max_stars_repo_name: Roger-luo/QuCmp

max_issues_repo_name: Roger-luo/QuCmp

max_forks_repo_name: Roger-luo/QuCmp


text:
@testset "benchmark_image" begin
    img = generate_surface(SynthSphere(50), 1.0, [0,0,1])
    img2 = generate_surface(SynthSphere(51), 1.0, [0,0,1])
    error = benchmark_image(img, img2)
    @test error ≈ 0.10657885454612316
end


text:
using Polytopes
using Colors
using GraphPlot
using LightGraphs

p = Polytope([[1, 2], [2, 3], [3, 1], [2, 3], [3, 4], [4, 2]], [[1, 2, 3], [4, 5, 6]])
g = p.graph
is = Polytopes.start_indices(p.connectivity)
colors = [colorant"red", colorant"turquoise", colorant"green"]
gplot(g, nodelabel=vertices(g), nodefillc=[colors[something(findfirst(x -> i ≤ x, is), 0) + 1] for i ∈ vertices(g)])


text:
export iauGst06
"""
Greenwich apparent sidereal time, IAU 2006, given the NPB matrix.

This function is part of the International Astronomical Union's
SOFA (Standards Of Fundamental Astronomy) software collection.

Status:  support function.

Given:
   uta,utb  double        UT1 as a 2-part Julian Date (Notes 1,2)
   tta,ttb  double        TT as a 2-part Julian Date (Notes 1,2)
   rnpb     double[3][3]  nutation x precession x bias matrix

Returned (function value):
            double        Greenwich apparent sidereal time (radians)

Notes:

   1. The UT1 and TT dates uta+utb and tta+ttb respectively, are both
   Julian Dates, apportioned in any convenient way between the
   argument pairs.  For example, JD=2450123.7 could be expressed in
   any of these ways, among others:

         Part A        Part B

      2450123.7           0.0       (JD method)
      2451545.0       -1421.3       (J2000 method)
      2400000.5       50123.2       (MJD method)
      2450123.5           0.2       (date & time method)

   The JD method is the most natural and convenient to use in
   cases where the loss of several decimal digits of resolution
   is acceptable (in the case of UT;  the TT is not at all critical
   in this respect).  The J2000 and MJD methods are good compromises
   between resolution and convenience.  For UT, the date & time
   method is best matched to the algorithm that is used by the Earth
   rotation angle function, called internally:  maximum precision is
   delivered when the uta argument is for 0hrs UT1 on the day in
   question and the utb argument lies in the range 0 to 1, or vice
   versa.

   2. Both UT1 and TT are required, UT1 to predict the Earth rotation
   and TT to predict the effects of precession-nutation.  If UT1 is
   used for both purposes, errors of order 100 microarcseconds
   result.

   3. Although the function uses the IAU 2006 series for s+XY/2, it is
   otherwise independent of the precession-nutation model and can in
   practice be used with any equinox-based NPB matrix.

   4. The result is returned in the range 0 to 2pi.

Called:
   iauBpn2xy    extract CIP X,Y coordinates from NPB matrix
   iauS06       the CIO locator s, given X,Y, IAU 2006
   iauAnp       normalize angle into range 0 to 2pi
   iauEra00     Earth rotation angle, IAU 2000
   iauEors      equation of the origins, given NPB matrix and s

Reference:

   Wallace, P.T. & Capitaine, N., 2006, Astron.Astrophys. 459, 981

This revision:  2013 June 18

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# double iauGst06(double uta, double utb, double tta, double ttb,
#                 double rnpb[3][3])

function iauGst06(uta::Real, utb::Real, tta::Real, ttb::Real,
                   rnpb::AbstractMatrix{<:Real})

   return ccall((:iauGst06, libsofa_c), Cdouble,
               (Cdouble, Cdouble, Cdouble, Cdouble, Ptr{Cdouble}),
               convert(Float64, uta), convert(Float64, utb),
               convert(Float64, tta), convert(Float64, ttb),
               convert(Array{Float64, 2}, rnpb'))
end


text:
module IRImodels
    using CategoricalArrays, DataArrays, DataFrames, DataStructures, Distributions,
        GLM, JSON, MixedModels, NamedArrays
    using StatsBase: CoefTable, counts
    export
        checksingularity,
        dropterm,
        drop1,
        GLMformula,
        GLMMform,
        profileθ,
        relevel,
        vif

    include("AoD.jl")

end # module

max_stars_repo_name: werbac/jlib

max_issues_repo_name: werbac/jlib

max_forks_repo_name: werbac/jlib


text:
# u ~ x ~ j
# v ~ y ~ i

function source_point_mat(X, Y)
    x1, x2, x3 = X
    y1, y2, y3 = Y
    F = [x1 x2 x3; y1 y2 y3; 1 1 1]
    F
end

function target_point_mat(U, V)
    u1, u2, u3 = U
    v1, v2, v3 = V
    T = [u1 u2 u3; v1 v2 v3]
    T
end

function affine_params(X, Y, U, V)
    F = source_point_mat(X, Y)
    T = target_point_mat(U, V)
    M = T * inv(F)
    M
end

warp_pixel(M, x::Float64, y::Float64) = M * [x, y, 1]

function pa_warp{N}(img::Array{Float64, N},
                    src::Matrix{Float64}, dst::Matrix{Float64},
                    trigs::Matrix{Int})
    warped = zeros(eltype(img), size(img))
    for t=1:size(trigs, 1)
        tr = slice(trigs, t, :)
        Y = src[tr, 1]
        X = src[tr, 2]

        V = dst[tr, 1]
        U = dst[tr, 2]

        # warp parameters from target (U, V) to source (X, Y)
        M = affine_params(U, V, X, Y)

        mask = poly2mask(U, V, size(img, 1), size(img, 2))
        vs, us = findn(mask)

        # for every pixel in target triangle we find corresponding pixel in source
        # and copy its value
        for i=1:length(vs)
            u, v = us[i], vs[i]
            x, y = warp_pixel(M, Float64(u), Float64(v))
            if 1 <= y && y <= size(img, 1) && 1 <= x && x <= size(img, 2)
                warped[v, u, :] = img[round(Int, y), round(Int, x), :]
            end
        end

    end
    warped
end


text:
function getNbors(M::Matrix{Int}, pos::CartesianIndex{2}) :: Vector{CartesianIndex{2}}
  s = size(M)
  nbors :: Vector{CartesianIndex{2}} = []
  (pos[1] != 1   ) && push!(nbors, pos + CartesianIndex(-1,0))
  (pos[1] != s[1]) && push!(nbors, pos + CartesianIndex(1,0))
  (pos[2] != 1   ) && push!(nbors, pos + CartesianIndex(0,-1))
  (pos[2] != s[2]) && push!(nbors, pos + CartesianIndex(0,1))
  (pos[1] != 1   ) && (pos[2] != 1   ) && push!(nbors, pos + CartesianIndex(-1,-1))
  (pos[1] != 1   ) && (pos[2] != s[2]) && push!(nbors, pos + CartesianIndex(-1, 1))
  (pos[1] != s[1]) && (pos[2] != 1   ) && push!(nbors, pos + CartesianIndex( 1,-1))
  (pos[1] != s[1]) && (pos[2] != s[2]) && push!(nbors, pos + CartesianIndex( 1, 1))
  return nbors
end

# Performs a flash cycle, returns total number of flashes.
function step!(M::Matrix{Int}) :: Int
  flashed :: Matrix{Bool} = zeros(size(M))

  # A) Everyone increases energy by 1
  M .+= 1

  # B) Octopuses flash
  charged = findall((M .> 9) .& .!flashed)
  while !isempty(charged)
    # Flash
    map(pos -> M[getNbors(M, pos)] .+= 1, charged)
    flashed[charged] .= true

    # Reset
    charged = findall((M .> 9) .& .!flashed)
  end

  # C) All octopuses that have flashed go back to 0
  M[flashed] .= 0
  return count(flashed)
end

function readInput() :: Matrix{Int}
  f = open("./inputs/day11.txt")
  data :: Matrix{Int} = zeros(10,10)

  i = 1
  for line in eachline(f)
    data[i,:] = map(x -> parse(Int,x), collect(line))
    i += 1
  end
  return data
end

problem1!(M::Matrix{Int}) = sum(map(i -> step!(M), 1:100))

function problem2!(M::Matrix{Int}) :: Int
  i :: Int = 1
  while step!(M) != 100
    i += 1
  end
  return i
end

println("Problem 1:")
M = readInput()
println(problem1!(M))

println("Problem 2:")
M = readInput()
println(problem2!(M))


text:
using SpecialFunctions
#Copyright (c) 2010, Xin-She Yang
#All rights reserved.
# Trasnlation to Julia language, Fredrik bagge Carlson
#
#Redistribution and use in source and binary forms, with or without
#modification, are permitted provided that the following conditions are met:
#
#* Redistributions of source code must retain the above copyright notice, this
#  list of conditions and the following disclaimer.
#
#* Redistributions in binary form must reproduce the above copyright notice,
#  this list of conditions and the following disclaimer in the documentation
#  and/or other materials provided with the distribution
#THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
#AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
#IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
#DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
#FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
#DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
#SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
#CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
#OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
#OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

# """
#     cuckoo_search(f,X0;
# `lb = -Inf`,
# `ub = Inf`,
# `n = Number of nests (or different solutions)
# `pa = 0.25` Discovery rate of alien eggs/solutions
# `tol = 1.0e-5`,
# `max_iter = 1e5`,
# `timeout = Inf`)
#
# Change this if you want to get better results
# Based on implementation by
# @inproceedings{yang2009cuckoo,
#   title={Cuckoo search via L{\'e}vy flights},
#   author={Yang, Xin-She and Deb, Suash},
#   booktitle={Nature \& Biologically Inspired Computing, 2009. NaBIC 2009. World Congress on},
#   pages={210--214},
#   year={2009},
#   organization={IEEE}
# }
#
# """
function cuckoo_search(f, X0;
    lb=-Inf,
    ub = Inf,
    n = 25,
    pa = 0.25,
    tol = 1.0e-5,
    max_iter = 1e5,
    timeout = Inf)

    nd = length(X0) # n dims

    if lb == -Inf
        lb = X0-0.99999*abs.(X0)
    end
    if ub == Inf
        ub = X0+0.99999*abs.(X0)
    end

    # Random initial solutions
    nest = [copy(X0) for _ = 1:n] # n = num_nests
    nest[1] = X0
    for i=2:n
        nest[i] = @. lb + (ub-lb)*rand()
    end
    # Get the current best
    fitness = fill(Inf,n)
    fmin, bestnest, nest, fitness = update_nests!(f,nest,nest,fitness)
    new_nest = deepcopy(nest)
    N_iter=0
    ## Starting iterations
    t0 = time()
    while fmin > tol && N_iter < max_iter
        # Generate new solutions (but keep the current best)
        move_cuckoos!(nest,new_nest,bestnest,lb,ub)
        fnew, best, nest, fitness = update_nests!(f,nest,new_nest,fitness)
        # Update the counter
        N_iter += n
        if fnew < fmin
            fmin = fnew
            bestnest = best
        end
        # Discovery and randomization
        new_nest = empty_nests(nest,lb,ub,pa)
        # Evaluate this set of solutions
        fnew, best, nest, fitness = update_nests!(f,nest,new_nest,fitness)
        # Update the counter again
        N_iter += n
        # Find the best objective so far
        if fnew < fmin
            fmin = fnew
            bestnest = best
        end
        if time()-t0 > timeout
            @info("Cuckoo search: timeout $(timeout)s reached ($(time()-t0)s)")
            break
        end
    end ## End of iterations
    ## Post-optimization processing
    ## Display all the nests
    bestnest, fmin

end
## --------------- All subfunctions are list below ------------------
## Get cuckoos by ramdom walk
function move_cuckoos!(nest,new_nest,best,lb,ub)
    # Levy flights
    n = length(nest)
    # Levy exponent and coefficient
    # For details, see equation (2.21), Page 16 (chapter 2) of the book
    # X. S. Yang, Nature-Inspired Metaheuristic Algorithms, 2nd Edition, Luniver Press, (2010).
    β = 3/2
    βi = 1/β
    σ = (gamma(1+β)*sin(pi*β/2)/(gamma((1+β)/2)*β*2^((β-1)/2)))^(1/β)
    u,v,step = similar(best),similar(best),similar(best)
    for j = 1:n
        s = nest[j]
        # This is a simple way of implementing Levy flights
        # For standard random walks, use step=1;
        ## Levy flights by Mantegna’s algorithm
        u .= σ.*randn.()
        v .= randn.()
        @. step = u/abs(v)^βi
        # In the next equation, the difference factor (s-best) means that
        # when the solution is the best solution, it remains unchanged.
        @. step = 0.01*step.*(s-best)
        # Here the factor 0.01 comes from the fact that L/100 should the typical
        # step size of walks/flights where L is the typical lenghtscale;
        # otherwise, Levy flights may become too aggresive/efficient,
        # which makes new solutions (even) jump out side of the design domain
        # (and thus wasting evaluations).
        # Now the actual random walks or flights
        s += step .* randn.()
        # Apply simple bounds/limits
        new_nest[j] = clamp(s,lb,ub)
    end
end
## Find the current best nest
function update_nests!(f,nest,new_nest,fitness)
    # Evaluating all new solutions
    for j = 1:length(nest)
        fnew = f(new_nest[j])
        if fnew <= fitness[j]
            fitness[j] = fnew
            nest[j] = new_nest[j]
        end
    end
    # Find the current best
    fmin,K = findmin(fitness)
    best = nest[K]
    fmin,best,nest,fitness
end

## Replace some nests by constructing new solutions/nests
function empty_nests(nest,lb,ub,pa)
    # A fraction of worse nests are discovered with a probability pa
    n = length(nest)
    # Discovered or not -- a status vector
    K = rand(size(nest)) .> pa
    # In the real world, if a cuckoo’s egg is very similar to a host’s eggs, then
    # this cuckoo’s egg is less likely to be discovered, thus the fitness should
    # be related to the difference in solutions. Therefore, it is a good idea
    # to do a random walk in a biased way with some random step sizes.
    ## New solution by biased/selective random walks
    stepsize = rand().*(nest[randperm(n)] .- nest[randperm(n)])
    new_nest = nest .+ stepsize .* K
    for j = 1:length(nest)
        new_nest[j] = clamp(new_nest[j],lb,ub)
    end
    return new_nest
end

function Base.clamp(x::AbstractVector,l::AbstractVector,u::AbstractVector)
    map(x,l,u) do x,l,u
        clamp(x,l,u)
    end
end

using BenchmarkTools, StaticArrays
# 78.619 ms (861606 allocations: 61.70 MiB)
# 6.321 ms (146181 allocations: 12.43 MiB)
rosenbrock(x) =  (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
cuckoo_search(rosenbrock, [10.,10], tol=1e-9)
@btime cuckoo_search(rosenbrock, [10.,10.], lb=[-20,-20], ub=[20,20], tol=1e-9)
@btime cuckoo_search(rosenbrock, @SVector([10.,10.]), lb=@SVector([-20,-20]), ub=@SVector([20,20]), tol=1e-9)

@profiler cuckoo_search(rosenbrock, @SVector([10.,10.]), lb=@SVector([-20,-20]), ub=@SVector([20,20]), tol=1e-9)


text:
module TopRequest

import ..Layer, ..request

export TopLayer

"""
    request(TopLayer, args...; kwargs...)

This layer is at the top of every stack, and does nothing.
It's useful for inserting a custom layer at the top of the stack.
"""
abstract type TopLayer{Next <: Layer} <: Layer{Next} end

request(::Type{TopLayer{Next}}, args...; kwargs...) where Next =
    request(Next, args...; kwargs...)

end


text:
@testset "FillDictionary" begin
    inds = Indices([2, 4, 6])
    dict = FillDictionary(inds, 100.0)

    @test keytype(dict) == Int
    @test eltype(dict) == Float64
    @test keys(dict) === inds
    @test dict[2] === 100.0
    @test_throws IndexError dict[3]
    @test dict == Dictionary([2, 4, 6], [100.0, 100.0, 100.0])
    @test dict == FillDictionary([2, 4, 6], 100.0)
    @test istokenizable(dict) == true
end


text:
#!/usr/bin/env julia

# This file is used to create the plots in Figure 2 of the main paper
using Plots
#reset past font scaling
Plots.scalefontsizes()
#scale fonts
Plots.scalefontsizes(2.0)

include("./optimizers_abstract.jl")
ε = 1e-6
γ = 1.0 #1.00
η = 0.20 #0.20
xmin = -2.0*2
xmax = 2.0*2
ymin = -2.0*2
ymax = 2.0*2

nIter = 50 #50
xHistGDA = zeros(1, nIter + 1)
yHistGDA = zeros(1, nIter + 1)
xHistLCGD = zeros(1, nIter + 1)
yHistLCGD = zeros(1, nIter + 1)
xHistConOpt = zeros(1, nIter + 1)
yHistConOpt = zeros(1, nIter + 1)
xHistOGDA = zeros(1, nIter + 1)
yHistOGDA = zeros(1, nIter + 1)
xHistCGDA = zeros(1, nIter + 1)
yHistCGDA = zeros(1, nIter + 1)

# α  in 1, 3, 6
α = 3
f(x,y) = α * x'*y
g(x,y) = - α * x'*y

x0 = [1.5]
y0 = [2.5]

xHistGDA[:,1] .= x0
yHistGDA[:,1] .= y0
xHistLCGD[:,1] .= x0
yHistLCGD[:,1] .= y0
xHistConOpt[:,1] .= x0
yHistConOpt[:,1] .= y0
xHistOGDA[:,1] .= x0
yHistOGDA[:,1] .= y0
xHistCGDA[:,1] .= x0
yHistCGDA[:,1] .= y0

∇x_f, ∇y_f, ∇x∇x_f, ∇x∇y_f, ∇y∇y_f, ∇y∇x_f = get_closures(f)
∇x_g, ∇y_g, ∇x∇x_g, ∇x∇y_g, ∇y∇y_g, ∇y∇x_g = get_closures(g)

ΔxC = zeros(1)
ΔyC = zeros(1)

ΔxO = zeros(1)
ΔyO = zeros(1)
for k = 1 : nIter
  global ΔxC
  global ΔyC
  global ΔxO
  global ΔyO
  x, y = xHistGDA[:, k], yHistGDA[:, k]
  GDA!(x,
       y,
       ∇x_f, ∇y_g, η)
  xHistGDA[:, k + 1], yHistGDA[:, k + 1] = x,y

  x, y = xHistOGDA[:, k], yHistOGDA[:, k]
  OGDA!(x,
        y,
        ΔxO, ΔyO, ∇x_f, ∇y_g, η)
  xHistOGDA[:, k + 1], yHistOGDA[:, k + 1] = x, y

  x, y = xHistCGDA[:, k], yHistCGDA[:, k]
  CG_CGDA!(x,
           y,
           ΔxC, ΔyC , ∇x_f, ∇x∇y_f, ∇y_g, ∇y∇x_g, η, ε)
  xHistCGDA[:, k + 1], yHistCGDA[:, k + 1] = x, y

  x, y = xHistLCGD[:, k], yHistLCGD[:, k]
  LCGDA!(x,
         y,
         ∇x_f, ∇x∇y_f, ∇y_g, ∇y∇x_g, η,)
  xHistLCGD[:, k + 1], yHistLCGD[:, k + 1] = x, y

  x, y = xHistConOpt[:, k], yHistConOpt[:, k]
  conOpt!(x,
          y,
          ∇x_f, ∇y_g, ∇x∇x_f, ∇x∇y_f, ∇y∇y_g, ∇y∇x_g, η, γ)
  xHistConOpt[:, k + 1], yHistConOpt[:, k + 1] = x, y
end

pl = plot(xlim=[xmin,xmax], ylim=[ymin,ymax], xlabel="\$x\$", ylabel="\$y\$",
          legend=:topleft)

# plot!(pl, xHistGDA[1:1:end], yHistGDA[1:1:end],
#       marker=:cross, markersize=5.0, label="GDA", linestyle=:dot,
#       linewidth=5.0)
plot!(pl, xHistLCGD[1:1:end], yHistLCGD[1:1:end],
      marker=:ltriangle, markersize=5.0, label="LCGD", linestyle=:dash,
      linewidth=5.0)
# plot!(pl, xHistConOpt[1:1:end], yHistConOpt[1:1:end],
#       marker=:star, markersize=5.0, label="SGA", linestyle=:dashdotdot,
#       linewidth=5.0)
# plot!(pl, xHistOGDA[1:1:end], yHistOGDA[1:1:end],
#       marker=:xcross, markersize=5.0, label="OGDA", linestyle=:dashdot,
#       linewidth=5.0)
plot!(pl, xHistCGDA[1:1:end], yHistCGDA[1:1:end],
      marker=:pentagon, markersize=5.0, label="CGD", linestyle=:solid,
      linewidth=5.0)

savefig(pl, "./out/exp1/test1/bilinear_strong_alpha$(α)_eta$(η)_gamma$(γ)")


text:
using DifferentialEquations
using ModelingToolkit

function rober(du, u, p, t)
    y1, y2, y3 = u
    k1, k2, k3 = p

    du[1] = -k1*y1 + k3*y2*y3
    du[2] = k1*y1 - k2*y2^2. - k3*y2*y3
    du[3] = k2*y2^2.

end

prob = ODEProblem(rober,[1.0,0.0,0.0],(0.0,1e2),(0.04,3e2,1e4))
de = modelingtoolkitize(prob)

jac = @eval ModelingToolkit.generate_jacobian(de)[2]

jac(zeros(3,3),[0.0184367966813563, 7.51141018910664e-8, 0.9815631282045398],(0.04,0.,0.), (0.0,1e5) )

function jacobian(u, p, jac)
    len = length(u)
    j = zeros(len, len)
    jac(j, u, p, nothing)
    return j
end

jacobian([0.,1.,2.], [0.33,0.3,0.1], jac)


text:
using NLsolve

"Holds the coefficients of a projected Gauss-Legendre Runge-Kutta method."
struct CoefficientsPGLRK{T} <: AbstractCoefficients{T}
    @HeaderCoefficientsRK
    @CoefficientsRK

    P::Matrix{T}
    Q::Matrix{T}
    X::Matrix{T}
    W::Matrix{T}
    A::Matrix{T}

    function CoefficientsPGLRK{T}(name,o,s,a,b,c,P,X,W) where {T}
        @assert T <: Real
        @assert isa(name, Symbol)
        @assert isa(o, Integer)
        @assert isa(s, Integer)
        @assert s ≥ 2 "Number of stages must be ≥ 2"
        @assert s==size(a,1)==size(a,2)==length(b)==length(c)
        @assert s==size(P,1)==size(P,2)
        @assert s==size(X,1)==size(X,2)
        @assert s==size(W,1)==size(W,2)

        Q = inv(P)
        A = zero(a)
        B = zero(a)

        simd_mult!(B, W, Q)
        simd_mult!(A, P, B)

        new(name,o,s,a,b,c,P,Q,X,W,A)
    end
end

function CoefficientsPGLRK(name::Symbol, order::Int, a::Matrix{T}, b::Vector{T}, c::Vector{T}, P::Matrix{T}, X::Matrix{T}, W::Matrix{T}) where {T}
    CoefficientsPGLRK{T}(name, order, length(c), a, b, c, P, X, W)
end

Base.hash(tab::CoefficientsPGLRK, h::UInt) = hash(tab.o, hash(tab.s, hash(tab.a, hash(tab.b, hash(tab.c, hash(tab.P, hash(tab.Q, hash(tab.X, hash(tab.W, hash(tab.A, h))))))))))

Base.:(==)(tab1::CoefficientsPGLRK, tab2::CoefficientsPGLRK) = (tab1.o == tab2.o
                                                             && tab1.s == tab2.s
                                                             && tab1.a == tab2.a
                                                             && tab1.b == tab2.b
                                                             && tab1.c == tab2.c
                                                             && tab1.P == tab2.P
                                                             && tab1.Q == tab2.Q
                                                             && tab1.X == tab2.X
                                                             && tab1.W == tab2.W
                                                             && tab1.A == tab2.A)

Base.isequal(tab1::CoefficientsPGLRK{T1}, tab2::CoefficientsPGLRK{T2}) where {T1,T2} = (tab1 == tab2 && T1 == T2 && typeof(tab1) == typeof(tab2))

"Print Runge-Kutta coefficients."
function Base.show(io::IO, tab::CoefficientsPGLRK)
    print(io, "Projected Gauss-Legendre Runge-Kutta Coefficients ", tab.name, " with ", tab.s, " stages and order ", tab.o)
    print(io, "  a = ", tab.a)
    print(io, "  b = ", tab.b)
    print(io, "  c = ", tab.c)
    print(io, "  P = ", tab.P)
    print(io, "  Q = ", tab.Q)
    print(io, "  X = ", tab.X)
    print(io, "  W = ", tab.W)
    print(io, "  A = ", tab.A)
end

function getTableauPGLRK(coeff::CoefficientsPGLRK{T}, λ, a::Matrix{T}) where {T}
    a .= coeff.a .+ λ .* coeff.A
end

function getTableauPGLRK(coeff::CoefficientsPGLRK{T}, λ) where {T}
    coeff.a .+ λ .* coeff.A
end

"Parameters for right-hand side function of projected Gauss-Legendre Runge-Kutta methods."
mutable struct ParametersPGLRK{DT,TT,D,S,ET} <: Parameters{DT,TT}
    equ::ET
    tab::CoefficientsPGLRK{TT}
    Δt::TT

    h₀::DT

    t̅::TT
    t::TT

    q̅::Vector{DT}
    λ::DT

    function ParametersPGLRK{DT,TT,D,S,ET}(equ, tab, Δt) where {DT,TT,D,S,ET}
        new(equ, tab, Δt, zero(TT), zero(TT), zero(DT), zeros(DT,D), zero(DT))
    end
end

mutable struct IntegratorCachePGLRK{DT,D,S} <: ODEIntegratorCache{DT,D}
    h::DT
    λ::DT

    q̃::Vector{DT}
    ṽ::Vector{DT}
    s̃::Vector{DT}

    Q::Vector{Vector{DT}}
    V::Vector{Vector{DT}}
    Y::Vector{Vector{DT}}

    function IntegratorCachePGLRK{DT,D,S}() where {DT,D,S}
        # create temporary vectors
        q̃ = zeros(DT,D)
        ṽ = zeros(DT,D)
        s̃ = zeros(DT,D)

        # create internal stage vectors
        Q = create_internal_stage_vector(DT, D, S)
        V = create_internal_stage_vector(DT, D, S)
        Y = create_internal_stage_vector(DT, D, S)

        new(0, 0,
            q̃, ṽ, s̃,
            Q, V, Y)
    end
end

"""
Projected Gauss-Legendre Runge-Kutta integrator.

    Reference: LUIGI BRUGNANO, FELICE IAVERNARO, AND DONATO TRIGIANTE.
        ENERGY- AND QUADRATIC INVARIANTS–PRESERVING INTEGRATORS BASED
        UPON GAUSS COLLOCATION FORMULAE.
        SIAM J. NUMER. ANAL. Vol. 50, No. 6, pp. 2897–2916, 2012.
"""
struct IntegratorPGLRK{DT, TT, PT <: ParametersPGLRK{DT,TT},
                               ST <: NonlinearSolver{DT},
                               IT <: InitialGuessODE{DT,TT}, N, D, S} <: IntegratorPRK{DT,TT}
    params::PT
    solver::ST
    iguess::IT
    cache::IntegratorCachePGLRK{DT,D,S}

    function IntegratorPGLRK(N, params::ParametersPGLRK{DT,TT,D,S,ET}, solver::ST, iguess::IT, cache) where {DT, TT, D, S, ET, ST, IT}
        new{DT, TT, typeof(params), ST, IT, N, D, S}(params, solver, iguess, cache)
    end
end

function IntegratorPGLRK(equation::ODE{DT,TT,VT,HT,N}, tableau::CoefficientsPGLRK{TT}, Δt::TT) where {DT,TT,VT,HT,N}
    D = equation.d
    M = equation.n
    S = tableau.s

    # create params
    params = ParametersPGLRK{DT,TT,D,S,typeof(equation)}(equation, tableau, Δt)

    # create solver
    solver  = create_nonlinear_solver(DT, D*S, params)

    # create initial guess
    iguess = InitialGuessODE(get_config(:ig_interpolation), equation, Δt)

    # create cache
    cache = IntegratorCachePGLRK{DT,D,S}()

    # create integrator
    IntegratorPGLRK(N, params, solver, iguess, cache)
end

@inline nstages(integrator::IntegratorPGLRK{DT,TT,PT,ST,IT,N,D,S}) where {DT,TT,PT,ST,IT,N,D,S} = S

function update_params!(params::ParametersPGLRK, sol::AtomicSolutionODE)
    # set time for nonlinear solver and copy previous solution
    params.t̅  = sol.t
    params.t  = sol.t + params.Δt
    params.q̅ .= sol.q
end

function compute_stages!(x::Vector{ST}, cache::IntegratorCachePGLRK{ST}, params::ParametersPGLRK{DT,TT,D,S}) where {ST,DT,TT,D,S}
    compute_stages!(x, cache.Q, cache.V, cache.Y, cache.q̃, cache.ṽ, params)
end

function compute_stages!(x::Vector{ST}, Q::Vector{Vector{ST}}, V::Vector{Vector{ST}}, Y::Vector{Vector{ST}},
                                        q::Vector{ST}, y::Vector{ST},
                                        params::ParametersPGLRK{DT,TT,D,S}) where {ST,DT,TT,D,S}

    local tᵢ::TT

    # copy x to Y
    for i in 1:S
        for k in 1:D
            Y[i][k] = x[D*(i-1)+k]
        end
    end

    # compute Q=q̅+Δt*Y
    for i in 1:S
        for k in 1:D
            Q[i][k] = params.q̅[k] + params.Δt * Y[i][k]
        end
    end

    # compute V=v(T,Q)
    for i in 1:S
        tᵢ = params.t̅ + params.Δt * params.tab.c[i]
        params.equ.v(tᵢ, Q[i], V[i])
    end

    # compute y=B*V
    y .= 0
    for k in 1:D
        for j in 1:S
            y[k] += params.tab.b[j] * V[j][k]
        end
    end

    # compute q=q̅+Δt*y
    q .= params.q̅ .+ params.Δt .* y
end

"Compute stages of projected Gauss-Legendre Runge-Kutta methods."
@generated function function_stages!(x::Vector{ST}, b::Vector{ST}, params::ParametersPGLRK{DT,TT,D,S}) where {ST,DT,TT,D,S}

    cache = IntegratorCachePGLRK{ST,D,S}()

    quote
        compute_stages!(x, $cache, params)

        a = getTableauPGLRK(params.tab, params.λ)

        # compute b = [Y-AV]
        for i in 1:S
            for k in 1:D
                b[D*(i-1)+k] = $cache.Y[i][k]
                for j in 1:S
                    b[D*(i-1)+k] -= a[i,j] * $cache.V[j][k]
                end
            end
        end
    end
end

function bisection(f::Function, λmin::DT, λmax::DT;
                   xtol::AbstractFloat=get_config(:nls_atol),
                   ftol::AbstractFloat=get_config(:nls_atol),
                   maxiter::Integer=get_config(:nls_nmax)) where {DT <: Number}
    a = λmin
    b = λmax
    fa = f(a)
    fb = f(b)

    # fa*fb ≤ 0 || error("Either no or multiple real roots in [λmin,λmax]")

    # local j = 0
    local λ = zero(DT)
    local fλ= zero(DT)

    for i in 1:maxiter
        # j += 1
        λ  = (a+b)/2
        fλ = f(λ)

        !isapprox(fλ, 0, atol=ftol) || break

        if fa*fλ > 0
            a  = λ  # Root is in the right half of [a,b].
            fa = fλ
        else
            b = λ  # Root is in the left half of [a,b].
        end

        abs(b-a) > xtol || break
    end

    # println(j, " bisection iterations, λ=", λ, ", f(λ)=", fλ, ", ftol=", ftol, ", abs(b-a)=", abs(b-a), ", xtol=", xtol)

    # i != maxiter || error("Max iteration number exceeded")

    return λ
end

function function_hamiltonian!(λ::Number, int::IntegratorPGLRK{DT,TT}) where {DT,TT}
    # copy λ to integrator parameters
    int.params.λ = λ

    # call nonlinear solver
    solve!(int.solver)

    # print solver status
    print_solver_status(int.solver.status, int.solver.params)

    # check if solution contains NaNs or error bounds are violated
    check_solver_status(int.solver.status, int.solver.params)

    # compute vector fields at internal stages
    compute_stages!(int.solver.x, int.cache, int.params)

    # compute h
    int.cache.h = int.params.equ.h(int.params.t, int.cache.q̃)

    # compute and return h₀-h
    return int.params.h₀ - int.cache.h
end

function function_hamiltonian!(λ::Vector, int::IntegratorPGLRK{DT,TT}) where {DT,TT}
    [function_hamiltonian!(λ[1], int)]
end

function initialize!(int::IntegratorPGLRK, sol::AtomicSolutionODE)
    sol.t̅ = sol.t - timestep(int)

    equation(int).v(sol.t, sol.q, sol.v)

    initialize!(int.iguess, sol.t, sol.q, sol.v,
                            sol.t̅, sol.q̅, sol.v̅)

    int.params.h₀ = int.params.equ.h(sol.t, sol.q)
end

function initial_guess!(int::IntegratorPGLRK{DT,TT}, sol::AtomicSolutionODE{DT,TT}) where {DT,TT}
    for i in eachstage(int)
        evaluate!(int.iguess, sol.q, sol.v, sol.q̅, sol.v̅,
                              int.cache.q̃, int.cache.ṽ,
                              tableau(int).c[i])

        for k in eachdim(int)
            int.solver.x[ndims(int)*(i-1)+k] = int.cache.ṽ[k]
        end
    end
end

"Integrate ODE with projected Gauss-Legendre Runge-Kutta integrator."
function integrate_step!(int::IntegratorPGLRK{DT,TT}, sol::AtomicSolutionODE{DT,TT}) where {DT,TT}
    # update nonlinear solver parameters from cache
    update_params!(int.params, sol)

    # compute initial guess
    initial_guess!(int, sol)

    # reset solution
    reset!(sol, timestep(int))

    # determine parameter λ
    λmin = -0.2^nstages(int)
    λmax = +0.2^nstages(int)

    int.params.λ = bisection(λ -> function_hamiltonian!(λ, int), λmin, λmax;
                    xtol=abs(λmax-λmin)*get_config(:nls_atol),
                    ftol=int.params.h₀*get_config(:nls_atol))
    # int.params.λ = nlsolve(λ -> function_hamiltonian!(λ, int), [zero(DT)];
    #             xtol=abs(λmax-λmin)*get_config(:nls_atol),
    #             ftol=int.params.h₀*get_config(:nls_atol)).zero[1]
    # println(int.params.λ)

    # compute final update
    update_solution!(sol.q, int.cache.V, tableau(int).b, timestep(int))

    # copy solution to initial guess
    update!(int.iguess, sol.t, sol.q, sol.v)
end


text:
using httpserver.jl

implicit server

do bind()
  # tbd
end


text:
module test_transientsol
using VoronoiFVM

function make_transientsol(;n=10,M=5,N=100,in_memory=true)
    makevec(k)=[ k+i*j for i=1:M,j=1:N]
    sol=TransientSolution(0,makevec(0),in_memory=in_memory)
    for k=1:n
        append!(sol,k,makevec(k))
    end
    sol
end

function test()
    msol=make_transientsol(in_memory=true)
    dsol=make_transientsol(in_memory=false)
    msol==dsol
end

end


text:
# Integrate the Rocket Landing example with MonteCarlo

# RUN:
# ] activate TrajOptSOCPs

using TrajOptSOCPs
using LinearAlgebra, SparseArrays
using Dates

include("../plotting/batchPlots.jl")
include("../plotting/batchPlots_MonteCarlo.jl")
pyplot()
runplots = true
saveplots = true

println("\n--------------------------------------------")
println("          Setting Up Rocket Landing ")
println("--------------------------------------------")

# We will keep the same rocket in all runs
mass = 1
isp = 1
grav = [0; 0; -9.81]
deltaTime = 1
rocket = rocket_simple(mass, isp, grav, deltaTime)

# Then ending point is always a soft landing at the origin.
# But, we want to vary the starting point.
minStateStart = [-10.0; -5.0; 20.0; -4.0; -4.0; -8.0]
maxStateStart = [  0.0;  5.0; 30.0;  4.0;  4.0;  0.0]
const numRuns = 50

# Stores each randomly generated starting case
startList = TrajOptSOCPs.montecarlo(minStateStart, maxStateStart, numRuns)

# Maintain the same end point
const rocketEnd = [0.0; 0.0; 0.0; 0.0; 0.0; 0.0]

# Initialize the hover control, which is the base control to prevent freefall.
uHover = mass * grav

# Number of time steps to discretize the trajectory
const NSteps = 60

# The following trajectories do not change with starting condition:
# ----------------------------
# Create the Ground Constraints
groundConstraint = makeGroundConstraint(NSteps, size(grav, 1), size(grav, 1))
groundLambda = zeros(size(groundConstraint.b, 1))

# Create the Max Thrust Constraints
thrustMax = 20.0
maxThrustConstraint = makeMaxThrustConstraint(NSteps, size(grav, 1), thrustMax)
maxThrustLambda = zeros(size(maxThrustConstraint.indicatorList))

# Create the Max Thrust Angle Constraints
thrustAngleMax = 10.0 # Deg
maxAngleConstraint = makeMaxAngleConstraint(NSteps, size(grav, 1),
                                                thrustAngleMax, true)
maxAngleLambda = zeros(size(maxAngleConstraint.indicatorList))

# Set up the LQR matrices
const lqrQMat = 0.0 * Diagonal(I, size(rocketEnd, 1))
const lqrRMat = 0.025 * Diagonal(I, Int64(size(rocketEnd, 1) / 2))

# Initialize the Augmented Lagrangian penalty on the constraints
penaltyStart = 1.0

# Next we select resonable solver parameters
currSolveParams = solverParams(0.1, 0.5,
                                8, 10,
                                10^-4,
                                10, 10^6,
                                2.5, 2, 0.2, 0.2, 0.4,
                                TrajOptSOCPs.NONE)
# ----------------------------

solParamPrint(currSolveParams)
println()
println()

# Arrays to store the results
trajPDFinals = []
trajPFinals = []

dynVio = []
constVio = []

# Now we loop through the different starting conditions and solve the
# optimization problem
for (ind, rocketStart) in enumerate(startList)
    # Initialize the trajectory with a line
    initTraj = initializeTraj(rocketStart, rocketEnd, uHover, uHover, NSteps)

    # Use a Linear Quadratic Regulator as the cost function
    costFun = makeLQR_TrajReferenced(lqrQMat, lqrRMat, NSteps, initTraj)

    # Create the Dynamics Constraints
    ADyn, BDyn = rocketDynamicsFull(rocket, rocketStart, rocketEnd, NSteps)
    dynConstraint = AL_AffineEquality(ADyn, BDyn)
    lambdaInit = -1 * ones(size(BDyn))

    # Create the constraint manager
    cMRocket = constraintManager_Dynamics(
                [groundConstraint, maxThrustConstraint, maxAngleConstraint],
                [groundLambda, maxThrustLambda, maxAngleLambda],
                dynConstraint, lambdaInit
                )

    # Initialize the primal-dual vector
    initTrajPD = [initTraj; lambdaInit]

    # Equiped with the constraint term and the objective term, I now build the
    # Augmented Lagrangian
    alRocket = augLag(costFun, cMRocket, penaltyStart)

    println("--------------------------------------------")
    println("             Beginning Solve $ind           ")
    println("--------------------------------------------")

    trajLambdaSolved, resArr = ALPrimalNewtonMain(initTrajPD, alRocket,
                                            currSolveParams)
    trajPD = parsePrimalDualVec(trajLambdaSolved[end], size(initTraj, 1))
    trajP = trajPD.primals

    push!(trajPDFinals, trajPD)
    push!(trajPFinals, trajP)

    # Next get the final constraint violations
    # (1) Dynamics Vioaltion
    push!(dynVio, safeNorm(norm(evalAffineEq(cMRocket, trajP))))
    # (2) Other Constraint Violation
    push!(constVio, safeNorm(evalConstraints(cMRocket, trajP, penaltyStart)))

end

if runplots
    println("--------------------------------------------")
    println("             Beginning Plotting             ")
    println("--------------------------------------------")

    allPlots = batchPlotMonteCarlo(trajPFinals, dynVio, constVio,
                            size(grav, 1), penaltyStart,
                            thrustMax, thrustAngleMax)
end

if runplots && saveplots
    curr = Dates.format(now(), "yyyymmdd_HH-MM-SS")
    header = "MonteCarlo_$(numRuns)Runs_" * curr * "_"
    saveBulk(allPlots, header)
end

max_stars_repo_name: danineamati/TrajOptSOCPs