text:
#### SVD, i.e. UDV decomposition
decompose_udv!(A::AbstractMatrix{<:Number}) = LinearAlgebra.LAPACK.gesvd!('A','A',A) # was gesdd! at some point
# decompose_udv!(A::AbstractMatrix{<:Number}) = LinearAlgebra.LAPACK.gesdd!('A',A) # was gesdd! at some point
# using JacobiSVD # experimental package
# decompose_udv!(A::AbstractMatrix{<:Number}) = JacobiSVD.gesvj!('G','U','V', A) # was gesdd! at some point
decompose_udv(A::AbstractMatrix{T}) where T<:Number = decompose_udv!(copy(A))

#### QR, i.e. UDT decomposition
function decompose_udt(A::AbstractMatrix{C}) where C<:Number
  F = qr(A, Val(true))
  p = F.p
  @views p[p] = collect(1:length(p))
  # D = abs.(real(diag(triu(R))))
  D = abs.(real(diag(F.R)))
  T = (Diagonal(1 ./ D) * F.R)[:, p]
  return Matrix(F.Q), D, T
end

function decompose_udt!(A::AbstractMatrix{C}, D) where C<:Number
  n = length(D)
  F = qr!(A, Val(true))
  R = F.R # F.R is of regular matrix type

  @views F.p[F.p] = 1:n

  @inbounds for i in 1:n
    D[i] = abs(real(R[i,i]))
  end

  # This (very!) strangely increases the runtime of measure_tdgfs! by a ton....
  # @inbounds for i in diagind(R)
  #   R[i] = sign(real(R[i]))
  # end

  lmul!(Diagonal(1 ./ D), R)

  return Matrix(F.Q), R[:, F.p] # Q, (D is modified in-place), T  # was full(F.Q) before upgrade
end

#### Other
# See https://discourse.julialang.org/t/asymmetric-speed-of-in-place-sparse-dense-matrix-product/10256/3
import LinearAlgebra.mul!
function mul!(C::StridedMatrix, X::StridedMatrix, A::SparseMatrixCSC)
    mX, nX = size(X)
    nX == A.m || throw(DimensionMismatch())
    fill!(C, zero(eltype(C)))
    rowval = A.rowval
    nzval = A.nzval
    @inbounds for  col = 1:A.n, k=A.colptr[col]:(A.colptr[col+1]-1)
        ki=rowval[k]
        kv=nzval[k]
        for multivec_row=1:mX
            C[multivec_row, col] += X[multivec_row, ki] * kv
        end
    end
    C
end

myrdiv!(dest, a, b) = copyto!(dest, adjoint(adjoint(b) \ adjoint(a)))

LinearAlgebra.Diagonal(I::UniformScaling{T}) where T <: Number = I.λ

function eye!(x)
    @assert ndims(x) == 2
    @assert size(x,1) == size(x,2)

    o = one(eltype(x))
    z = zero(eltype(x))

    for j in 1:size(x,2)
        for i in 1:size(x,1)
            @inbounds x[i, j] = i == j ? o : z
        end
    end
    x
end

"""
Kronecker Delta
"""
@inline kd(a,b) = a == b

##############################################################
#
#                   QR / UDT
#
##############################################################
"""
Safely multiply two UDT decompositions
"""
function multiply_safely(Ul,Dl,Tl, Ur,Dr,Tr)
  mat = Tl * Ur
  lmul!(Diagonal(Dl), mat)
  rmul!(mat, Diagonal(Dr))
  U, D, T = decompose_udt(mat)
  return Ul*U, D, T*Tr
end

"""

  inv_udt(U, D, T) -> result

Calculate (UDT)^(-1).
"""
function inv_udt(U,D,T)
  res = similar(U)
  inv_udt!(res, U, D, T)
  res
end

"""

  inv_udt!(res, U, D, T) -> nothing

Calculate (UDT)^(-1) and store result in `res`
"""
function inv_udt!(res, U,D,T)
  tmp = similar(U)
  ldiv!(tmp, lu(T), Diagonal(1 ./ D))
  mul!(res, tmp, U')
  nothing
end

"""

  inv_one_plus_udt(U, D, T) -> result

Stable calculation of [1 + UDT]^(-1):

  * Use one intermediate UDT decomposition.

Faster but less accurate than the scalettar approach.

Consider `inv_one_plus_udt!` as an efficient (not one-to-one) replacement.
"""
function inv_one_plus_udt(U,D,T)
  @warn "Calling somewhat inefficient and potentially inaccurate `inv_one_plus_udt`"

  m = U' / T
  m[diagind(m)] .+= D
  u,d,t = decompose_udt(m)
  u = U * u
  t = t * T
  ldiv!(m, lu!(t), Diagonal(1 ./ d))
  m * u'
end

"""

  inv_one_plus_udt!(mc, res, U, D, T) -> nothing

Stable calculation of [1 + UDT]^(-1):

  * Use one intermediate UDT decomposition.

Uses preallocated memory in `mc`. Writes the result into `res`.

Much faster (~50%) than `inv_one_plus_udt_scalettar!` but less accurate.
"""
function inv_one_plus_udt!(mc, res, U,D,T)
  @warn "Calling potentially inaccurate `inv_one_plus_udt!`"

  d = mc.s.d
  u = mc.s.tmp
  t = mc.s.tmp2

  m = U' / T
  m[diagind(m)] .+= D

  utmp,ttmp = decompose_udt!(m, d)
  mul!(u, U, utmp)
  mul!(t, ttmp, T)

  ldiv!(m, lu!(t), Diagonal(1 ./ d))

  mul!(res, m, u')
  nothing
end

"""

  inv_one_plus_two_udts!(mc, U, D, T, Ul, Dl, Tl, Ur, Dr, Tr) -> nothing

Stable calculation of [1 + UlDlTl(UrDrTr)^†]^(-1).

Uses preallocated memory in `mc`. Writes the result into `U`, `D`, and `T`.
"""
function inv_one_plus_two_udts!(mc, U,D,T, Ul,Dl,Tl, Ur,Dr,Tr)
  s = mc.s
  tmp = mc.s.tmp
  tmp2 = mc.s.tmp2
  tmp3 = mc.s.curr_U

  mul!(tmp, Tl, adjoint(Tr))
  rmul!(tmp, Diagonal(Dr))
  lmul!(Diagonal(Dl), tmp)
  U1, T1 = decompose_udt!(tmp, s.D)

  mul!(tmp3, Ul, U1)
  mul!(tmp2, T1, adjoint(Ur))
  mul!(tmp, adjoint(tmp3), inv(tmp2))

  tmp .+= Diagonal(s.D)

  u, t = decompose_udt!(tmp, D)

  mul!(tmp, t, tmp2)
  copyto!(U, inv(tmp))

  mul!(tmp, tmp3, u)
  copyto!(T, adjoint(tmp))

  D .= 1 ./ D

  nothing
end

"""

  inv_one_plus_two_udts!(mc, res, Ul, Dl, Tl, Ur, Dr, Tr) -> nothing

Stable calculation of [1 + UlDlTl(UrDrTr)^†]^(-1).

Uses preallocated memory in `mc`. Writes the result into `res`.
"""
function inv_one_plus_two_udts!(mc, res, Ul,Dl,Tl, Ur,Dr,Tr)
  inv_one_plus_two_udts!(mc, mc.s.U, mc.s.d, mc.s.T, Ul, Dl, Tl, Ur, Dr, Tr)
  rmul!(mc.s.U, Diagonal(mc.s.d))
  mul!(res, mc.s.U, mc.s.T)
  nothing
end

"""

  inv_one_plus_udt_scalettar!(mc, res, Ua, Da, Ta, Ub, Db, Tb) -> nothing

Stable calculation of [1 + UDT]^(-1):

  * Separate scales larger and smaller than unity
  * Use two intermediate UDT decompositions.

Uses preallocated memory in `mc`. Writes the result into `res`.
"""
function inv_one_plus_udt_scalettar!(mc, res, U,D,T)
  d = mc.s.d
  r = mc.s.tmp
  l = mc.s.tmp2

  Dp = max.(D,1.)
  Dm = min.(D,1.)

  Dp .\= 1
  Dpinv = Dp # renaming

  ldiv!(l, lu(T), Diagonal(Dpinv)) # Don't use lu!(T) because T is input

  mul!(r, U, Diagonal(Dm))
  r .+= l

  u, t = decompose_udt!(r, d)

  ldiv!(r, lu!(t), Diagonal(1 ./ d))
  mul!(l, r, u')

  lmul!(Diagonal(Dpinv), l)
  u, t = decompose_udt!(l, d)

  ldiv!(l, lu(T), u)

  rmul!(l, Diagonal(d))
  mul!(res, l, t)
  nothing
end

"""

  inv_sum_udts(Ua, Da, Ta, Ub, Db, Tb) -> U, D, T

Stable calculation of [UaDaTa + UbDbTb]^(-1):

  * Use one intermediate UDT decompositions.

Faster but less accurate than `inv_sum_udts_scalettar`.

Consider `inv_sum_udts!` as an efficient (not one-to-one) replacement.
"""
function inv_sum_udts(Ua,Da,Ta,Ub,Db,Tb)
  @warn "Calling somewhat inefficient and potentially inaccurate `inv_sum_udts`"

  m1 = Ta / Tb
  lmul!(Diagonal(Da), m1)

  m2 = Ua' * Ub
  rmul!(m2, Diagonal(Db))

  u,d,t = decompose_udt(m1 + m2)

  mul!(m1, Ua, u)
  mul!(m2, t, Tb)

  return inv(m2), 1 ./ d, m1'
end

"""

  inv_sum_udts!(mc, res, Ua, Da, Ta, Ub, Db, Tb) -> nothing

Stable calculation of [UaDaTa + UbDbTb]^(-1):

  * Use one intermediate UDT decompositions.

Uses preallocated memory in `mc`. Writes the result into `res`.

Much faster (~40%) than `inv_sum_udts_scalettar!` but less accurate.
"""
function inv_sum_udts!(mc, res, Ua,Da,Ta,Ub,Db,Tb)
  @warn "Calling potentially inaccurate `inv_sum_udts!`"

  d = mc.s.d
  m1 = mc.s.tmp
  m2 = mc.s.tmp2
  tmp = mc.s.U

  m1 = Ta / Tb
  lmul!(Diagonal(Da), m1)

  mul!(m2, Ua', Ub)
  rmul!(m2, Diagonal(Db))

  u,t = decompose_udt!(m1 + m2, d)

  mul!(m1, Ua, u)
  mul!(m2, t, Tb)

  ldiv!(tmp, lu!(m2), Diagonal(1 ./ d))
  mul!(res, tmp, m1')

  nothing
end

"""

  inv_sum_udts_scalettar(Ua, Da, Ta, Ub, Db, Tb) -> U, D, T

Stable calculation of [UaDaTa + UbDbTb]^(-1):

  * Separate scales larger and smaller than unity
  * Use two intermediate UDT decompositions.

Consider `inv_sum_udts_scalettar!` as an efficient (not one-to-one) replacement.
"""
function inv_sum_udts_scalettar(Ua, Da, Ta, Ub, Db, Tb)
    @warn "Calling somewhat inefficient `inv_sum_udts_scalettar`"

    d=length(Da)

    # separating scales larger and smaller than unity
    Dap = max.(Da,1.)
    Dam = min.(Da,1.)
    Dbp = max.(Db,1.)
    Dbm = min.(Db,1.)

    # mat1 = Dam * Vda * Vdb' / Dbp
    mat1 = Ta * inv(Tb)
    @inbounds for j in 1:d, k in 1:d
        mat1[j,k]=mat1[j,k] * Dam[j]/Dbp[k]
    end

    # mat2 = 1/(Dap) * Ua' * Ub * Dbm
    mat2 = Ua' * Ub
    @inbounds for j in 1:d, k in 1:d
        mat2[j,k]=mat2[j,k] * Dbm[k]/Dap[j]
    end

    #mat1 = mat1 + mat2
    mat1 = mat1 + mat2

    # decompose mat1: U, D, T
    U, D, T = decompose_udt(mat1)

    # invert inner part: mat1 = (U D T)^(-1) = mat1^(-1)
    # was UDT_to_mat!(mat1, U, D, T, invert=true)
    lmul!(Diagonal(D), T)
    ldiv!(mat1, lu!(T), U') # mat1 = T \ (U')

    # mat1 = 1/Dbp * mat1 /Dap
    @inbounds for j in 1:d, k in 1:d
        mat1[j,k]=mat1[j,k] / Dbp[j] / Dap[k]
    end

    #mat1 = U D T
    U, D, T = decompose_udt(mat1)

    # U = Tb^(-1) * U , T = T * Ua'
    mul!(mat1, inv(Tb), U)
    mul!(mat2, T, Ua')
    U = mat1
    T = mat2

    return U, D, T
end

"""

  inv_sum_udts_scalettar!(mc, res, Ua, Da, Ta, Ub, Db, Tb) -> nothing

Stable calculation of [UaDaTa + UbDbTb]^(-1):

  * Separate scales larger and smaller than unity
  * Use two intermediate UDT decompositions.

Uses preallocated memory in `mc`. Writes the result into `res`.
"""
function inv_sum_udts_scalettar!(mc, res, Ua, Da, Ta, Ub, Db, Tb)
    # optimization rounds: 1
    mat1 = mc.s.tmp
    mat2 = mc.s.U
    D = mc.s.d
    to = mc.a.to

    d=length(Da)

    # separating scales larger and smaller than unity
    Dap = max.(Da,1.)
    Dam = min.(Da,1.)
    Dbp = max.(Db,1.)
    Dbm = min.(Db,1.)

    # mat1 = Dam * Vda * Vdb' / Dbp
    mat1 = Ta / Tb
    @inbounds for j in 1:d, k in 1:d
        mat1[j,k]=mat1[j,k] * Dam[j]/Dbp[k]
    end

    # mat2 = 1/(Dap) * Ua' * Ub * Dbm
    mul!(mat2, Ua', Ub)
    @inbounds for j in 1:d, k in 1:d
        mat2[j,k]=mat2[j,k] * Dbm[k]/Dap[j]
    end

    # mat1 = mat1 + mat2
    mat1 .+= mat2

    # decompose mat1: U, D, T
    U, T = decompose_udt!(mat1, D)

    # invert and combine inner part: mat1 = (U D T)^(-1)
    lmul!(Diagonal(D), T)
    ldiv!(mat1, lu!(T), U') # mat1 = T \ (U')

    # mat1 = 1/Dbp * mat1 /Dap
    @inbounds for j in 1:d, k in 1:d
        mat1[j,k]=mat1[j,k] / Dbp[j] / Dap[k]
    end

    #mat1 = U D T
    U, T = decompose_udt!(mat1, D)

    # U = Tb^(-1) * U , T = T * Ua'
    ldiv!(mat1, lu!(Tb), U) # mat1 = Tb \ U
    mul!(mat2, T, Ua')
    U = mat1
    T = mat2

    # combine UDT into res
    rmul!(U, Diagonal(D))
    mul!(res, U, T)
    nothing
end

"""
  UDT_to_mat!(mat, U, D, T[; invert=false]) -> nothing

Combine given UDT to a matrix. Result will be stored in `mat`.

If `invert = true` (UDT)^(-1) is calculated.

Note that if parts of the input UDT may be overwritten it's
more efficient to inline this function and remove unnecessary copies!
"""
function UDT_to_mat!(mat, U, D, T; invert=false)
    # TODO: Check where this is called and inline manually for better speed. (remove copy())
    @warn "UDT_to_mat! probably shouldn't be called here" # TODO remove once checked
    if !invert
        mat1 = copy(U)
        rmul!(mat1, Diagonal(D))
        mul!(mat,mat1,T)
    else # (DT)^-1 * U^dagger
        mat1 = copy(T)
        lmul!(Diagonal(D), mat1)
        mat .= mat1 \ U'
    end
    nothing
end

function UDT_to_mat(U, D, T; kw...)
  res = copy(U)
  UDT_to_mat!(res, U, D, T; kw...)
  res
end

##############################################################
#
#                   SVD / UDV(t)
#
##############################################################
# multiplies two UDVds -> UDVd
function multiply_safely_udv(Ul,Dl,Vdl,Ur,Dr,Vdr)
  tmp = adjoint(Vdl) * Ur
  rmul!(tmp, Diagonal(Dr))
  lmul!(Diagonal(Dl), tmp)
  U, D, Vd = decompose_udv!(tmp)
  U = Ul * U
  Vd = Vd * Vdr
  return U, D, Vd
end

# Calculates (UDVd)^-1, where U, D, Vd come from SVD decomp.
function inv_udv(U,D,Vd)
  m = copy(Vd')
  rmul!(m, Diagonal(1 ./ D))
  res = similar(m)
  mul!(res, m, U')
  res
end

# Calculates (UDVd)^-1, where U, D, Vd come from SVD decomp.
function inv_udv!(res, U,D,Vd)

  # copy here isn't necessary but Vd would be overwritten if we drop it
  # also, m could be preallocated in stack if necessary
  m = copy(Vd')
  rmul!(m, Diagonal(1 ./ D))
  mul!(res, m, U')
  nothing
end

# Calculates (1 + UDVd)^-1, where U, D, Vd come from SVD decomp.
# !! Breaks down for large spread in D (i.e. low temperatures).
function inv_one_plus_udv(U,D,Vd)
  inner = copy(Vd')
  inner .+= U * Diagonal(D)
  I = decompose_udv!(inner)
  u = copy(adjoint(I[3] * Vd))
  d = 1 ./ I[2]
  vd = adjoint(I[1])

  rmul!(u,Diagonal(d))
  u * vd
end

# same as inv_one_plus_udt but separating both U AND Vd from D
# !! Breaks down for large spread in D (i.e. low temperatures). Slightly better than normal version.
function inv_one_plus_udv_alt(U,D,Vd)
  inner = copy((Vd*U)')
  inner[diagind(inner)] .+= D
  u, d, vd = decompose_udv!(inner)

  t1 = adjoint(vd*Vd)
  t2 = adjoint(U*u)
  rmul!(t1, Diagonal(1 ./ d))
  t1*t2
end

# Calculates (1 + UDVd)^-1, where U, D, Vd come from SVD decomp.
# More controlled handling of scales, however also slower.
# This is Yoni's `stable_invert_B`
# Seems to work even for large spread in D (i.e. low temperatures)!
function inv_one_plus_udv_scalettar(U,D,Vd)
  Dp = max.(D,1.)
  Dm = min.(D,1.)
  Dpinv = 1 ./ Dp

  l = copy(Vd')
  rmul!(l, Diagonal(Dpinv))

  r = copy(U)
  rmul!(r, Diagonal(Dm))

  u, d, vd = decompose_udv!(l+r)

  m = inv_udv(u,d,vd)
  lmul!(Diagonal(Dpinv), m)
  u, d, vd = decompose_udv!(m)

  mul!(m, Vd', u)
  # return m, d, vd
  rmul!(m, Diagonal(d))
  m*vd
end

# TODO: Optimize!
# I only made the function overwrite res. Otherwise it's unchanged compared to the one above.
function inv_one_plus_udv_scalettar!(mc, res, U,D,Vd)
  # all similars here could go into mc.s
  Dp = similar(D)
  Dm = similar(D)
  l = similar(Vd)
  l .= Vd'
  r = similar(U)
  r .= U
  # u, d, vd below could be preallocated, but does this matter for speed?
  # m could be preallocated, I guess

  Dp .= max.(D, 1)
  Dm .= min.(D, 1)

  Dp .\= 1 # Dp now Dpinv!!!

  rmul!(l, Diagonal(Dp))
  rmul!(r, Diagonal(Dm))

  u, d, vd = decompose_udv!(l + r)

  m = inv_udv(u,d,vd) # TODO: optimize
  lmul!(Diagonal(Dp), m)
  u, d, vd = decompose_udv!(m)

  mul!(m, Vd', u)
  # return m, d, vd
  rmul!(m, Diagonal(d))
  # res .= m*vd
  mul!(res, m, vd)
  nothing
end

function UDV_to_mat!(mat, U, D, Vd; invert=false)
    if !invert
        mat1 = copy(U)
        rmul!(mat1, Diagonal(D))
        mul!(mat, mat1, Vd)
    else #V D^(-1) Ud = (D^-1 *Vd)^(dagger) *Ud
        mat1 = copy(Vd)
        lmul!(Diagonal(1 ./ D), mat1)
        mul!(mat, mat1', U')
    end
    nothing
end

function UDV_to_mat(U, D, Vd; kw...)
  res = copy(U)
  UDV_to_mat!(res, U, D, Vd; kw...)
  res
end

# Calculates (UaDaVda + UbDbVdb)^-1
function inv_sum_udvs(Ua, Da, Vda, Ub, Db, Vdb)

    d=length(Da)

    # separating scales larger and smaller than unity
    Dap = max.(Da,1.)
    Dam = min.(Da,1.)
    Dbp = max.(Db,1.)
    Dbm = min.(Db,1.)

    # mat1 = Dam * Vda * Vdb' / Dbp

    mat1 = Vda * adjoint(Vdb)
    for j in 1:d, k in 1:d
        mat1[j,k]=mat1[j,k] * Dam[j]/Dbp[k]
    end

    # mat2 = 1/(Dap) * Ua' * Ub * Dbm
    mat2 = adjoint(Ua) * Ub
    for j in 1:d, k in 1:d
        mat2[j,k]=mat2[j,k] * Dbm[k]/Dap[j]
    end

    #mat1 = mat1 + mat2
    mat1 = mat1 + mat2

    # decompose mat1: U, D, Vd
    U, D, Vd = decompose_udv!(mat1)

    # invert inner part: mat1 = (U D Vd)^(-1) = mat1^(-1)
    UDV_to_mat!(mat1, U, D, Vd, invert=true)

    # mat1 = 1/Dbp * mat1 /Dap
    for j in 1:d, k in 1:d
        mat1[j,k]=mat1[j,k] / Dbp[j] / Dap[k]
    end

    #mat1 = U D Vd
    U, D, Vd = decompose_udv!(mat1)

    # U = Vdb' * U , Vd = Vd * Ua'
    mul!(mat1, Vdb', U)
    mul!(mat2, Vd, Ua')
    U = mat1
    Vd = mat2

    return U, D, Vd
end

# Calculates (UaDaVda + UbDbVdb)^-1
# TODO: Optimize!
# I only made the function overwrite res. Otherwise it's unchanged compared to the one above.
function inv_sum_udvs!(mc, res, Ua, Da, Vda, Ub, Db, Vdb)

    d=length(Da)

    #DXp = max (X%D, 1) ,DXm = min (X%D, 1) and similarly for Y.
    Dap = max.(Da,1.)
    Dam = min.(Da,1.)
    Dbp = max.(Db,1.)
    Dbm = min.(Db,1.)

    #mat1= DXm * X%Vd * (Y%Vd)^dagger /DYp
    mat1 = Vda * adjoint(Vdb)
    for j in 1:d, k in 1:d
        mat1[j,k]=mat1[j,k] * Dam[j]/Dbp[k]
    end

    #mat2 = 1/(DXp) * (X%U)^dagger * Y%U * DYm
    mat2 = adjoint(Ua) * Ub
    for j in 1:d, k in 1:d
        mat2[j,k]=mat2[j,k] * Dbm[k]/Dap[j]
    end

    #mat1 = mat1+mat2
    mat1 = mat1 + mat2

    #invert mat1: mat1=mat1^(-1)
    U, D, Vd = decompose_udv!(mat1)
    UDV_to_mat!(mat1, U, D, Vd, invert=true)

    #mat1 = 1/DYp * mat1 /DXp
    for j in 1:d, k in 1:d
        mat1[j,k]=mat1[j,k] / Dbp[j] / Dap[k]
    end

    #mat1 = U D Vd
    U, D, Vd = decompose_udv!(mat1)

    # U = (Y%Vd)^dagger * U , Vd = Vd * (X%U)^dagger
    mul!(mat1,adjoint(Vdb),U)
    mul!(mat2,Vd,adjoint(Ua))
    U=mat1
    Vd=mat2

    # return U, D, Vd
    UDV_to_mat!(res, U, D, Vd, invert=false)
    nothing
end

max_stars_repo_name: crstnbr/dqmc

max_issues_repo_name: crstnbr/dqmc

max_forks_repo_name: crstnbr/dqmc


text:
function posre(topdir)
 # posre.itp file
  file   = open("$topdir/posre.itp","r+")
  file_new = open("./posre.itp","w")
  for line in eachline(file)
    println(file_new,line)
  end
  close(file)
  close(file_new)
end

### functoin to print mdp_files

function mdp_files_solute()

  # minimization
  file = open("mim.mdp", "w")
  println(file,"""
  ;Parameters describing what to do, when to stop and what to save
  integrator  = steep         ; Algorithm (steep = steepest descent minimization)
  emtol       = 100.0        ; Stop minimization when the maximum force < 100.0 kJ/mol/nm
  emstep      = 0.01          ; Minimization step size
  nsteps      = 50000         ; Maximum number of (minimization) steps to perform

  ; Parameters describing how to find the neighbors of each atom and how to calculate the interactions
  nstlist         = 1         ; Frequency to update the neighbor list and long range forces
  cutoff-scheme   = Verlet    ; Buffered neighbor searching
  ns_type         = grid      ; Method to determine neighbor list (simple, grid)
  coulombtype     = PME       ; Treatment of long range electrostatic interactions
  rcoulomb        = 1.0       ; Short-range electrostatic cut-off
  rvdw            = 1.0       ; Short-range Van der Waals cut-off
  pbc             = xyz       ; Periodic Boundary Conditions in all 3 dimensions
  """)
  close(file)

  # NVT equilibration
  file=open("nvt.mdp","w")
  println(file,"""
  ; title                   = OPLS Ubiquitin solvated by water with a ionicliquid being a cosolute NVT equilibration
  ;define                  = -DPOSRES  ; position restrain the protein
  ; Run parameters
  integrator              = md        ; leap-frog integrator
  nsteps                  = 500000     ; 2 * 50000 = 1000 ps
  dt                      = 0.002     ; 2 fs
  ; Output control
  nstxout                 = 5000       ; save coordinates every 1.0 ps
  nstvout                 = 5000       ; save velocities every 1.0 ps
  nstenergy               = 5000       ; save energies every 1.0 ps
  nstlog                  = 5000       ; update log file every 1.0 ps
  ; Bond parameters

  continuation            = no        ; first dynamics run
  constraint_algorithm    = lincs     ; holonomic constraints
  constraints             = h-bonds   ; bonds involving H are constrained
  lincs_iter              = 1         ; accuracy of LINCS
  lincs_order             = 4         ; also related to accuracy
  ; Nonbonded settings
  cutoff-scheme           = Verlet    ; Buffered neighbor searching
  ns_type                 = grid      ; search neighboring grid cells
  nstlist                 = 10        ; 20 fs, largely irrelevant with Verlet
  rcoulomb                = 1.0       ; short-range electrostatic cutoff (in nm)
  rvdw                    = 1.0       ; short-range van der Waals cutoff (in nm)
  DispCorr                = EnerPres  ; account for cut-off vdW scheme
  ; Electrostatics
  coulombtype             = PME       ; Particle Mesh Ewald for long-range electrostatics
  pme_order               = 4         ; cubic interpolation
  fourierspacing          = 0.16      ; grid spacing for FFT
  ; Temperature coupling is on
  tcoupl                  = V-rescale             ; modified Berendsen thermostat
  tc-grps                 = Protein Non-Protein   ; two coupling groups - more accurate
  tau_t                   = 0.1     0.1           ; time constant, in ps
  ref_t                   = 300     300           ; reference temperature, one for each group, in K
  ; Pressure coupling is off
  pcoupl                  = no        ; no pressure coupling in NVT
  ; Periodic boundary conditions
  pbc                     = xyz       ; 3-D PBC
  ; Velocity generation
  gen_vel                 = yes       ; assign velocities from Maxwell distribution
  gen_temp                = 300       ; temperature for Maxwell distribution
  gen_seed                = -1        ; generate a random seed

  freezegrps              = Protein
  freezedim               = X Y Z
  """)
  close(file)

  # NPT equilibration
  file=open("npt.mdp","w")
  println(file,"""
  ; title                   =OPLSsystem  NPT equilibration
  define                  = -DPOSRES  ; position restrain the protein
  ; Run parameters
  integrator              = md        ; leap-frog integrator
  nsteps                  = 2500000    ; 2 * 2500000 = 5000 ps
  dt                      = 0.002     ; 2 fs
  ; Output control
  nstxout                 = 0       ; save coordinates every 1.0 ps
  nstvout                 = 0      ; save velocities every 1.0 ps
  nstenergy               = 50000       ; save energies every 1.0 ps
  nstlog                  = 50000       ; update log file every 1.0 ps
  ; Bond parameters
  continuation            = yes       ; Restarting after NVT
  constraint_algorithm    = lincs     ; holonomic constraints
  constraints             = all-bonds   ; bonds involving H are constrained
  lincs_iter              = 1         ; accuracy of LINCS
  lincs_order             = 4         ; also related to accuracy
  ; Nonbonded settings
  cutoff-scheme           = Verlet    ; Buffered neighbor searching
  ns_type                 = grid      ; search neighboring grid cells
  nstlist                 = 10        ; 20 fs, largely irrelevant with Verlet scheme
  rcoulomb                = 1.0       ; short-range electrostatic cutoff (in nm)
  rvdw                    = 1.0       ; short-range van der Waals cutoff (in nm)
  DispCorr                = EnerPres  ; account for cut-off vdW scheme
  ; Electrostatics
  coulombtype             = PME       ; Particle Mesh Ewald for long-range electrostatics
  pme_order               = 4         ; cubic interpolation
  fourierspacing          = 0.16      ; grid spacing for FFT
  ; Temperature coupling is on
  tcoupl                  = V-rescale             ; modified Berendsen thermostat
  tc-grps                 = Protein Non-Protein   ; two coupling groups - more accurate
  tau_t                   = 0.1     0.1           ; time constant, in ps
  ref_t                   = 300     300           ; reference temperature, one for each group, in K
  ; Pressure coupling is on
  pcoupl                  = Parrinello-Rahman     ; Pressure coupling on in NPT
  pcoupltype              = isotropic             ; uniform scaling of box vectors
  tau_p                   = 2.0                   ; time constant, in ps
  ref_p                   = 1.0                   ; reference pressure, in bar
  compressibility         = 4.5e-5                ; isothermal compressibility of water, bar^-1
  refcoord_scaling        = com
  ; Periodic boundary conditions
  pbc                     = xyz       ; 3-D PBC
  ; Velocity generation
  gen_vel                 = no        ; Velocity generation is off """)
  close(file)

  # free simulation
  file=open("md.mdp","w")
  println(file,"""
  ; Run parameters
  integrator              = md        ; leap-frog integrator
  nsteps                  = 500000    ; 2 * 500000 = 1000 ps (1 ns)
  dt                      = 0.002     ; 2 fs
  ; Output control
  nstxout                 = 0         ; suppress bulky .trr file by specifying
  nstvout                 = 0         ; 0 for output frequency of nstxout,
  nstfout                 = 0         ; nstvout, and nstfout
  nstenergy               = 5000      ; save energies every 10.0 ps
  nstlog                  = 5000      ; update log file every 10.0 ps
  nstxout-compressed      = 5000      ; save compressed coordinates every 10.0 ps
  compressed-x-grps       = System    ; save the whole system
  ; Bond parameters
  continuation            = yes       ; Restarting after NPT
  constraint_algorithm    = lincs     ; holonomic constraints
  constraints             = h-bonds   ; bonds involving H are constrained
  lincs_iter              = 1         ; accuracy of LINCS
  lincs_order             = 4         ; also related to accuracy
  ; Neighborsearching
  cutoff-scheme           = Verlet    ; Buffered neighbor searching
  ns_type                 = grid      ; search neighboring grid cells
  nstlist                 = 10        ; 20 fs, largely irrelevant with Verlet scheme
  rcoulomb                = 1.0       ; short-range electrostatic cutoff (in nm)
  rvdw                    = 1.0       ; short-range van der Waals cutoff (in nm)
  ; Electrostatics
  coulombtype             = PME       ; Particle Mesh Ewald for long-range electrostatics
  pme_order               = 4         ; cubic interpolation
  fourierspacing          = 0.16      ; grid spacing for FFT
  ; Temperature coupling is on
  tcoupl                  = V-rescale             ; modified Berendsen thermostat
  tc-grps                 = Protein Non-Protein   ; two coupling groups - more accurate
  tau_t                   = 0.1     0.1           ; time constant, in ps
  ref_t                   = 300     300           ; reference temperature, one for each group, in K
  ; Pressure coupling is on
  pcoupl                  = Parrinello-Rahman     ; Pressure coupling on in NPT
  pcoupltype              = isotropic             ; uniform scaling of box vectors
  tau_p                   = 2.0                   ; time constant, in ps
  ref_p                   = 1.0                   ; reference pressure, in bar
  compressibility         = 4.5e-5                ; isothermal compressibility of water, bar^-1
  ; Periodic boundary conditions
  pbc                     = xyz      ; 3-D PBC
  ; Dispersion correction
  DispCorr                = EnerPres  ; account for cut-off vdW scheme
  ; Velocity generation
  gen_vel                 = no        ; Velocity generation is off """)
  close(file)

  # production
  file=open("md_prod.mdp","w")
  println(file,"""
  ; Run parameters
  integrator              = md        ; leap-frog integrator
  nsteps                  = 5000000    ; 2 *5000000 =10000 ps (10 ns)
  dt                      = 0.002     ; 2 fs
  tinit                   = 1000      ; we have already simulated 1000 ps
  ; Output control
  nstxout                 = 0         ; suppress bulky .trr file by specifying
  nstvout                 = 0         ; 0 for output frequency of nstxout,
  nstfout                 = 0         ; nstvout, and nstfout
  nstenergy               = 5000      ; save energies every 100.0 ps
  nstlog                  = 5000      ; update log file every 100.0 ps
  nstxout-compressed      = 5000      ; save compressed coordinates every 100.0 ps
  compressed-x-grps       = System    ; save the whole system
  ; Bond parameters
  continuation            = yes       ; Restarting after NPT
  constraint_algorithm    = lincs     ; holonomic constraints
  constraints             = h-bonds   ; bonds involving H are constrained
  lincs_iter              = 1         ; accuracy of LINCS
  lincs_order             = 4         ; also related to accuracy
  ; Neighborsearching
  cutoff-scheme           = Verlet    ; Buffered neighbor searching
  ns_type                 = grid      ; search neighboring grid cells
  nstlist                 = 10        ; 20 fs, largely irrelevant with Verlet scheme
  rcoulomb                = 1.0       ; short-range electrostatic cutoff (in nm)
  rvdw                    = 1.0       ; short-range van der Waals cutoff (in nm)
  ; Electrostatics
  coulombtype             = PME       ; Particle Mesh Ewald for long-range electrostatics
  pme_order               = 4         ; cubic interpolation
  fourierspacing          = 0.16      ; grid spacing for FFT
  ; Temperature coupling is on
  tcoupl                  = V-rescale             ; modified Berendsen thermostat
  tc-grps                 = Protein Non-Protein   ; two coupling groups - more accurate
  tau_t                   = 0.1     0.1           ; time constant, in ps
  ref_t                   = 300     300           ; reference temperature, one for each group, in K
  ; Pressure coupling is on
  pcoupl                  = Parrinello-Rahman     ; Pressure coupling on in NPT
  pcoupltype              = isotropic             ; uniform scaling of box vectors
  tau_p                   = 2.0                   ; time constant, in ps
  ref_p                   = 1.0                   ; reference pressure, in bar
  compressibility         = 4.5e-5                ; isothermal compressibility of water, bar^-1
  ; Periodic boundary conditions
  pbc                     = xyz       ; 3-D PBC
  ; Dispersion correction
  DispCorr                = EnerPres  ; account for cut-off vdW scheme
  ; Velocity generation
  gen_vel                 = no        ; Velocity generation is off

  ;energygrps = Protein SOL BMI BF4 NO3""")
  close(file)

end

function mdp_files()
  # minimization
  file = open("mim.mdp", "w")
  println(file,"""
  ;Parameters describing what to do, when to stop and what to save
  integrator  = steep         ; Algorithm (steep = steepest descent minimization)
  emtol       = 100.0        ; Stop minimization when the maximum force < 100.0 kJ/mol/nm
  emstep      = 0.01          ; Minimization step size
  nsteps      = 50000         ; Maximum number of (minimization) steps to perform

  ; Parameters describing how to find the neighbors of each atom and how to calculate the interactions
  nstlist         = 1         ; Frequency to update the neighbor list and long range forces
  cutoff-scheme   = Verlet    ; Buffered neighbor searching
  ns_type         = grid      ; Method to determine neighbor list (simple, grid)
  coulombtype     = PME       ; Treatment of long range electrostatic interactions
  rcoulomb        = 1.0       ; Short-range electrostatic cut-off
  rvdw            = 1.0       ; Short-range Van der Waals cut-off
  pbc             = xyz       ; Periodic Boundary Conditions in all 3 dimensions
  """)
  close(file)

  # NVT equilibration
  file=open("nvt.mdp","w")
  println(file,"""
  ; title                   = OPLS Ubiquitin solvated by water with a ionicliquid being a cosolute NVT equilibration
  ; Run parameters
  integrator              = md        ; leap-frog integrator
  nsteps                  = 500000     ; 2 * 50000 = 1000 ps
  dt                      = 0.002     ; 2 fs
  ; Output control
  nstxout                 = 5000       ; save coordinates every 1.0 ps
  nstvout                 = 5000       ; save velocities every 1.0 ps
  nstenergy               = 5000       ; save energies every 1.0 ps
  nstlog                  = 5000       ; update log file every 1.0 ps
  ; Bond parameters

  continuation            = no        ; first dynamics run
  constraint_algorithm    = lincs     ; holonomic constraints
  constraints             = h-bonds   ; bonds involving H are constrained
  lincs_iter              = 1         ; accuracy of LINCS
  lincs_order             = 4         ; also related to accuracy
  ; Nonbonded settings
  cutoff-scheme           = Verlet    ; Buffered neighbor searching
  ns_type                 = grid      ; search neighboring grid cells
  nstlist                 = 10        ; 20 fs, largely irrelevant with Verlet
  rcoulomb                = 1.0       ; short-range electrostatic cutoff (in nm)
  rvdw                    = 1.0       ; short-range van der Waals cutoff (in nm)
  DispCorr                = EnerPres  ; account for cut-off vdW scheme
  ; Electrostatics
  coulombtype             = PME       ; Particle Mesh Ewald for long-range electrostatics
  pme_order               = 4         ; cubic interpolation
  fourierspacing          = 0.16      ; grid spacing for FFT
  ; Temperature coupling is on
  tcoupl                  = V-rescale             ; modified Berendsen thermostat
  tc-grps                 = System
  tau_t                   = 0.1
  ref_t                   = 300
  ; Pressure coupling is off
  pcoupl                  = no        ; no pressure coupling in NVT
  ; Periodic boundary conditions
  pbc                     = xyz       ; 3-D PBC
  ; Velocity generation
  gen_vel                 = yes       ; assign velocities from Maxwell distribution
  gen_temp                = 300       ; temperature for Maxwell distribution
  gen_seed                = -1        ; generate a random seed """)
  close(file)

  # NPT equilibration
  file=open("npt.mdp","w")
  println(file,"""
  ; Run parameters
  integrator              = md        ; leap-frog integrator
  nsteps                  = 2500000    ; 2 * 2500000 = 5000 ps
  dt                      = 0.002     ; 2 fs
  ; Output control
  nstxout                 = 0       ; save coordinates every 1.0 ps
  nstvout                 = 0      ; save velocities every 1.0 ps
  nstenergy               = 50000       ; save energies every 1.0 ps
  nstlog                  = 50000       ; update log file every 1.0 ps
  ; Bond parameters
  continuation            = yes       ; Restarting after NVT
  constraint_algorithm    = lincs     ; holonomic constraints
  constraints             = all-bonds   ; bonds involving H are constrained
  lincs_iter              = 1         ; accuracy of LINCS
  lincs_order             = 4         ; also related to accuracy
  ; Nonbonded settings
  cutoff-scheme           = Verlet    ; Buffered neighbor searching
  ns_type                 = grid      ; search neighboring grid cells
  nstlist                 = 10        ; 20 fs, largely irrelevant with Verlet scheme
  rcoulomb                = 1.0       ; short-range electrostatic cutoff (in nm)
  rvdw                    = 1.0       ; short-range van der Waals cutoff (in nm)
  DispCorr                = EnerPres  ; account for cut-off vdW scheme
  ; Electrostatics
  coulombtype             = PME       ; Particle Mesh Ewald for long-range electrostatics
  pme_order               = 4         ; cubic interpolation
  fourierspacing          = 0.16      ; grid spacing for FFT
  ; Temperature coupling is on
  tcoupl                  = V-rescale             ; modified Berendsen thermostat
  tc-grps                 = System
  tau_t                   = 0.1
  ref_t                   = 300
  ; Pressure coupling is on
  pcoupl                  = Parrinello-Rahman     ; Pressure coupling on in NPT
  pcoupltype              = isotropic             ; uniform scaling of box vectors
  tau_p                   = 2.0                   ; time constant, in ps
  ref_p                   = 1.0                   ; reference pressure, in bar
  compressibility         = 4.5e-5                ; isothermal compressibility of water, bar^-1
  refcoord_scaling        = com
  ; Periodic boundary conditions
  pbc                     = xyz       ; 3-D PBC
  ; Velocity generation
  gen_vel                 = no        ; Velocity generation is off """)
  close(file)

  # free simulation
  file=open("md.mdp","w")
  println(file,"""
  ; Run parameters
  integrator              = md        ; leap-frog integrator
  nsteps                  = 500000    ; 2 * 500000 = 1000 ps (1 ns)
  dt                      = 0.002     ; 2 fs
  ; Output control
  nstxout                 = 0         ; suppress bulky .trr file by specifying
  nstvout                 = 0         ; 0 for output frequency of nstxout,
  nstfout                 = 0         ; nstvout, and nstfout
  nstenergy               = 5000      ; save energies every 10.0 ps
  nstlog                  = 5000      ; update log file every 10.0 ps
  nstxout-compressed      = 5000      ; save compressed coordinates every 10.0 ps
  compressed-x-grps       = System    ; save the whole system
  ; Bond parameters
  continuation            = yes       ; Restarting after NPT
  constraint_algorithm    = lincs     ; holonomic constraints
  constraints             = h-bonds   ; bonds involving H are constrained
  lincs_iter              = 1         ; accuracy of LINCS
  lincs_order             = 4         ; also related to accuracy
  ; Neighborsearching
  cutoff-scheme           = Verlet    ; Buffered neighbor searching
  ns_type                 = grid      ; search neighboring grid cells
  nstlist                 = 10        ; 20 fs, largely irrelevant with Verlet scheme
  rcoulomb                = 1.0       ; short-range electrostatic cutoff (in nm)
  rvdw                    = 1.0       ; short-range van der Waals cutoff (in nm)
  ; Electrostatics
  coulombtype             = PME       ; Particle Mesh Ewald for long-range electrostatics
  pme_order               = 4         ; cubic interpolation
  fourierspacing          = 0.16      ; grid spacing for FFT
  ; Temperature coupling is on
  tcoupl                  = V-rescale             ; modified Berendsen thermostat
  tc-grps                 = System
  tau_t                   = 0.1
  ref_t                   = 300
  ; Pressure coupling is on
  pcoupl                  = Parrinello-Rahman     ; Pressure coupling on in NPT
  pcoupltype              = isotropic             ; uniform scaling of box vectors
  tau_p                   = 2.0                   ; time constant, in ps
  ref_p                   = 1.0                   ; reference pressure, in bar
  compressibility         = 4.5e-5                ; isothermal compressibility of water, bar^-1
  ; Periodic boundary conditions
  pbc                     = xyz      ; 3-D PBC
  ; Dispersion correction
  DispCorr                = EnerPres  ; account for cut-off vdW scheme
  ; Velocity generation
  gen_vel                 = no        ; Velocity generation is off """)
  close(file)

  # production
  file=open("md_prod.mdp","w")
  println(file,"""
  ; Run parameters
  integrator              = md        ; leap-frog integrator
  nsteps                  = 5000000    ; 2 *5000000 =10000 ps (10 ns)
  dt                      = 0.002     ; 2 fs
  tinit                   = 1000      ; we have already simulated 1000 ps
  ; Output control
  nstxout                 = 0         ; suppress bulky .trr file by specifying
  nstvout                 = 0         ; 0 for output frequency of nstxout,
  nstfout                 = 0         ; nstvout, and nstfout
  nstenergy               = 5000      ; save energies every 100.0 ps
  nstlog                  = 5000      ; update log file every 100.0 ps
  nstxout-compressed      = 5000      ; save compressed coordinates every 100.0 ps
  compressed-x-grps       = System    ; save the whole system
  ; Bond parameters
  continuation            = yes       ; Restarting after NPT
  constraint_algorithm    = lincs     ; holonomic constraints
  constraints             = h-bonds   ; bonds involving H are constrained
  lincs_iter              = 1         ; accuracy of LINCS
  lincs_order             = 4         ; also related to accuracy
  ; Neighborsearching
  cutoff-scheme           = Verlet    ; Buffered neighbor searching
  ns_type                 = grid      ; search neighboring grid cells
  nstlist                 = 10        ; 20 fs, largely irrelevant with Verlet scheme
  rcoulomb                = 1.0       ; short-range electrostatic cutoff (in nm)
  rvdw                    = 1.0       ; short-range van der Waals cutoff (in nm)
  ; Electrostatics
  coulombtype             = PME       ; Particle Mesh Ewald for long-range electrostatics
  pme_order               = 4         ; cubic interpolation
  fourierspacing          = 0.16      ; grid spacing for FFT
  ; Temperature coupling is on
  tcoupl                  = V-rescale             ; modified Berendsen thermostat
  tc-grps                 = System
  tau_t                   = 0.1
  ref_t                   = 300
  ; Pressure coupling is on
  pcoupl                  = Parrinello-Rahman     ; Pressure coupling on in NPT
  pcoupltype              = isotropic             ; uniform scaling of box vectors
  tau_p                   = 2.0                   ; time constant, in ps
  ref_p                   = 1.0                   ; reference pressure, in bar
  compressibility         = 4.5e-5                ; isothermal compressibility of water, bar^-1
  ; Periodic boundary conditions
  pbc                     = xyz       ; 3-D PBC
  ; Dispersion correction
  DispCorr                = EnerPres  ; account for cut-off vdW scheme
  ; Velocity generation
  gen_vel                 = no        ; Velocity generation is off
  ;energygrps = Protein SOL BMI BF4 NO3""")
  close(file)

end

ext: jl


text:
struct PCEulerConstantCache <: StochasticDiffEqConstantCache end

alg_cache(alg::PCEuler,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,::Type{Val{false}}) = PCEulerConstantCache()

@cache struct PCEulerCache{uType,rateType,rateNoiseType,rateNoiseCollectionType} <: StochasticDiffEqMutableCache
  utmp::uType
  ftmp::rateType
  gtmp::rateNoiseType
  gdWtmp::rateNoiseCollectionType
  bbprimetmp::rateType
end

function alg_cache(alg::PCEuler,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,::Type{Val{true}})
  utmp = zero(u); ftmp = zero(rate_prototype);
  gtmp = zero(noise_rate_prototype)
  bbprimetmp = zero(ftmp)
  if is_diagonal_noise(prob)
    gdWtmp = gtmp
  else
    gdWtmp = zero(rate_prototype)
  end
  PCEulerCache(utmp,ftmp,gtmp,gdWtmp,bbprimetmp)
end


text:
# Copyright 2019 Tobias Frilling
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

module DataIO

using DelimitedFiles

include("utils.jl")

include("CLS.jl")
export readCLS

include("LRN.jl")
export LRNData, writeLRN, readLRN

include("UMX.jl")
export readUMX

end # module

max_stars_repo_name: ckafi/DataIO

max_issues_repo_name: ckafi/DataIO

max_forks_repo_name: ckafi/DataIO


text:
using MultistartOptimization
using Test
using NLopt: NLopt

include("test_functions.jl")

@testset "test function sanity checks" begin
    for F in TEST_FUNCTIONS
        @test F(minimum_location(F, 10)) ≈ 0
    end
end

@testset "global optimization" begin
    for F in setdiff(TEST_FUNCTIONS, (RASTRIGIN, )) # Rastrigin disabled for now
        n = 10
        P = MinimizationProblem(F, lower_bounds(F, n), upper_bounds(F, n))
        local_method = NLoptLocalMethod(NLopt.LN_BOBYQA)
        multistart_method = TikTak(100)
        p = multistart_minimization(multistart_method, local_method, P)
        x₀ = minimum_location(F, n)
        @test p.location ≈ x₀ atol = 1e-5
        @test p.value ≈ F(x₀) atol = 1e-10
    end
end

@testset "custom local minimization and infeasibility" begin
    P = MinimizationProblem(x -> sum(abs2, x), [-1, -1], [1, 1])
    x0 = [0.0, 0.0]
    function _local_method(minimization_problem, x)
        if sum(abs2, x) < 0.25
            (location = x0, value = 0.0)
        else
            nothing
        end
    end
    p = multistart_minimization(TikTak(100), _local_method, P)
    @test p.location == x0
    @test p.value == 0.0
end


text:
int_rules_5_5_2 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 5.5.2*Miscellaneous*inverse*secant =#
    @apply_utils Antiderivative(asec(~c + ~(d') * ~x), ~x) => ((~c + ~d * ~x) * asec(~c + ~d * ~x)) / ~d - Antiderivative(1 / ((~c + ~d * ~x) * sqrt(1 - 1 / (~c + ~d * ~x) ^ 2)), ~x) <-- FreeQ([~c, ~d], ~x)
    @apply_utils Antiderivative(acsc(~c + ~(d') * ~x), ~x) => ((~c + ~d * ~x) * acsc(~c + ~d * ~x)) / ~d + Antiderivative(1 / ((~c + ~d * ~x) * sqrt(1 - 1 / (~c + ~d * ~x) ^ 2)), ~x) <-- FreeQ([~c, ~d], ~x)
    @apply_utils Antiderivative((~(a') + ~(b') * asec(~c + ~(d') * ~x)) ^ ~(p'), ~x) => (1 / ~d) * Subst(Antiderivative((~a + ~b * asec(~x)) ^ ~p, ~x), ~x, ~c + ~d * ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && IGtQ(~p, 0)
    @apply_utils Antiderivative((~(a') + ~(b') * acsc(~c + ~(d') * ~x)) ^ ~(p'), ~x) => (1 / ~d) * Subst(Antiderivative((~a + ~b * acsc(~x)) ^ ~p, ~x), ~x, ~c + ~d * ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && IGtQ(~p, 0)
    @apply_utils Antiderivative((~(a') + ~(b') * asec(~c + ~(d') * ~x)) ^ ~p, ~x) => Unintegrable((~a + ~b * asec(~c + ~d * ~x)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~p], ~x) && Not(IGtQ(~p, 0))
    @apply_utils Antiderivative((~(a') + ~(b') * acsc(~c + ~(d') * ~x)) ^ ~p, ~x) => Unintegrable((~a + ~b * acsc(~c + ~d * ~x)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~p], ~x) && Not(IGtQ(~p, 0))
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * (~(a') + ~(b') * asec(~c + ~(d') * ~x)) ^ ~(p'), ~x) => (1 / ~d) * Subst(Antiderivative(((~f * ~x) / ~d) ^ ~m * (~a + ~b * asec(~x)) ^ ~p, ~x), ~x, ~c + ~d * ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m], ~x) && (EqQ(~d * ~e - ~c * ~f, 0) && IGtQ(~p, 0))
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * (~(a') + ~(b') * acsc(~c + ~(d') * ~x)) ^ ~(p'), ~x) => (1 / ~d) * Subst(Antiderivative(((~f * ~x) / ~d) ^ ~m * (~a + ~b * acsc(~x)) ^ ~p, ~x), ~x, ~c + ~d * ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m], ~x) && (EqQ(~d * ~e - ~c * ~f, 0) && IGtQ(~p, 0))
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * (~(a') + ~(b') * asec(~c + ~(d') * ~x)) ^ ~(p'), ~x) => (1 / (~d) ^ (~m + 1)) * Subst(Antiderivative((~a + ~b * ~x) ^ ~p * sec(~x) * tan(~x) * ((~d * ~e - ~c * ~f) + ~f * sec(~x)) ^ ~m, ~x), ~x, asec(~c + ~d * ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (IGtQ(~p, 0) && IntegerQ(~m))
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * (~(a') + ~(b') * acsc(~c + ~(d') * ~x)) ^ ~(p'), ~x) => (-1 / (~d) ^ (~m + 1)) * Subst(Antiderivative((~a + ~b * ~x) ^ ~p * csc(~x) * cot(~x) * ((~d * ~e - ~c * ~f) + ~f * csc(~x)) ^ ~m, ~x), ~x, acsc(~c + ~d * ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (IGtQ(~p, 0) && IntegerQ(~m))
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * (~(a') + ~(b') * asec(~c + ~(d') * ~x)) ^ ~(p'), ~x) => (1 / ~d) * Subst(Antiderivative(((~d * ~e - ~c * ~f) / ~d + (~f * ~x) / ~d) ^ ~m * (~a + ~b * asec(~x)) ^ ~p, ~x), ~x, ~c + ~d * ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m], ~x) && IGtQ(~p, 0)
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * (~(a') + ~(b') * acsc(~c + ~(d') * ~x)) ^ ~(p'), ~x) => (1 / ~d) * Subst(Antiderivative(((~d * ~e - ~c * ~f) / ~d + (~f * ~x) / ~d) ^ ~m * (~a + ~b * acsc(~x)) ^ ~p, ~x), ~x, ~c + ~d * ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m], ~x) && IGtQ(~p, 0)
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * (~(a') + ~(b') * asec(~c + ~(d') * ~x)) ^ ~p, ~x) => Unintegrable((~e + ~f * ~x) ^ ~m * (~a + ~b * asec(~c + ~d * ~x)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m, ~p], ~x) && Not(IGtQ(~p, 0))
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * (~(a') + ~(b') * acsc(~c + ~(d') * ~x)) ^ ~p, ~x) => Unintegrable((~e + ~f * ~x) ^ ~m * (~a + ~b * acsc(~c + ~d * ~x)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m, ~p], ~x) && Not(IGtQ(~p, 0))
    @apply_utils Antiderivative(~(u') * asec(~(c') / (~(a') + ~(b') * (~x) ^ ~(n'))) ^ ~(m'), ~x) => Antiderivative(~u * acos(~a / ~c + (~b * (~x) ^ ~n) / ~c) ^ ~m, ~x) <-- FreeQ([~a, ~b, ~c, ~n, ~m], ~x)
    @apply_utils Antiderivative(~(u') * acsc(~(c') / (~(a') + ~(b') * (~x) ^ ~(n'))) ^ ~(m'), ~x) => Antiderivative(~u * asin(~a / ~c + (~b * (~x) ^ ~n) / ~c) ^ ~m, ~x) <-- FreeQ([~a, ~b, ~c, ~n, ~m], ~x)
    @apply_utils Antiderivative(~(u') * (~f) ^ (~(c') * asec(~(a') + ~(b') * ~x) ^ ~(n')), ~x) => (1 / ~b) * Subst(Antiderivative(ReplaceAll(~u, (~x->begin
                                        #= none:1 =#
                                        #= none:3 =#
                                        -(~a) / ~b + sec(~x) / ~b
                                    end)) * (~f) ^ (~c * (~x) ^ ~n) * sec(~x) * tan(~x), ~x), ~x, asec(~a + ~b * ~x)) <-- FreeQ([~a, ~b, ~c, ~f], ~x) && IGtQ(~n, 0)
    @apply_utils Antiderivative(~(u') * (~f) ^ (~(c') * acsc(~(a') + ~(b') * ~x) ^ ~(n')), ~x) => (-1 / ~b) * Subst(Antiderivative(ReplaceAll(~u, (~x->begin
                                        #= none:1 =#
                                        #= none:3 =#
                                        -(~a) / ~b + csc(~x) / ~b
                                    end)) * (~f) ^ (~c * (~x) ^ ~n) * csc(~x) * cot(~x), ~x), ~x, acsc(~a + ~b * ~x)) <-- FreeQ([~a, ~b, ~c, ~f], ~x) && IGtQ(~n, 0)
    @apply_utils Antiderivative(asec(~u), ~x) => ~x * asec(~u) - (~u / sqrt((~u) ^ 2)) * Antiderivative(SimplifyIntegrand((~x * D(~u, ~x)) / (~u * sqrt((~u) ^ 2 - 1)), ~x), ~x) <-- InverseFunctionFreeQ(~u, ~x) && Not(FunctionOfExponentialQ(~u, ~x))
    @apply_utils Antiderivative(acsc(~u), ~x) => ~x * acsc(~u) + (~u / sqrt((~u) ^ 2)) * Antiderivative(SimplifyIntegrand((~x * D(~u, ~x)) / (~u * sqrt((~u) ^ 2 - 1)), ~x), ~x) <-- InverseFunctionFreeQ(~u, ~x) && Not(FunctionOfExponentialQ(~u, ~x))
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * (~(a') + ~(b') * asec(~u)), ~x) => ((~c + ~d * ~x) ^ (~m + 1) * (~a + ~b * asec(~u))) / (~d * (~m + 1)) - ((~b * ~u) / (~d * (~m + 1) * sqrt((~u) ^ 2))) * Antiderivative(SimplifyIntegrand(((~c + ~d * ~x) ^ (~m + 1) * D(~u, ~x)) / (~u * sqrt((~u) ^ 2 - 1)), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m], ~x) && (NeQ(~m, -1) && (InverseFunctionFreeQ(~u, ~x) && (Not(FunctionOfQ((~c + ~d * ~x) ^ (~m + 1), ~u, ~x)) && Not(FunctionOfExponentialQ(~u, ~x)))))
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * (~(a') + ~(b') * acsc(~u)), ~x) => ((~c + ~d * ~x) ^ (~m + 1) * (~a + ~b * acsc(~u))) / (~d * (~m + 1)) + ((~b * ~u) / (~d * (~m + 1) * sqrt((~u) ^ 2))) * Antiderivative(SimplifyIntegrand(((~c + ~d * ~x) ^ (~m + 1) * D(~u, ~x)) / (~u * sqrt((~u) ^ 2 - 1)), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m], ~x) && (NeQ(~m, -1) && (InverseFunctionFreeQ(~u, ~x) && (Not(FunctionOfQ((~c + ~d * ~x) ^ (~m + 1), ~u, ~x)) && Not(FunctionOfExponentialQ(~u, ~x)))))
    @apply_utils Antiderivative(~v * (~(a') + ~(b') * asec(~u)), ~x) => With([w = IntHide(~v, ~x)], Dist(~a + ~b * asec(~u), w, ~x) - ((~b * ~u) / sqrt((~u) ^ 2)) * Antiderivative(SimplifyIntegrand((w * D(~u, ~x)) / (~u * sqrt((~u) ^ 2 - 1)), ~x), ~x) <-- InverseFunctionFreeQ(w, ~x)) <-- FreeQ([~a, ~b], ~x) && (InverseFunctionFreeQ(~u, ~x) && Not(MatchQ(~v, (~(c') + ~(d') * ~x) ^ ~(m') <-- FreeQ([c, d, m], ~x))))
    @apply_utils Antiderivative(~v * (~(a') + ~(b') * acsc(~u)), ~x) => With([w = IntHide(~v, ~x)], Dist(~a + ~b * acsc(~u), w, ~x) + ((~b * ~u) / sqrt((~u) ^ 2)) * Antiderivative(SimplifyIntegrand((w * D(~u, ~x)) / (~u * sqrt((~u) ^ 2 - 1)), ~x), ~x) <-- InverseFunctionFreeQ(w, ~x)) <-- FreeQ([~a, ~b], ~x) && (InverseFunctionFreeQ(~u, ~x) && Not(MatchQ(~v, (~(c') + ~(d') * ~x) ^ ~(m') <-- FreeQ([c, d, m], ~x))))
end

max_stars_repo_path: src/rules/5 Inverse trig functions/5.5 Inverse secant/5.5.2 Miscellaneous inverse secant.jl


text:
module Romeo

using GLWindow
using GLAbstraction
using ModernGL
using ImmutableArrays
using Reactive
using GLFW
using Images
using Quaternions
using GLText
using Compat
using Color
using FixedPointNumbers
import Mustache

const sourcedir = Pkg.dir("Romeo", "src")
const shaderdir = joinpath(sourcedir, "shader")

function maxper(v0::Vector3, v1::Vector3)
	return Vector3(max(v0[1], v1[1]),
            max(v0[2], v1[2]),
            max(v0[3], v1[3]))
end
function minper(v0::Vector3, v1::Vector3)
	return Vector3(min(v0[1], v1[1]),
            min(v0[2], v1[2]),
            min(v0[3], v1[3]))
end

Base.minimum{T, NDIM}(x::Array{Vector3{T},NDIM}) = reduce(minper, x)
Base.maximum{T, NDIM}(x::Array{Vector3{T},NDIM}) = reduce(maxper, x)

include(joinpath(     sourcedir, "utils.jl"))
include(joinpath(     sourcedir, "types.jl"))
include_all(joinpath( sourcedir, "display"))
include(joinpath(     sourcedir, "color.jl"))
include_all(joinpath( sourcedir, "share"))
include_all(joinpath( sourcedir, "edit"))
include_all(joinpath( sourcedir, "visualize"))
include(joinpath(     sourcedir, "visualize_interface.jl"))
include(joinpath(     sourcedir, "edit_interface.jl"))

export visualize    # Visualize an object
export edit         # Edit an object

export RGBAU8       # typealias for RGBA ufixed 8 value
export rgba         # function for creating a rgba Float32 color
export rgbaU8       # function for creating a rgba Ufixed8 color

export tohsva       # Convert to HSVA
export torgba       # Converts to RGBA

# Surface Rendering
export mix      # mix colors
export SURFACE  # function that generates a Surface primitive for every datapoint, with an optional gap between the surfaces
export CIRCLE   # function that generates Circular surface primitive for every datapoint
export CUBE     # function that generates Cube primitives for every datapoint
export POINT    # function that generates Point primitives for every datapoint

end # module


text:
"""$(TYPEDEF)
iLQR is an unconstrained indirect method for trajectory optimization that parameterizes only the controls and enforces strict dynamics feasibility at every iteration by simulating forward the dynamics with an LQR feedback controller.
The main algorithm consists of two parts:
1) a backward pass that uses Differential Dynamic Programming to compute recursively a quadratic approximation of the cost-to-go, along with linear feedback and feed-forward gain matrices, `K` and `d`, respectively, for an LQR tracking controller, and
2) a forward pass that uses the gains `K` and `d` to simulate forward the full nonlinear dynamics with feedback.
"""
struct StaticiLQRSolver{T,I<:QuadratureRule,L,O,n,m,L1,D,F,E1,E2,A} <: iLQRSolver{T}
    # Model + Objective
    model::L
    obj::O

    # Problem info
    x0::MVector{n,T}
    xf::SVector{n,T}
    tf::T
    N::Int

    opts::iLQRSolverOptions{T}
    stats::iLQRStats{T}

    # Primal Duals
    Z::Vector{KnotPoint{T,n,m,L1}}
    Z̄::Vector{KnotPoint{T,n,m,L1}}

    # Data variables
    # K::Vector{SMatrix{m,n̄,T,L2}}  # State feedback gains (m,n,N-1)
    K::Vector{A}  # State feedback gains (m,n,N-1)
    d::Vector{SVector{m,T}} # Feedforward gains (m,N-1)

    ∇F::Vector{D} # discrete dynamics jacobian (block) (n,n+m+1,N)
    G::Vector{F}  # state difference jacobian (n̄, n)

    S::E1  # Optimal cost-to-go expansion trajectory
    Q::E2  # cost-to-go expansion trajectory

    ρ::Vector{T} # Regularization
    dρ::Vector{T} # Regularization rate of change

    grad::Vector{T} # Gradient

    logger::SolverLogger

    function StaticiLQRSolver{T,I}(model::L, obj::O, x0, xf, tf, N, opts, stats,
            Z::Vector{KnotPoint{T,n,m,L1}}, Z̄, K::Vector{A}, d,
            ∇F::Vector{D}, G::Vector{F}, S::E1, Q::E2, ρ, dρ, grad,
            logger) where {T,I,L,O,n,m,L1,D,F,E1,E2,A}
        new{T,I,L,O,n,m,L1,D,F,E1,E2,A}(model, obj, x0, xf, tf, N, opts, stats, Z, Z̄, K, d,
            ∇F, G, S, Q, ρ, dρ, grad, logger)
    end
end

function StaticiLQRSolver(prob::Problem{I,T}, opts=SolverOptions{T}()) where {I,T}

    # Init solver statistics
    stats = iLQRStats{T}() # = Dict{Symbol,Any}(:timer=>TimerOutput())

    # Init solver results
    n,m,N = size(prob)
    n̄ = state_diff_size(prob.model)

    x0 = SVector{n}(prob.x0)
    xf = SVector{n}(prob.xf)

    Z = prob.Z
    # Z̄ = Traj(n,m,Z[1].dt,N)
    Z̄ = copy(prob.Z)

    if m*n̄ > MAX_ELEM
		K  = [zeros(T,m,n̄) for k = 1:N-1]
	else
		K  = [@SMatrix zeros(T,m,n̄) for k = 1:N-1]
	end
    d  = [@SVector zeros(T,m)   for k = 1:N-1]

	if n*(n+m+1) > MAX_ELEM
		∇F = [zeros(T,n,n+m+1) for k = 1:N-1]
	else
		∇F = [@SMatrix zeros(T,n,n+m+1) for k = 1:N-1]
	end
    ∇F = [@SMatrix zeros(T,n,n+m+1) for k = 1:N-1]
    G = [state_diff_jacobian(prob.model, x0) for k = 1:N]

    S = CostExpansion(n̄,m,N)
    Q = CostExpansion(n̄,m,N)

    ρ = zeros(T,1)
    dρ = zeros(T,1)

    grad = zeros(T,N-1)

    logger = default_logger(opts.verbose)

	opts_ilqr = iLQRSolverOptions(opts)

    solver = StaticiLQRSolver{T,I}(prob.model, prob.obj, x0, xf, prob.tf, N, opts_ilqr, stats,
        Z, Z̄, K, d, ∇F, G, S, Q, ρ, dρ, grad, logger)

    reset!(solver)
    return solver
end


text:
Pkg.update()

# Install required Julia packages
Pkg.add("SpecialFunctions")
Pkg.add("FastGaussQuadrature")
Pkg.add("GSL")
Pkg.add("IterativeSolvers")
Pkg.add("LinearMaps")
Pkg.add("PyPlot")
Pkg.add("SymPy")
Pkg.add("PyCall")
Pkg.add("LowRankApprox")
Pkg.add("MATLAB")
Pkg.add("NullableArrays")
Pkg.add("Revise")
Pkg.add("FMMLIB2D")
Pkg.add("FINUFFT")
Pkg.add("LaTeXStrings")
Pkg.add("JLD")
Pkg.add("Glob")
Pkg.add("HDF5")
Pkg.add("FileIO")
Pkg.add("Humanize")
Pkg.add("Interpolations")


text:
# This file contains code for native implementations of some of the
# benchmarks in NLEIGS

function nlevp_native_gun()
    gunbase=joinpath(dirname(@__FILE__()),
      "converted_nlevp", "gun_")
    K=read_sparse_matrix(gunbase * "K.txt", Int64)
    M=read_sparse_matrix(gunbase * "M.txt", Int64)
    W1=read_sparse_matrix(gunbase * "W1.txt", Int64)
    W2=read_sparse_matrix(gunbase * "W2.txt", Int64)
    # The gun problem is a sum of a PEP and a problem containing square roots.
    pep=PEP([K,-M])
    sqrt1op= S -> 1im*sqrt(S)
    sqrt2op= S -> 1im*sqrt(S-108.8774^2*one(S))
    sqrtnep=SPMF_NEP([W1,W2],[sqrt1op,sqrt2op])
    nep=SumNEP(pep,sqrtnep)
    return nep
end

function nlevp_native_cd_player()
    cdbase=joinpath(dirname(@__FILE__()),
      "converted_nlevp", "cd_player_")
    K = Matrix(read_sparse_matrix(cdbase * "K.txt", Int64))
    C = Matrix(read_sparse_matrix(cdbase * "C.txt", Int64))
    M = one(K)
    nep = PEP([K, C, M])
    return nep
end

function nlevp_native_fiber()
    # Since the bessel functions are not available as matrix functions
    # we rely on interpolation (of "denominator" and "numerator" separately)
     # Construct the complicated third function
    L::Float64=2400;
    besselkp= (m,z)->  - besselk(m-1,z) - m*besselk(m,z)/z;  # Derivative
    numer= x-> ((L+0.5)/L^2)*x/(besselk(1,ComplexF64(x))^2)
    denom= x-> 1/(besselkp(1,ComplexF64(x))*besselk(1,ComplexF64(x)));
    f_sqrt= x ->  numer(x)/denom(x)
    s3=λ -> f_sqrt(sqrt(λ)*L)  # This is the complicated function in "fiber"

    m=10;
    TT=Complex{BigFloat}; # Do interpolation very accurately
    ### Different interpolation approaches ####
    ## Disc
    #phi=2im*pi*range(0.0;length=m+1,stop=1);
    #interp_points=0.1*exp.(phi[1:end-1]); #
    ## Uniform on an interval
    interp_points=0.01.+range(0;length=m,stop=3);; # This choice is done with eigenvalue information in mind (x approx 2.027892679678583)
    interp_points=Vector{TT}(interp_points) # High precision

    # Do the interpolation on numerator and denominator
    (Newton_Matrix,fnum)=construct_newton_matrix(TT,numer,interp_points)
    (Newton_Matrix,fdenom)=construct_newton_matrix(TT,denom,interp_points)
    num_coeffs=Newton_Matrix\fnum;
    denom_coeffs=Newton_Matrix\fdenom;

    # Now go back to ComplexF64
    interp_points_c64=Vector{ComplexF64}(interp_points);
    num_coeffs_c64=Vector{ComplexF64}(num_coeffs);
    denom_coeffs_c64=Vector{ComplexF64}(denom_coeffs);

    # Interpolated numerator and denominators
    numer_new_c64=x->newton_eval(num_coeffs_c64,x,interp_points_c64)
    denom_new_c64=x->newton_eval(denom_coeffs_c64,x,interp_points_c64)

    f_new_c64=x->  denom_new_c64(x)\numer_new_c64(x)
    s3_new= λ -> f_new_c64(sqrt(λ)*L); # This is the new function. Works for matrices and functions
    n=2400;

    A2=sparse([n],[n],[1.0])  # This is a rank-one matrix!
    A1=one(A2)

    # Create the A0-matrix takes a bit more work
    eta_cl = 1.4969;  # "physical params"
    alpha=25; ell=1.1;
    gam=0.003; delta=0.01;

    k_cl = 2*pi*eta_cl/ell;
    n_c = 400; n = 6*n_c;
    r = (1:n+1)*delta;
    mm = 1;

    inc = (1:n_c);
    i_n = (n_c+1:n-1);
    e = ones(n_c);

    # Helper functions. Note C(r) is r-independent, according to NLEVP.
    C = sqrt.( (1 .- 2*gam*(inc/n_c).^alpha) / (1 - 2*gam) ) .- 1;
    eta0 = r-> (eta_cl .+ 1.4201*C)
    k = r -> 2*pi*eta0(r)/ell;

    # setup the vectors in the diagonal
    y1 = -2*e - mm^2*(e ./ inc.^2) + delta^2*(k(r[1:n_c]).^2 .- k_cl^2);
    e = ones(size(i_n));
    y2 = -2*e - mm^2*(e ./ i_n.^2);
    y = [y1; y2; -1 + 1/(2*n) - mm^2/n^2];
    i = 1:n-1;
    z = (i.+0.5) ./ sqrt.( i.*(i.+1) );

    A0=spdiagm(0  => y[1:n], 1=> z[1:n-1], -1 => z[1:n-1])

    # The functions
    f1=S-> one(S);
    f2=S-> -S;
    f3=s3_new; # Interpolated function
    # f3=s3 # Use this if you want the exact formula

    nep=SPMF_NEP([A0,A1,A2],[f1,f2,f3])

    return nep;
end

# Helper functions for Newton interpolation (mainly for nlevp_native_fiber)
function construct_newton_matrix(T,ff,interp_points)
    m=size(interp_points,1);
    Newton_Matrix=zeros(T,m,m);
    Newton_Matrix[:,1] .= 1
    for col=2:m
        for row=col:m
            Newton_Matrix[row,col]=Newton_Matrix[row,col-1]*(interp_points[row]-interp_points[col-1])
        end
    end
    f=zeros(T,m);
    for k=1:m
        f[k]=ff(interp_points[k]);
    end
    return Newton_Matrix,f
end
function newton_eval(coeffs,S,interp_points)  # This works for λ::Number and λ::Matrix
    F=coeffs[1]*one(S);
    prod=one(S);
    for k=2:size(coeffs,1)
        prod = prod*(S-interp_points[k-1]*one(S))
        F += prod*coeffs[k];
    end
    return F
end

function nlevp_native_hadeler(α=100,n=8)
    i=1:n
    I2 = ones(n)*i'  # matrix with constant columns 1,2,3,4...n
    # Identity matrix
    II = Matrix{typeof(α)}(I,n,n);
    # Matrices
    A0= α*II
    A2 = n*II+1 ./(I2 + I2');
    B = ((n+1) .- max.(I2',I2)) .* (i*i');
    # Functions
    fv= [S->-one(S) ; S->S^2 ; S->(exp(S)-one(S))]
    nep=SPMF_NEP([A0, A2, B], fv);
    return nep
end

function nlevp_native_pdde_stability(n=15)

    # default values
    a0 = 2;
    b0 = .3;
    a1 = -2;
    b1 = .2;
    a2 = -2;
    b2 = -.3;
    t1 = -pi/2;

    # Construct a DEP
    h = pi/(n+1);  x = (1:n)*h;
    e = ones(n);
    A0 = spdiagm(-1 => e[1:end-1], 0 => -2*e, 1 => e[1:end-1])/(h^2)
    A0 = A0 + spdiagm(0 => (a0 .+ b0*sin.(x)))
    A1 = spdiagm(0 => (a1 .+ b1*x.*(1 .- exp.(x .- pi))))
    A2 = spdiagm(0 => (a2 .+ b2*x.*(pi .- x)))

    # The stabiltiy of the DEP given by -λI+A0+A1exp(-t1 λ)+A2exp(- TT λ) can be
    # as a function of TT can be analyzed with a QEP
    II=complex(sparse(1.0*I,n,n))
    E = kron(II,A2);
    gamma = exp(1im*t1); gamma = gamma/abs(gamma);
    F = kron(II,A0-gamma*A1) + kron(A0+gamma*A1,II);

    # Compute a permutation vector (faster than multiplication)
    p = vec(copy(reshape(1:n^2,n,n)'));

    Av = [E[p,p],F,E];

    return PEP(Av)
end

function toeplitz(v::AbstractVector)
    # This can be replaced with the types in ToeplitzMatrices
    # once the package is more stable (e.g. matrix matrix products works in general).
    n=length(v);
    T=zeros(eltype(v),n,n);
    for i=1:n
        for j=1:(n-i+1)
            T[i,j+i-1]=v[j]
            T[j+i-1,i]=v[j]
        end
    end
    if (issparse(v))
        return SparseMatrixCSC(T)
    else
        return T;
    end

end

function nlevp_native_loaded_string(n=20,kappa=1,m=1)
    # See Gallery for documentation
    A0 = toeplitz(SparseVector([2.0*n;-n;zeros(n-2)]));
    A1 = zeros(n,n); A1[n,n] = n-A0[n,n];
    B0 = toeplitz(SparseVector([4/(6*n);1/(6*n);zeros(n-2)]));
    B1 = zeros(n,n); B1[n,n] = 2/(6*n)-B0[n,n]
    C = zeros(n,n); C[n,n] = kappa;
    σ=kappa/m;
    fv=[S->one(S), S->-S, S -> (S-σ*one(S))\S]

    spmf1=SPMF_NEP([A0,B0],fv[1:2]);
    spmf2=SPMF_NEP([A1,B1,C],fv);  # This can be a Low-Rank SPMF

    nep=SPMFSumNEP(spmf1,spmf2)
    return nep;

end


text:
#=
testODE:
- Julia version:
- Author: ymocquar
- Date: 2019-11-13
=#
include("../src/interface.jl")
using LinearAlgebra
using Plots
using Plots.PlotMeasures
function fctMain(n_tau)

    #   u0 =[big"0.12345678", big"0.1209182736", big"0.1290582671", big"0.1239681094" ]
    tab_eps = zeros(BigFloat, 7)
    epsilon = big"0.15"
    for i = 1:size(tab_eps, 1)
        tab_eps[i] = epsilon
        epsilon /= 10
    end
    u0 = BigFloat.([-34 // 100, 78 // 100, 67 // 100, -56 // 10])
    B =
        BigFloat.([
            12//100 -78//100 91//100 34//100
            -45//100 56//100 3//100 54//100
            -67//100 09//100 18//100 89//100
            -91//100 -56//100 11//100 -56//100
        ])
    alpha = BigFloat.([12 // 100, -98 // 100, 45 // 100, 26 // 100])
    beta = BigFloat.([-4 // 100, 48 // 100, 23 // 100, -87 // 100])
    fct = (u, p, t) -> B * u + t * p[1] + p[2]
    t_max = big"1.0"
    nbMaxTest = 12
    order = 6
    y = ones(Float64, nbMaxTest, size(tab_eps, 1))
    x = zeros(Float64, nbMaxTest)
    ind = 1
    A = [0 0 1 0; 0 0 0 0; -1 0 0 0; 0 0 0 0]
    for epsilon in tab_eps
        prob = HiOscODEProblem(fct, u0, (big"0.0", t_max), (alpha, beta), A, epsilon, B)
        eps_v = convert(Float32, epsilon)
        println("epsilon = $eps_v")
        nb = 100
        indc = 1
        labels = Array{String,2}(undef, 1, ind)
        while indc <= nbMaxTest
            res = solve(prob, nb_tau = n_tau, order = order, nb_t = nb, dense = false)
            sol = res[end]
            solref = getexactsol(res.par_u0.parphi, u0, t_max)
            println("solref=$solref")
            println("nb=$nb sol=$sol")
            diff = solref - sol
            x[indc] = t_max / nb
            println("nb=$nb dt=$(1.0/nb) normInf=$(norm(diff,Inf)) norm2=$(norm(diff))")
            ni = norm(diff, Inf)
            y[indc, ind] = (ni < 1) ? ni : NaN
            println("epsilon=$epsilon result=$y")
            println("epsilon=$epsilon reslog2=$(log2.(y))")
            nb *= 2
            indc += 1
        end
        for i = 1:ind
            labels[1, i] = " ε = $(convert(Float32,tab_eps[i])) "
        end
        p = Plots.plot(
            x,
            view(y, :, 1:ind),
            xlabel = "Δt",
            xaxis = :log,
            ylabel = "error",
            yaxis = :log,
            legend = :bottomright,
            label = labels,
            marker = 2,
        )
        pp = Plots.plot(
            x,
            view(y, :, 1:ind),
            xlabel = "Δt",
            xaxis = :log,
            ylabel = "error",
            yaxis = :log,
            legend = :bottomright,
            label = labels,
            marker = 2,
            bottom_margin = 30px,
            left_margin = 60px,
        )
        prec_v = precision(BigFloat)
        Plots.savefig(p, "out/r10p_$(prec_v)_$(eps_v)_$(order)_$(n_tau)_epsilon.pdf")
        Plots.savefig(p, "out/r10p_$(prec_v)_$(eps_v)_$(order)_$(n_tau)_epsilon.png")
        Plots.savefig(pp, "out/r10pp_$(prec_v)_$(eps_v)_$(order)_$(n_tau)_epsilon.pdf")
        Plots.savefig(pp, "out/r10pp_$(prec_v)_$(eps_v)_$(order)_$(n_tau)_epsilon.png")
        ind += 1
    end
end

# testODESolverEps()

# for i=3:9
#     fctMain(2^i)
# end
# setprecision(512)
fctMain(32)


text:
struct RockPaperScissors end

n_agents(simpleGame::RockPaperScissors) = 2

ordered_actions(simpleGame::RockPaperScissors, i::Int) = [:rock, :paper, :scissors]
ordered_joint_actions(simpleGame::RockPaperScissors) = vec(collect(Iterators.product([ordered_actions(simpleGame, i) for i in 1:n_agents(simpleGame)]...)))

n_joint_actions(simpleGame::RockPaperScissors) = length(ordered_joint_actions(simpleGame))
n_actions(simpleGame::RockPaperScissors, i::Int) = length(ordered_actions(simpleGame, i))

function reward(simpleGame::RockPaperScissors, i::Int, a)
    if i == 1
        noti = 2
    else
        noti = 1
    end

    if a[i] == a[noti]
        r = 0.0
    elseif a[i] == :rock && a[noti] == :paper
        r = -1.0
    elseif a[i] == :rock && a[noti] == :scissors
        r = 1.0
    elseif a[i] == :paper && a[noti] == :rock
        r = 1.0
    elseif a[i] == :paper && a[noti] == :scissors
        r = -1.0
    elseif a[i] == :scissors && a[noti] == :rock
        r = -1.0
    elseif a[i] == :scissors && a[noti] == :paper
        r = 1.0
    end

    return r
end

function joint_reward(simpleGame::RockPaperScissors, a)
    return [reward(simpleGame, i, a) for i in 1:n_agents(simpleGame)]
end

function SimpleGame(simpleGame::RockPaperScissors)
    return SimpleGame(
        0.9,
        vec(collect(1:n_agents(simpleGame))),
        [ordered_actions(simpleGame, i) for i in 1:n_agents(simpleGame)],
        (a) -> joint_reward(simpleGame, a)
    )
end


text:
@testset "obl06" begin
    obl = obl06(2400000.5, 54388.0)
    @test obl ≈ 0.4090749229387258204 atol=1e-16
end

@testset "obl80" begin
    obl = obl80(2400000.5, 54388.0)
    @test obl ≈ 0.409075134764381621 atol=1e-16
end


text:
using Finch, BenchmarkTools

SUITE = BenchmarkGroup()

include("rawcode.jl")

max_stars_repo_name: paralab/Finch

max_issues_repo_name: paralab/Finch

max_forks_repo_name: paralab/Finch


text:
using LinearAlgebra, SparseArrays, DifferentialEquations

#=
    Scenario 1: Nominal actuators
    Scenario 2: Realistic actuators
=#
const scenario = 1;

# Configuration of physical properties λ and c
function property(temperature :: T , a₀ :: T, a₁ :: T) where T <: Real
    Θ = temperature

    return a₀ + a₁ * Θ
end

# Emitted heat flux on boundary
function boundary_heatflux(boundary_temperature :: T, ambient_temperature ::T, transfer_coefficient :: T, radiation_coefficient :: T) where T <: Real
    θbound = boundary_temperature;
    θamb = ambient_temperature
    h = transfer_coefficient
    k = radiation_coefficient

    return -1*h *(θbound - θamb) -  k *(θbound^4 - θamb^4)
end

# Emitted heat flux on boundary
function boundary_heatflux(boundary_temperature :: Array{T}, ambient_temperature :: Array{T}, transfer_coefficient :: T, radiation_coefficient :: T) where T <: Real
    θbound = boundary_temperature;
    θamb = ambient_temperature

    if size(θbound) != size(θamb)
        error("Boundary temperature array and ambient temperature array have different sizes: size(boundary_temperature) = $(size(θbound)) ≠ $(size(θamb)) = size(ambient_temperature)!")
    end

    h = transfer_coefficient
    k = radiation_coefficient

    return -1*h *(θbound .- θamb) -  k *(θbound.^4 .- θamb.^4)
end

# Spatial approximation
const L = 300.0*10^(-3) # Length of rod
const H = 10.0*10^(-3)  # Height of rod

const N₁ = 100; # Number of cells in x₁ direction
const N₂ = 40;  # Number of cells in x₂ direction
const Δx₁ = L/N₁; # Discretizaion
const Δx₂ = H/N₂;
const x1span = Δx₁/2 : Δx₁ : L; # Coordinates of cell centers
const x2span = Δx₂/2 : Δx₂ : H;

# Heat transfer + radiation
const h = 10.0             # Heat transfer coefficient
const ϵ = 0.6;             # Emissivity
const sb = 5.67*10^(-8);   # Stefan-Boltzmann constant
const k = ϵ*sb;            # Heat radiation coefficient

const ρ = 7800.0;       # Density
const ρinv = inv(ρ);
const Id12 = sparse(I, N₁*N₂, N₁*N₂) # Template matrix for specific heat capacity

const θamb = 300.0; # Ambient temperature

# Initial noisy data
const a₀ = 3;   # Scaling
const a₁ = 10;
const a₂ = 5;
const noise = a₀ * cos.(a₁*2*pi/L * x1span) .* transpose(cos.(a₂ *2*pi/H * x2span))
θinit = ones(N₁, N₂) * 300.0 .+ noise;
θinit = reshape(θinit, N₁*N₂, 1 )

# Input / Output
const Nu = 5; # Number of inputs
const Ny = 5; # Number of outputs

b_in = spzeros(N₁, Nu)          # Actuators' spatial characterization
g_out = spzeros(N₁, Ny)         # Sensors' spatial characterization

const β_in = spzeros(Int64,N₁, Nu)    # Input partition
const γ_out = spzeros(Int64, N₁, Ny   )# Output partition

const xspan = Δx₁/2  : Δx₁ : L;           # Centers of cells
const xc_in = (L/Nu)*((1/2) : 1 : Nu)     # Centers of actuators
const xc_out = (L/Ny)*((1/2) : 1 : Ny)    # Centers of sensors

# Decision of scenario: actuators' spatial characterization
if scenario == 1
    in_coeff = (m = 1.0, M = 0.0, ν = 4.0 )
else
    in_coeff = (m = 1.0, M = 30.0, ν = 4.0 )
end

# Sensors' spatial characterization
out_coeff = (m = 1.0, M = 10.0, ν = 4.0 )

function characterization(x, xcenter, m::Number=1.0, M::Number=1.0, ν::Number=1.0)
    return  m*exp(-1 * norm( M * (x - xcenter) )^ν )
end

# Input characterization
for i = 1 : Nu
    β_in[round(Int, (i-1)*N₁/Nu)+1:round(Int, i*N₁/Nu), i] .= 1;
    b_in[:,i] = characterization.(xspan, xc_in[i], in_coeff.m, in_coeff.M , in_coeff.ν)
end

b_in = β_in .* b_in;

# Output characterization
for i = 1 : Ny
    γ_out[round(Int, (i-1)*N₁/Ny)+1:round(Int, i*N₁/Ny), i] .= 1;
    g_out[:,i] = characterization.(xspan, xc_in[i], out_coeff.m, out_coeff.M , out_coeff.ν)
end

g_out = γ_out .* g_out;
const g_weight = inv.(sum(transpose(g_out), dims=2))

# Controller
const Kp = ones(Nu) * 10^(4);     # Controller gain
const u_in  = zeros(Nu)           # System input signals
const y_out = zeros(Ny)           # System output signals
const y_ref = ones(Ny) * 400.0;   # Reference temperatures

# History of input signals
global u_hist
u_hist = zeros(Nu,0);

# History of measurements
global y_hist
y_hist = zeros(Ny,0)

# History of time steps
const t_hist = Real[];

# System dynamics: approximated spatial derivatives
const Mx1_lo = spdiagm(-1 => ones(N₁*N₂-1), 0 => -1*ones(N₁*N₂))
const Mx1_hi = spdiagm( 0 => -1*ones(N₁*N₂), 1 => ones(N₁*N₂-1))

const Mx2_lo = spdiagm(-1*N₁ => ones( N₁*(N₂-1) ), 0 => -1*ones(N₁*N₂))
const Mx2_hi = spdiagm( 0 => -1*ones(N₁*N₂), N₁ => ones( N₁*(N₂-1) ))

# Boundary value description
const bound_under = sparse(I,N₁*N₂, N₁)
const bound_top = spzeros(N₁*N₂, N₁)
bound_top[(N₂-1)*N₁+1 : N₁*N₂,1:N₁] = sparse(I, N₁, N₁)

const bound_left = spzeros(N₁*N₂, N₂)
const bound_right = spzeros(N₁*N₂, N₂)

for k = 1 : N₂
    bound_left[1 + (k-1)*N₁, k] = 1
    bound_right[k*N₁, k] = 1
end

# Approximated PDE as large-scale ODE
function heat_eq(dθ,θ,p,t)
    global y_hist
    global u_hist
    push!(t_hist, t);

    # Reshape temperature vector to matrix (2-dimensional temperature field)
    θ = reshape(θ, N₁, N₂)

    # Computation of system output  signals
    y_out = g_weight .* (transpose(g_out) * θ[:,end]);

    # Controller
    u = Kp .* (y_ref - y_out)
    u = u .* (0.5*(sign.(u) .+ 1.0));

    # Saving input and output signals
    y_hist = hcat(y_hist, Matrix(y_out))
    u_hist = hcat(u_hist, u)

    # Averaging temperatures at internal cell boundaries
    θ1m = (θ[1:end-1,:] + θ[2:end,:])/2.0;          # x₁-direction
    θ2m = (θ[:, 1:end-1] + θ[:, 2:end])/2.0;        # x₂-direction

    # Specific heat capacity
    c = property.(θ, 330.0, 0.4)
    Cinv = reshape(inv.(c), N₁* N₂,1) .* Id12

    # Thermal conductivity
    λ₁ = property.( θ1m , 10., 0.1); # x₁-direction
    λ₂ = property.( θ2m , 10., 0.1); # x₂-direction

    # Add addional "virtual" cell at left and right boundary
    λ₁_ext1 = cat( zeros(1,N₂), λ₁, dims=1)
    λ₁_ext2 = cat( λ₁, zeros(1,N₂), dims=1)

    # Reshape λ matrix of 2-dimensional temperature field to vector
    λ₁_ext1 = reshape(λ₁_ext1, N₁*N₂, 1 )
    λ₁_ext2 = reshape(λ₁_ext2, N₁*N₂, 1 )

    # Add addional "virtual" cell at underside and topside boundary
    λ₂_ext1 = cat( zeros(N₁,1), λ₂, dims=2 )
    λ₂_ext2 = cat( λ₂, zeros(N₁,1), dims=2 )

    # Reshape λ matrix of 2-dimensional temperature field to vector
    λ₂_ext1 = reshape(λ₂_ext1, N₁*N₂, 1 )
    λ₂_ext2 = reshape(λ₂_ext2, N₁*N₂, 1 )

    # Calculate diffusion
    diff_x1 = 1/(Δx₁^2) * (λ₁_ext1 .* Mx1_lo +  λ₁_ext2 .* Mx1_hi)
    diff_x2 = 1/(Δx₂^2) * (λ₂_ext1 .* Mx2_lo +  λ₂_ext2 .* Mx2_hi)

    # Emitted heat flux
    ϕout1 = bound_left * boundary_heatflux.(θ[1,:], θamb, h, k) + bound_right * boundary_heatflux.(θ[end,:], θamb, h, k) # x₁-direction
    ϕout2 = bound_top * boundary_heatflux.(θ[:,end], θamb, h, k) # x₂-direction

    # Induced heat flux
    ϕin = bound_under * b_in * u

    # Boundary conditions as vectors
    in_out_x1 = 1/Δx₁ * ϕout1
    in_out_x2 = 1/Δx₂ * (ϕout2 + ϕin)

    # Reshape matrix (2D temperature field) to vector
    θ = reshape(θ, N₁*N₂, 1 )

    # Calculate temporal derivative
    dθ .=  ρinv * Cinv * ( (diff_x1 + diff_x2) * θ + in_out_x1 + in_out_x2 )
end

# Sampling time and final time
const Δt = 10^(-3);
const Tfinal = 10.;
const tspan = (0, Tfinal)

# Solve the ODE
prob = ODEProblem( heat_eq, θinit,tspan )
sol = solve(prob,Euler(),dt=Δt,progress=true,save_everystep=false,save_start=true)

# Save and export results
import DelimitedFiles
basepath        = "results/";

input_filename  = basepath * "input_history_scenario_" * string(scenario) * ".txt";
output_filename = basepath *  "output_history_scenario_" * string(scenario) * ".txt";
temperature2D_filename = basepath *  "temperature2D_scenario_" * string(scenario) * ".txt";

# Save signals with time steps (t_hist)
u_io_data = [t_hist transpose(u_hist)]
y_io_data = [t_hist transpose(y_hist)]

# Save temperature as matrix (2D temperature field)
temperature2D = reshape(sol.u[2], N₁, N₂);

# Export files
DelimitedFiles.writedlm(input_filename, u_io_data, "\t");
DelimitedFiles.writedlm(output_filename, y_io_data, "\t");
DelimitedFiles.writedlm(temperature2D_filename, temperature2D, "\t");


text:
include("bondad.jl")
include("rna.jl")

using Random
using Random:seed!
using ScikitLearn;
@sk_import svm: SVC
@sk_import tree: DecisionTreeClassifier
@sk_import neighbors: KNeighborsClassifier

function crossvalidation(N::Int64, k::Int64)
    indices = repeat(1:k, Int64(ceil(N/k)));
    indices = indices[1:N];
    shuffle!(indices);
    return indices;
end;

function modelCrossValidation(modelType::Symbol, modelHyperparameters::Dict, inputs::Array{Float64,2}, targets::Array{Any,1}, numFolds::Int64)

    # Comprobamos que el numero de patrones coincide
    @assert(size(inputs,1)==length(targets));

    # Que clases de salida tenemos
    # Es importante calcular esto primero porque se va a realizar codificacion one-hot-encoding varias veces, y el orden de las clases deberia ser el mismo siempre
    classes = unique(targets);

    # Primero codificamos las salidas deseadas en caso de entrenar RR.NN.AA.
    if modelType==:ANN
        targets = oneHotEncoding(targets, classes);
    end;

    # Creamos los indices de crossvalidation
    crossValidationIndices = crossvalidation(size(inputs,1), numFolds);

    # Creamos los vectores para las metricas que se vayan a usar
    # En este caso, solo voy a usar precision y F1, en otro problema podrían ser distintas
    testAccuracies = Array{Float64,1}(undef, numFolds);
    testF1         = Array{Float64,1}(undef, numFolds);

    #max_value=0;

    # Para cada fold, entrenamos
    for numFold in 1:numFolds

        # Si vamos a usar unos de estos 3 modelos
        if (modelType==:SVM) || (modelType==:DecisionTree) || (modelType==:kNN)

            # Dividimos los datos en entrenamiento y test
            trainingInputs    = inputs[crossValidationIndices.!=numFold,:];
            testInputs        = inputs[crossValidationIndices.==numFold,:];
            trainingTargets   = targets[crossValidationIndices.!=numFold];
            testTargets       = targets[crossValidationIndices.==numFold];

            if modelType==:SVM
                model = SVC(kernel=modelHyperparameters["kernel"], degree=modelHyperparameters["kernelDegree"], gamma=modelHyperparameters["kernelGamma"], C=modelHyperparameters["C"]);
            elseif modelType==:DecisionTree
                model = DecisionTreeClassifier(max_depth=modelHyperparameters["maxDepth"], random_state=1);
            elseif modelType==:kNN
                model = KNeighborsClassifier(modelHyperparameters["numNeighbors"]);
            end;

            # Entrenamos el modelo con el conjunto de entrenamiento
            model = fit!(model, trainingInputs, trainingTargets);

            # Pasamos el conjunto de test
            testOutputs =convert(Array{Any,1},predict(model, testInputs));

            # Calculamos las metricas correspondientes con la funcion desarrollada en la practica anterior
            (acc, _, _, _, _, _, F1, _) = confusionMatrix(testOutputs, testTargets);
            #=
            if (max_value<acc)
                max_value = acc;
                printConfusionMatrix(testOutputs, testTargets);
            end;
            =#
        else

            # Vamos a usar RR.NN.AA.
            @assert(modelType==:ANN);

            # Dividimos los datos en entrenamiento y test
            trainingInputs    = inputs[crossValidationIndices.!=numFold,:];
            testInputs        = inputs[crossValidationIndices.==numFold,:];
            trainingTargets   = targets[crossValidationIndices.!=numFold,:];
            testTargets       = targets[crossValidationIndices.==numFold,:];

            # Como el entrenamiento de RR.NN.AA. es no determinístico, hay que entrenar varias veces, y
            #  se crean vectores adicionales para almacenar las metricas para cada entrenamiento
            testAccuraciesEachRepetition = Array{Float64,1}(undef, modelHyperparameters["numExecutions"]);
            testF1EachRepetition         = Array{Float64,1}(undef, modelHyperparameters["numExecutions"]);

            # Se entrena las veces que se haya indicado
            for numTraining in 1:modelHyperparameters["numExecutions"]

                if modelHyperparameters["validationRatio"]>0

                    # Para el caso de entrenar una RNA con conjunto de validacion, hacemos una división adicional:
                    #  dividimos el conjunto de entrenamiento en entrenamiento+validacion
                    #  Para ello, hacemos un hold out
                    (trainingIndices, validationIndices) = holdOut(size(trainingInputs,1), modelHyperparameters["validationRatio"]*size(trainingInputs,1)/size(inputs,1));
                    # Con estos indices, se pueden crear los vectores finales que vamos a usar para entrenar una RNA

                    # Entrenamos la RNA, teniendo cuidado de codificar las salidas deseadas correctamente
                    ann, = trainClassANN(modelHyperparameters["topology"],
                        trainingInputs[trainingIndices,:],   trainingTargets[trainingIndices,:],
                        trainingInputs[validationIndices,:], trainingTargets[validationIndices,:],
                        testInputs,                          testTargets;
                        maxEpochs=modelHyperparameters["maxEpochs"], learningRate=modelHyperparameters["learningRate"], maxEpochsVal=modelHyperparameters["maxEpochsVal"]);

                else

                    # Si no se desea usar conjunto de validacion, se entrena unicamente con conjuntos de entrenamiento y test,
                    #  teniendo cuidado de codificar las salidas deseadas correctamente
                    ann, = trainClassANN(modelHyperparameters["topology"],
                        trainingInputs, trainingTargets,
                        testInputs,     testTargets;
                        maxEpochs=modelHyperparameters["maxEpochs"], learningRate=modelHyperparameters["learningRate"]);

                end;

                # Calculamos las metricas correspondientes con la funcion desarrollada en la practica anterior
                (testAccuraciesEachRepetition[numTraining], _, _, _, _, _, testF1EachRepetition[numTraining], _) = confusionMatrix(collect(ann(testInputs')'), testTargets);
                #=
                if (max_value<testAccuraciesEachRepetition[numTraining])
                    max_value = testAccuraciesEachRepetition[numTraining];
                    printConfusionMatrix(collect(ann(testInputs')'), testTargets);
                end;
                =#
            end;

            # Calculamos el valor promedio de todos los entrenamientos de este fold
            acc = mean(testAccuraciesEachRepetition);
            F1  = mean(testF1EachRepetition);

        end;

        # Almacenamos las 2 metricas que usamos en este problema
        testAccuracies[numFold] = acc;
        testF1[numFold]         = F1;

        #println("Results in test in fold ", numFold, "/", numFolds, ": accuracy: ", 100*testAccuracies[numFold], " %, F1: ", 100*testF1[numFold], " %");

    end; # for numFold in 1:numFolds
    println()
    println(modelType, ": Average test accuracy on a ", numFolds, "-fold crossvalidation: ", 100*mean(testAccuracies), ", with a standard deviation of ", 100*std(testAccuracies));
    println(modelType, ": Average test F1 on a ", numFolds, "-fold crossvalidation: ", 100*mean(testF1), ", with a standard deviation of ", 100*std(testF1));

    return (mean(testAccuracies), std(testAccuracies), mean(testF1), std(testF1));

end;


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule cif_api_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("cif_api")
JLLWrappers.@generate_main_file("cif_api", UUID("6fcef0ae-1c05-5fc1-b206-1cf994addbad"))
end  # module cif_api_jll


text:
@testset "Learning" begin
  @testset "Tasks" begin
    t = RegressionTask(:x,:y)
    @test inputvars(t) == (:x,)
    @test outputvars(t) == (:y,)
    @test features(t) == (:x,)
    @test label(t) == :y
    @test sprint(show, t) == "Regression x → y"
    t = RegressionTask([:x,:y], :z)
    @test inputvars(t) == (:x,:y)
    @test features(t) == (:x,:y)
    @test sprint(show, t) == "Regression (x, y) → z"

    t = ClassificationTask(:x,:y)
    @test inputvars(t) == (:x,)
    @test outputvars(t) == (:y,)
    @test features(t) == (:x,)
    @test label(t) == :y
    @test sprint(show, t) == "Classification x → y"
    t = ClassificationTask([:x,:y], :z)
    @test inputvars(t) == (:x,:y)
    @test features(t) == (:x,:y)
    @test sprint(show, t) == "Classification (x, y) → z"

    t = ClusteringTask(:x,:y)
    @test inputvars(t) == (:x,)
    @test outputvars(t) == (:y,)
    @test features(t) == (:x,)
    @test sprint(show, t) == "Clustering x → y"
    t = ClusteringTask([:x,:y], :z)
    @test inputvars(t) == (:x,:y)
    @test features(t) == (:x,:y)
    @test sprint(show, t) == "Clustering (x, y) → z"
  end
end


text:
type LAMMPSDAT <: AtomFileType end

DAT = LAMMPSDAT

file_exts(::Type{LAMMPSDAT}) = (".dat",)

function Atoms(::Type{LAMMPSDAT}, num_atoms::Int)
    Atoms(num_atoms,
        :molecule => Vector{Int}(num_atoms),
        :type => Vector{Symbol}(num_atoms),
        :charge => Vector{Float64}(num_atoms),
        :coords => Matrix{Float64}(3, num_atoms)
    )
end

function read_frame_header{T}(::Type{T}, ::Type{LAMMPSDAT}, chunk::Chunk)
    num_atoms_regex = r"(\d+)\s*atoms"
    num_bonds_regex = r"(\d+)\s*bonds"
    x_dims_regex = r"(\S+)\s+(\S+)\s+xlo\s+xhi"
    y_dims_regex = r"(\S+)\s+(\S+)\s+ylo\s+yhi"
    z_dims_regex = r"(\S+)\s+(\S+)\s+zlo\s+zhi"
    masses_regex = r"Masses"
    atoms_regex = r"Atoms"
    bonds_regex = r"Bonds"

    frame_start = chunk.i

    num_atoms = parse(Int, read_through_regex(chunk, num_atoms_regex)[1])
    num_bonds = parse(Int, read_through_regex(chunk, num_bonds_regex)[1])

    x_dims = read_through_regex(chunk, x_dims_regex)
    x_lo = parse(Float64, x_dims[1])
    x_hi = parse(Float64, x_dims[2])

    y_dims = read_through_regex(chunk, y_dims_regex)
    y_lo = parse(Float64, y_dims[1])
    y_hi = parse(Float64, y_dims[2])

    z_dims = read_through_regex(chunk, z_dims_regex)
    z_lo = parse(Float64, z_dims[1])
    z_hi = parse(Float64, z_dims[2])

    read_through_regex(chunk, atoms_regex)
    read_newline(chunk)

    atom_data_start = chunk.i
    read_newlines(chunk, num_atoms)
    atom_data_size = chunk.i - atom_data_start

    read_through_regex(chunk, bonds_regex)
    read_newline(chunk)

    bond_data_start = chunk.i
    read_newlines(chunk, num_bonds)
    bond_data_size = chunk.i - bond_data_start

    frame_end = endof(chunk.data)
    frame_size = frame_end - frame_start + 1
    chunk.i = frame_end

    frame_header = FrameHeader(
        frame_start,
        frame_size,
        num_atoms,
        atom_data_start,
        atom_data_size,
        num_bonds,
        bond_data_start,
        bond_data_size
    )

    frame = Frame(
        :box_min => (x_lo, y_lo, z_lo),
        :box_max => (x_hi, y_hi, z_hi)
    )

    return frame_header, frame
end

function read_atoms{T}(::Type{T}, ::Type{LAMMPSDAT}, chunk::Chunk, atoms::Atoms, indices::UnitRange{Int})
    molecules = atoms.props[:molecule]::Vector{Int}
    types = atoms.props[:type]::Vector{Symbol}
    charges = atoms.props[:charge]::Vector{T}
    coords = atoms.props[:coords]::Matrix{T}

    for i in indices
        read_whitespace(chunk)
        read_nonwhitespace(chunk)
        read_whitespace(chunk)

        molecule = read_int(chunk)
        read_whitespace(chunk)

        atom_type = read_symbol(chunk)
        read_whitespace(chunk)

        charge = read_float(T, chunk)
        read_whitespace(chunk)

        x = read_float(T, chunk)
        read_whitespace(chunk)

        y = read_float(T, chunk)
        read_whitespace(chunk)

        z = read_float(T, chunk)
        read_newline(chunk)

        molecules[i] = molecule
        types[i] = atom_type
        charges[i] = charge
        coords[1, i] = x
        coords[2, i] = y
        coords[3, i] = z
    end
end

function read_bonds(::Type{LAMMPSDAT}, chunk::Chunk, bonds::Matrix{Int}, indices::UnitRange{Int})
    num_bonds = size(bonds, 2)

    for i in indices
        read_whitespace(chunk)
        read_nonwhitespace(chunk)
        read_whitespace(chunk)

        read_nonwhitespace(chunk)
        read_whitespace(chunk)

        atom1 = read_int(chunk)
        read_whitespace(chunk)

        atom2 = read_int(chunk)
        read_newline(chunk)

        bonds[1, i] = atom1
        bonds[2, i] = atom2
    end
end

function write_atoms(::Type{LAMMPSDAT}, io::IO, frame::Frame)
    haskey(frame, :atoms) || throw(ArgumentError("Need frame atoms"))
    haskey(frame, :box_min) || throw(ArgumentError("Need frame box min"))
    haskey(frame, :box_max) || throw(ArgumentError("Need frame box max"))
    atoms = frame[:atoms]

    haskey(atoms, :coords) || throw(ArgumentError("Need atom coords"))
    haskey(atoms, :molecule) || throw(ArgumentError("Need atom molecule IDs)"))
    haskey(atoms, :type) || throw(ArgumentError("Need atom types"))
    haskey(atoms, :charge) || throw(ArgumentError("Need atom charges)"))
    haskey(atoms, :coords) || throw(ArgumentError("Need atom coords"))

    num_atoms = length(atoms)
    unique_elements = sort!(unique(atoms[:type]))
    element_ids = Dict{Symbol, Int}(zip(unique_elements, eachindex(unique_elements)))
    element_masses = Dict(:C => 12.0107, :Cl => 35.453, :H => 1.00794, :Li => 6.941, :N => 14.0067, :O => 15.9994, :P => 30.973762, :S => 32.065)
    xlo, ylo, zlo = frame[:box_min]
    xhi, yhi, zhi = frame[:box_max]

    println(io,
        """
        LAMMPS data file
        $num_atoms atoms
        $(length(unique_elements)) atom types
        $xlo $xhi xlo xhi
        $ylo $yhi ylo yhi
        $zlo $zhi zlo zhi

        Masses
        """
    )

    for (id, element) in enumerate(unique_elements)
        println(io, id, ' ', element_masses[element], " # ", element)
    end
    println(io)

    println(io, "Atoms")
    println(io)

    for (i, atom) in enumerate(atoms)
        print(io,
            i, ' ',
            atom[:molecule], ' ',
            element_ids[atom[:type]], ' ',
            atom[:charge], ' ',
            atom[:coords][1], ' ',
            atom[:coords][2], ' ',
            atom[:coords][3], '\n'
        )
    end
end

max_stars_repo_name: NickMcNutt/Alpine

max_issues_repo_name: NickMcNutt/Alpine

max_forks_repo_name: NickMcNutt/Alpine


text:
# Physical constants
const ħ = 1.05457162825e-34; # Reduced Planck's constant (kg m^2 s^{-1})
const eV = 1.602176487e-19; # Electron Volt (kg m^2 s^{-2})
const m_e = 9.10938188e-31; # Electron Mass (kg)
const k_B = 1.3806504e-23; # Boltzmann's constant (kg m^2 K^{-1} s^{-2})
const ϵ_0 = 8.854e-12 # Dielectric constant (C^2 N^{-1} m^{-2})
const c = 2.99792458e8 # Speed of light (m s^{-1})
const amu = 1.66053906660e-27 # Atomic Mass Unit (kg)

function ϵ_ionic_mode(phonon_mode_freq, ir_activity, volume)
    ω_ν = 2π * phonon_mode_freq * 1e12
    ϵ_mode = eV^2 * ir_activity / (3 * volume * ω_ν^2 * amu)
    return ϵ_mode / ϵ_0
end

function frohlich_α_ν(ϵ_optic, ϵ_ionic, ϵ_total, phonon_mode_freq, m_eff)
    Ry = eV^4 * m_e / (2 * ħ^2)
    ω = 2π * 1e12 * phonon_mode_freq
    ϵ_static = ϵ_optic + ϵ_total
    (m_eff * Ry / (ħ * ω))^(1 / 2) * ϵ_ionic / (4π * ϵ_0 * ϵ_optic * ϵ_static)
end

MAPI= [
# 96.20813558773261 0.4996300522819191
# 93.13630357703363 1.7139631746083817
# 92.87834578121567 0.60108592692181
# 92.4847918585963 0.0058228799414729
# 92.26701437594754 0.100590086574602
# 89.43972834606603 0.006278895133832249
# 46.89209141511332 0.2460894564364346
# 46.420949316788 0.14174282581124137
# 44.0380222871706 0.1987196948553428
# 42.89702947649343 0.011159939465770681
# 42.67180170168193 0.02557751102757614
# 41.46971205834201 0.012555230726601503
# 37.08982543385215 0.00107488277468418
# 36.53555265689563 0.02126940080871224
# 30.20608114002676 0.009019481779712388
# 27.374810898415028 0.03994453721421388
# 26.363055017011728 0.05011922682554448
# 9.522966890022039 0.00075631870522737
# 4.016471586720514 0.08168931020200264
# 3.887605410774121 0.006311654262282101
# 3.5313112232401513 0.05353548710183397
2.755392921480459 0.021303020776321225
2.4380741812443247 0.23162784335484837
2.2490917637719408 0.2622203718355982
2.079632190634424 0.23382298607799906
# 2.0336707697261187 0.0623239656843172
# 1.5673011873879714 0.0367465760261409
# 1.0188379384951798 0.0126328938653956
# 1.0022960504442775 0.006817361620021601
# 0.9970130778462072 0.0103757951973341
# 0.9201781906386209 0.01095811116040592
# 0.800604081794174 0.0016830270365341532
# 0.5738689505255512 0.00646428491253749
# 0.022939578929507105 8.355742795827834e-06   # Acoustic modes!
# 0.04882611767873102 8.309858592685e-06
# 0.07575149723846182 2.778248540373041e-05
]

phonon_freqs = MAPI[:, 1]
ir_activity = MAPI[:, 2]
ir_total = sum(ir_activity)

function ϵ_total(phonon_freqs, ir_activity, volume)
    result = 0.0
    for (f, r) in zip(phonon_freqs, ir_activity)
        result += ϵ_ionic_mode(f, r, volume)
    end
    return result
end

α = 0
s = 0
ϵ_t = ϵ_total(phonon_freqs, ir_activity, (6.29E-10)^3)
for (f, r) in eachrow(MAPI)
    global α
    global s
    ϵ_ionic = ϵ_ionic_mode(f, r, (6.29E-10)^3)
    α_ν = frohlich_α_ν(4.5, ϵ_ionic, ϵ_t, f, 0.12)
    # println("f: $f, ir: $r, α_ν: $α_ν, α: $α, ϵ: $ϵ_ionic")
    α += α_ν
    s += ϵ_ionic / ϵ_t
end
println("α = $(α), $s")

using QuadGK

function κ_i(i, v, w)
    κ = v[i]^2 - w[i]^2
    if length(v) > 1
        for j in 1:length(v)
            if j != i
                κ *= (v[j]^2 - w[i]^2) / (w[j]^2 - w[i]^2)
            end
        end
    end
    return κ
end

function h_i(i, v, w)
    h = v[i]^2 - w[i]^2
    if length(v) > 1
        for j in 1:length(v)
            if j != i
                h *= (w[j]^2 - v[i]^2) / (v[j]^2 - v[i]^2)
            end
        end
    end
    return h
end

function C_ij(i, j, v, w)
    C = w[i] * κ_i(i, v, w) * h_i(j, v, w) / (4 * (v[j]^2 - w[i]^2))
    return C
end

function D_ν(τ, β, v, w)
    D = τ * (1 - τ / β)
    for i in 1:length(v)
        D += (h_i(i, v, w) / v[i]^2) * (2 * sinh(v[i] * τ / 2) * sinh(v[i] * (β - τ) / 2) / (v[i] * sinh(v[i] * β / 2)) - τ * (1 - τ / β))
    end
    return D
end

function multi_free_energy(v, w, T, ϵ_optic, ϵ_static, phonon_freqs, m_eff, ir_activity, volume)

    ϵ_tot = ϵ_total(phonon_freqs, ir_activity, volume)

    S_integrand(τ, β) = cosh(β / 2 - τ) / (sinh(β / 2) * sqrt(D_ν(τ, β, v, w)))
    S(β, α) = β * α / √π * QuadGK.quadgk(τ -> S_integrand(τ, β), 0, β / 2)[1]

    C(i, j) = C_ij(i, j, v, w)

    function S_0(β)
        s = 0.0
        for i in 1:length(v)
            for j in 1:length(v)
                s += C(i, j) / (v[j] * w[i]) * (coth(β * v[i] / 2) - 2 / (β * v[j]))
            end
        end
        -3 * s * β
    end

    function log_Z_0(β)
        s = 0.0
        for i in 1:length(v)
            s += log(v[i] * sinh(β * w[i] / 2) / (w[i] * sinh(β * v[i] / 2))) - log(2π * β) / 2
        end
        3 * s
    end

    F = 0.0
    for ν in 1:length(phonon_freqs)
        ω = 2π * 1e12 * phonon_freqs[ν]
        βred = ħ * ω / (k_B * T)
        β = ħ * 2π * 1e12 / (k_B * T)
        ϵ_ionic = ϵ_ionic_mode(phonon_freqs[ν], ir_activity[ν], volume)
        α = frohlich_α_ν(ϵ_optic, ϵ_ionic, ϵ_tot, phonon_freqs[ν], m_eff)
        println("$α, $β, $βred, $(phonon_freqs[ν]), $(-log_Z_0(β) / β), $(-S(βred, α) / β), $(S_0(β) / β), $(-(log_Z_0(β) + S(βred, α) - S_0(β)) / β)")
        F += -(log_Z_0(β) + S(βred, α) - S_0(β)) / β * 2π * 1e12
    end
    return F * 1e3 * ħ / eV
end

F = multi_free_energy([19.6], [16.3], 300, 4.5, 24.1, phonon_freqs, 0.12, ir_activity, (6.29E-10)^3)
@show(F)
# 2.2490917637719408 0.2622203718355982

using Optim
using LineSearches

function multi_variation(T, ϵ_optic, ϵ_static, phonon_freqs, m_eff, ir_activity, volume; N = 1)

    # Intial guess for v and w.
    initial = sort(rand(2 * N)) .* 10.0

    # Limits of the optimisation.
    lower = repeat([0.0], 2 * N)
    upper = repeat([Inf], 2 * N)

    # @show(initial, [initial[2 * n - 1] for n in 1:Int(N)], [initial[2 * n] for n in 1:Int(N)])

    # Osaka Free Energy function to minimise.
    f(x) = multi_free_energy([x[2 * n] for n in 1:Int(N)], [x[2 * n - 1] for n in 1:Int(N)], T, ϵ_optic, ϵ_static, phonon_freqs, m_eff, ir_activity, volume)

    err = false
    # Use Optim to optimise the free energy function w.r.t v and w.
    solution = try
        Optim.optimize(
            Optim.OnceDifferentiable(f, initial; autodiff = :forward),
            lower,
            upper,
            initial,
            Fminbox(BFGS(; linesearch = LineSearches.BackTracking())),
        )
    catch
        err = true
    end

    # # If optimisation does not converge or if v ≤ w, pick new random starting guesses for v and w between 1.0 and 11.0. Repeat until the optimisation converges with v > w.
    if err
        while err
            initial = sort(rand(2 * N)) .* 10.0
            err = false
            solution = try
                Optim.optimize(
                    Optim.OnceDifferentiable(f, initial; autodiff = :forward),
                    lower,
                    upper,
                    initial,
                    Fminbox(BFGS(; linesearch = LineSearches.BackTracking())),
                )
            catch
                err = true
            end
        end
    end

    # Return variational parameters that minimise the free energy.
    return Optim.minimizer(solution)
end

# N = 1
# rs = sort(multi_variation(300, 4.5, 24.1, phonon_freqs, 0.12, ir_activity, (6.29E-10)^3; N = N))
# F = multi_free_energy([r[2 * n] for n in 1:N], [r[2 * n - 1] for n in 1:N], 300, 4.5, 24.1, phonon_freqs, 0.12, ir_activity, (6.29E-10)^3)
# println("$rs, $(F)")