text:
# functions that do not call the PUMI shared library, but are related to those
# function

export getAdjacentFull, resetAllIts2, countDownwards, countAllNodes, printEdgeVertNumbers, printFaceVertNumbers,  getValues2, getLocalGradients2, getJacobian2, getNodeEntities, getEntityLocalNumber, printElementVertNumbers, isShared_sharing

export VERTEX, EDGE, TRIANGLE, QUAD, TET, HEX, PRISM, PYRAMIX, simplexTypes, getTypeDimension, getDimension

export ElementNodeEntities

# declare the enums
global const VERTEX=0
global const EDGE=1
global const TRIANGLE=2
global const QUAD=3
global const TET=4
global const HEX=5
global const PRISM=6
global const PYRAMID=7

global const simplexTypes = [VERTEX, EDGE, TRIANGLE, TET]

_tmp1, _tmp2, _tmp3 = getTopologyMaps()
global const tri_edge_verts = _tmp1 + 1
global const tet_edge_verts = _tmp2 + 1
global const tet_tri_verts = _tmp3 + 1

# this is not exported, use the getter function
global const typeDimension = Array{Int}(8)
typeDimension[VERTEX + 1]   = 0
typeDimension[EDGE + 1]     = 1
typeDimension[TRIANGLE + 1] = 2
typeDimension[QUAD + 1]     = 2
typeDimension[TET + 1]      = 3
typeDimension[HEX + 1]      = 3
typeDimension[PRISM + 1]    = 3
typeDimension[PYRAMID + 1]  = 3

"""
  Returns the dimension of a given type (one of the apf::Type enums)
"""
function getTypeDimension(typ::Integer)

  return typeDimension[typ + 1]
end

function getDimension(m_ptr::Ptr{Void}, entity::Ptr{Void})

  typ = getType(m_ptr, entity)
  return getTypeDimension(typ)
end

function getAdjacentFull(m_ptr, entity, dimension::Integer)
# returns an array with the adjacencies of meshentity entity of specified dimension, and the number of entries in the array
# this encasuplates the calls to countAdjacent and getAdjacent

  n = countAdjacent(m_ptr, entity, dimension)
  adj = getAdjacent(n)

  return adj, n

end

function countDownwards(m_ptr, entity)
# gets an array containing the number of downward adjacencies of each type for
# the given entity, in ascending order
# ie. array[1] = # of vertices, array[3] = # of faces

  entity_type = getType(m_ptr, entity)

  if entity_type == 0  # vertex
    return [0]  # return a non empty array

  elseif entity_type == 1  # edge
    return [2]

  elseif entity_type == 2 # triangle
    return [3 3]
  elseif entity_type == 3 # quadralaterial
    return [4, 4]
  elseif entity_type == 4 # tet
    return [4, 6, 4]
  elseif entity_type == 5 # hex (wedge?)
    return [6, 9, 5]
  elseif entity_type == 6 # prism (brick?)
    return [8, 12, 6]
  elseif entity_type == 7 #pyramid
    return [5, 8, 5]
  else
    println("unrecognized entity type")
    return [0]
  end

end

function countAllNodes(mshape_ptr, entity_type::Integer)
# gets the total number of nodes on an entity type, including downward dependencies

  eshape_ptr = getEntityShape(mshape_ptr, entity_type)
  i = countNodes(eshape_ptr)
  return i
end

function printEdgeVertNumbers(edgeN_ptr, vertN_ptr; fstream=STDOUT)
# prints the number of the verteces on the edge
# edgeN_ptr is a pointer to the numbering of edges
# vertN_ptr is a pointer to the number of vertices
# numEdges is the number of edges in the mesh
# fstream is the stream the data is written to, default STDOUT
# this function uses iterators

  m_ptr = getNumberingMesh(edgeN_ptr)  # get pointer to mesh
  it = MeshIterator(m_ptr, 1)
  m = countJ(m_ptr, 1)  # count number of edges on the mesh
  println("m = ", m)

  for i=1:m  # loop over edges
    edge_i = iterate(it)
    edge_num = getNumberJ(edgeN_ptr, edge_i, 0, 0)

    (verts, nverts) = getDownward(m_ptr, edge_i, 0)
    n1 = getNumberJ(vertN_ptr, verts[1], 0, 0)
    n2 = getNumberJ(vertN_ptr, verts[2], 0, 0)
    println(fstream, n1, " ", n2)
  end

  free(it)

end

function printFaceVertNumbers(edges::AbstractArray{Ptr{Void}}, edgeN_ptr, vertN_ptr; fstream=STDOUT)

  m_ptr = getNumberingMesh(edgeN_ptr)
  m = length(edges)

  for i=1:m
    edge_i = edges[i]
    edge_num = getNumberJ(edgeN_ptr, edge_i, 0, 0)
#    print(fstream, "face ", edge_num, " has vertices ")
    (verts, nverts) = getDownward(m_ptr, edge_i, 0)
    for j=1:nverts
      n = getNumberJ(vert_Nptr, verts[j], 0, 0)
      print(fstream, n1, " ")
    end
    print(fstream, "\n")
  end

end

function printFaceVertNumbers(faceN_ptr, vertN_ptr; fstream=STDOUT)
#print the numbers of the vertices that compose a face
  m_ptr = getNumberingMesh(faceN_ptr)
  it = MeshIterator(m_ptr, 2)
  m = countJ(m_ptr, 2)  # count number of faces on the mesh

 for i=1:m  # loop over faces
   face_i = iterate(it)
   face_num = getNumberJ(faceN_ptr, face_i, 0, 0)
   (verts, nverts) = getDownward(m_ptr, face_i, 0)
   vertnums = zeros(Int, nverts)
   for j=1:nverts
     vertnums[j] = getNumberJ(vertN_ptr, verts[j], 0, 0)
     print(fstream, vertnums[j], " ")
   end
   print(fstream, "\n")
 end

 free(it)

 return nothing
end

function printElementVertNumbers(el_Nptr, vert_Nptr; fstream=STDOUT)
  m_ptr = getNumberingMesh(el_Nptr)
  it = MeshIterator(m_ptr, 3)
  m = countJ(m_ptr, 3)  # count number of element

  for i=1:m
    el_i = iterate(it)
    elnum = getNumberJ(el_NPtr, el_i, 0, 0)
    (verts, nverts) = getDownward(m_ptr, el_i, 0)
    vertnums = zeros(Int, nverts)
    for j=1:nverts
      vertnums[j] = getNumberJ(vert_Nptr, verts[j], 0, 0)
      print(fstream, vertnums[j], " ")
    end
    print(fstream, "\n")
  end

  free(it)

  return nothing
end

function getValues2(eshape_ptr, coords::Array{Float64,1})
# gets the values of the shape functions at coords
# output is a vector whose length is the number of nodes affecting the entity used to get eshape_ptr

  numN = countNodes(eshape_ptr)
  vals = getValues(eshape_ptr, coords, numN)
  return vals
end

function getLocalGradients2(eshape_ptr, coords)
# get the shape function derivative values at coords
# the output is a 3 x numN matrix, where numN is the number of nodes affecting the entity usued to get eshape_ptr

   numN = countNodes(eshape_ptr)
   vals = getLocalGradients(eshape_ptr, coords, numN)
   return vals
end

function getJacobian2(m_ptr, entity, coords)
# creates a mesh element, then gets its jacobian at the coordinates

   mel_ptr = createMeshElement(m_ptr, entity)
   jac = getJacobian(mel_ptr, coords)
   return jac
end

"""
  Struct that contains all the MeshEntities that have nodes on them, as
  well as the number of nodes on each one

  It might be helpful to generalize this so it works on any MeshEntity,
  not just an element, but that would require more topology information

  **Public Fields**

   * entities: array of MeshEntities that have nodes on them, length `nentities`
   * nodecounts: array containing the number of nodes on each entity,
                 length `nentities`
   * nentities: number of entities
   * nnodes: total number of nodes on all the entities
"""
struct ElementNodeEntities
  # public fields
  entities::Vector{Ptr{Void}}
  nodecounts::Vector{Int}  # number of nodes on each entity
  dimensions::Vector{Int}  # dimension of each entity
  nentities::Int  # = length of above arrays
  nnodes::Int # total number of nodes

  # private fields
  m_ptr::Ptr{Void}
  has_nodes_in_dim::Vector{Bool}
  dim::Int

  # temporary arrays
  down_entities::Array{Ptr{Void}, 1}
end

function ElementNodeEntities(m_ptr::Ptr{Void}, mshape_ptr::Ptr{Void},
                             dim::Integer)

  if dim == 2  # triangle
    num_type_per_entity = [3, 3, 1]
    dim = 2
  elseif dim == 3  # tetrahedron
    num_type_per_entity = [4, 6, 4, 1]
    numRegions = num_type_per_entity[4]
    dim = 3
  else
    throw(ErrorException("unsupported dimension in getNodeEntites, dim = $dim"))
  end

  num_nodes_per_entity = zeros(Int, 4)
  for i=1:3
    num_nodes_per_entity[i] = countNodesOn(mshape_ptr, i-1)
  end
  num_nodes_per_entity[4] = countNodesOn(mshape_ptr, 4)  # tetrahedron

  has_nodes_in_dim = Array{Bool}(4)
  for i=1:4
    has_nodes_in_dim[i] = num_nodes_per_entity[i] != 0
  end

  # compute number of entities
  nentities = 0
  nnodes = 0
  for i=0:dim
    if num_nodes_per_entity[i+1] > 0
      nentities += num_type_per_entity[i+1]
      nnodes += num_type_per_entity[i+1]*num_nodes_per_entity[i+1]
    end
  end

  entities = Array{Ptr{Void}}(nentities)
  nodecounts = Array{Int}(nentities)
  dimensions = Array{Int}(nentities)
  down_entities = Array{Ptr{Void}}(12)  # apf::Downward

  # get the nodecounts
  idx = 1
  for i=0:dim
    if num_nodes_per_entity[i+1] > 0
      for j=1:num_type_per_entity[i+1]
        nodecounts[idx] = num_nodes_per_entity[i+1]
        dimensions[idx] = i
        idx += 1
      end
    end
  end

  @assert idx == nentities + 1

  return ElementNodeEntities(entities, nodecounts, dimensions,  nentities,
                             nnodes, m_ptr,
                             has_nodes_in_dim, dim, down_entities)
end

"""
  Gets the apf::MeshEntities of the entities that have nodes on them.
  The entities are retrieve in order: vertices, edges, faces, and
  regions.  Entities of each dimension are ordered consistently with
  Pumi's `getDownward` function.

  **Inputs**

   * obj: an [`ElementNodeEntities`](@ref) object.  The `entities` field
          will be overwritten.
   * entity: an apf::MeshEntity* for the element.  The MeshEntities of
             the downward adjacent entities that have nodes on them will
             be retrieve
"""
function getNodeEntities(obj::ElementNodeEntities, entity::Ptr{Void})

  # entity must be of the same type passed to the constructor of obj

  # get entities in order of increasing dimension
  idx = 1
  for i=1:(obj.dim+1)
    if obj.has_nodes_in_dim[i]
      n = getDownward(obj.m_ptr, entity, i-1, obj.down_entities)
      for j=1:n
        obj.entities[idx] = obj.down_entities[j]
        idx += 1
      end
    end
  end

  return nothing
end

function getNodeEntities(m_ptr, mshape_ptr, entity)

  entity_type = getType(m_ptr, entity)
  eshape_ptr = getEntityShape(mshape_ptr, entity_type)
  nnodes = countNodes(eshape_ptr)
  downward_entities = Array{Ptr{Void}}(nnodes)  # holds mesh entities
  getNodeEntities(m_ptr, mshape_ptr, entity, downward_entities)

  return downward_entities
end

global const _getNodeEntities_retrieved_entities = Array{Ptr{Void}}(12)  # reusable storage
function getNodeEntities(m_ptr, mshape_ptr, entity,
                         downward_entities::AbstractArray{Ptr{Void}})
# get the meshentities that have nodes on them
# duplicates entities that have multiple nodes
# this only works for simplex elements

  entity_type = getType(m_ptr, entity)
  eshape_ptr = getEntityShape(mshape_ptr, entity_type)
  nnodes = countNodes(eshape_ptr)
#  downward_entities = Array{Ptr{Void}}(nnodes)  # holds mesh entities
  num_vert_nodes = countNodesOn(mshape_ptr, 0)
  num_edge_nodes = countNodesOn(mshape_ptr, 1)
  num_face_nodes = countNodesOn(mshape_ptr, 2)
  num_region_nodes = countNodesOn(mshape_ptr, 4) # tetrahedron

  has_vert_nodes = (num_vert_nodes != 0)
  has_edge_nodes = (num_edge_nodes != 0)
  has_face_nodes = (num_face_nodes != 0)
  has_region_nodes = (num_region_nodes != 0)

  if entity_type == 2  # triangle
    num_type_per_entity = [3, 3, 1]
  elseif entity_type == 4  # tetrahedron
    num_type_per_entity = [4, 6, 4, 1]
    numRegions = num_type_per_entity[4]
  else
    throw(ErrorException("unsupported entity type in getNodeEntites, entity type = $entity_type"))
  end

  numVerts = num_type_per_entity[1]
  numEdges = num_type_per_entity[2]
  numFaces = num_type_per_entity[3]
#=

  if entity_type == 2  # triangle
    num_entities = 3*has_vert_nodes + 3*has_edge_nodes + has_face_nodes
  elseif entity_type == 4 # tetrahedron
    num_entities = 4*has_vert_nodes + 6*has_edge_nodes + 4*has_face_nodes + has_region_nodes
  else
    println("element type not supported by getNodeEntities")
  end

=#
  retrieved_entities = _getNodeEntities_retrieved_entities
  vert_offset = 0
  if has_vert_nodes  # vertices
    vert_offset = getDownward(m_ptr, entity, 0, downward_entities)
  end

  if has_edge_nodes # edges
    retrieved_entities[1] = Ptr{Void}(42)
    numEdges = getDownward(m_ptr, entity, 1, retrieved_entities)
        for i=1:numEdges
          insertN(downward_entities, retrieved_entities[i], vert_offset+num_edge_nodes*(i-1) + 1, num_edge_nodes)
        end
  end

  if has_face_nodes
     numFaces = getDownward(m_ptr, entity, 2, retrieved_entities)
    for i=1:numFaces
      insertN(downward_entities, retrieved_entities[i], vert_offset + num_edge_nodes*numEdges +num_edge_nodes*(i-1) + 1, num_face_nodes)
    end
  end

  if has_region_nodes
    for i=1:numRegions
      insertN(downward_entities, entity, vert_offset + num_edge_nodes*numEdges + num_face_nodes*numFaces + num_region_nodes*(i-1) + 1, num_region_nodes)
    end
  end

  return nothing
#  return downward_entities
end

function insertN(vec::AbstractArray{T}, element::T,  index::Integer, n::Integer) where T
# insert element into the vector n times, starting at index

  for i=index:(index+n-1)
    vec[i] = element
  end

  return nothing
end

function getEntityLocalNumber(mesh::Ptr{Void}, entity::Ptr{Void}, parent::Ptr{Void},  entity_dim::Integer, parent_dim::Integer)
# get the local index of the entity among all the entities of the same dimension that
# belong to the same parent
# entity_dim is the dimension of the entity
# parent_dim is the dimension of the parent
# zero based indexing

  if parent_dim <= entity_dim
    println("parent dimension must be > entity dimension")
    return -1
  end

  (down_entities, numdown) = getDownward(mesh, parent, entity_dim)

  local_num = -1
  for i=1:numdown
    if (down_entities[i] == entity)
      local_num = i-1
    end
  end

  return local_num
end

# check if an entity is shared according to the supplied sharing
# this is a bit expensive because it has to actually get the
# shared entities to find out how many there are
function isShared_sharing(shr::Ptr{Void},  entity::Ptr{Void})
  return  countCopies(shr, entity) != 0
end

#------------------------------------------------------------------------------
# Iteration over entities with nodes on them (according to a FieldShape)

import Base: start, next, done, iteratorsize, iteratoreltype, eltype, length, size
using Base: HasLength, HasEltype

"""
  This iterator iterates over MeshEntities that have nodes on them in a given
  FieldShape.  The iterator returns the MeshEntity* and its dimension.  The
  iterator does not need to be freed unless the loop is exited early using a
  `break` statement.  In this case, call `apf.free` on the iterator.

  **Constructor**

   * m_ptr: apf::Mesh*
   * fshape_ptr: apf::FieldShape*

  **Usage**

   ```
     for (e, e_dim) in apf.FieldEntityIt(m_ptr, fshape_ptr)
       # do stuff
     end
   ```
"""
mutable struct FieldEntityIt
  m_ptr::Ptr{Void}
  fshape_ptr::Ptr{Void}
  dim::Int  # current dimension
  lastdim::Int  # greatest dimension entity that has nodes on it
  mesh_dim::Int  # mesh dimension
  its::Vector{MeshIterator}
  has_nodes_in::Vector{Bool}
  len::Int  # number of entities
end

mutable struct FieldEntityItState
  e::Ptr{Void}
  dim::Int
end

using ODLCommonTools

function FieldEntityIt(m_ptr::Ptr{Void}, fshape_ptr::Ptr{Void})

  dim = 0
  lastdim = 0
  mesh_dim = getDimension(m_ptr)
  its = Vector{MeshIterator}(mesh_dim + 1)
  has_nodes_in = Vector{Bool}(mesh_dim + 1)
  for i=0:mesh_dim
    its[i+1] = MeshIterator(m_ptr, i)
    has_nodes_in[i+1] = hasNodesIn(fshape_ptr, i)
  end

  # find first dimension with nodes
  for i=1:length(has_nodes_in)
    if has_nodes_in[i]
      dim = i-1
      break
    end
  end

  for i=length(has_nodes_in):-1:1
    if has_nodes_in[i]
      lastdim = i - 1
      break
    end
  end

  len = 0
  for i=1:length(has_nodes_in)
    if has_nodes_in[i]
      len += apf.countJ(m_ptr, i-1)
    end
  end

  return FieldEntityIt(m_ptr, fshape_ptr, dim, lastdim, mesh_dim, its,
                       has_nodes_in, len)
end

function free(iter::FieldEntityIt)

  for it in iter.its
    free(it)
  end

  return nothing
end

function start(iter::FieldEntityIt)

  dim = iter.dim
  e = iterate(iter.its[dim+1])

  return FieldEntityItState(e, dim)
end

function next(iter::FieldEntityIt, state::FieldEntityItState)

  # the state always leads the current element by 1, so we can do the done
  # check
  new_val = (state.e, state.dim)
  state.e = iterate(iter.its[iter.dim+1])
  state.dim = iter.dim
  return new_val, state
end

function done(iter::FieldEntityIt, state)

  e = state.e
  isdone = e == C_NULL && iter.dim == (iter.lastdim)
  if e == C_NULL && !isdone
    # start iterating over next dimension entities
    iter.dim += 1
    state.e = iterate(iter.its[iter.dim+1])
    state.dim = iter.dim
  end

  if isdone
    free(iter)
  end

  return isdone
end

function iteratorsize(::Type{FieldEntityIt})
  return HasLength()
end

function length(iter::FieldEntityIt)
  return iter.len
end

function iteratoreltype(::Type{FieldEntityIt})
  return HasEltype()
end

function eltype(::Type{FieldEntityIt})
  return Tuple{Ptr{Void}, Int}
end


text:
module cont_opt

using Random
using Plots
using Statistics
using BlackBoxOptimizationBenchmarking
using LinearAlgebra
using ForwardDiff
include("plot_utils.jl")

const MAX_GEN = 2000
const DIM = 10
const POP_SIZE = 100
const NELITES = 20
const GENS_PER_FRAME = 100
const X_LIM = 5.0

const MUT_PROB = 1.0
const CROSS_PROB = 1.0
const ROBOT_SWAP = 0.1

#mutation variance
const MUT_VAR = 0.5

const NEXPERIMENTS = 10
const EXPERIMENT_NAME = "adapt"
const COMPARE_TO = "default"

mutable struct Individial
    fitness::Float64
    obj_val::Float64
    mut_var::Float64
    x::Vector{Float64}
    robot::Bool
end

const Population = Vector{Individial}

function randu(lo::Float64, hi::Float64)
    return (hi - lo)*rand() + lo
end

function randu(lo::Float64, hi::Float64, dim::Integer)
    return [randu(lo, hi) for _ in 1:dim]
end

function spawn()::Individial
    return Individial(-Inf, Inf, MUT_VAR, randu(-X_LIM, X_LIM, DIM), true)
end

function eval!(ind::Individial, fun::BBOBFunction)
    val = fun(ind.x)
    ind.fitness = -val
    ind.obj_val = val - fun.f_opt + 1e-15
end

function mutate!(ind::Individial)
    ind.robot = (rand() < ROBOT_SWAP ? !ind.robot : ind.robot)
    ind.mut_var = clamp(ind.mut_var*exp(MUT_VAR*randn()), 0.0, MUT_VAR)
    if rand() < MUT_PROB && !ind.robot
        for i in 1:DIM
            ind.x[i] = clamp(ind.x[i] + ind.mut_var*randn(), -X_LIM, X_LIM)
        end
    end
end

function climb!(ind::Individial, Dfun::Any)
    if ind.robot
        way = -Dfun(ind.x)
        way *= ind.mut_var/norm(way)
        ind.x += way
    end
end

function breed!(parents::Population, childs::Population)
    npairs = div(POP_SIZE, 2)
    Threads.@threads for n in 1:npairs
        parent1 = tournament_selection(parents, 2)
        parent2 = tournament_selection(parents, 2)
        child1 = childs[n]
        child2 = childs[n + npairs]
        copy!(child1, parent1)
        copy!(child2, parent2)
        cross!(parent1, parent2, child1, child2)
    end
end

function copy!(dest::Individial, src::Individial)
    dest.x .= src.x
    dest.robot = src.robot
    #dest.mut_var = src.mut_var
end

function cross!(
        parent1::Individial, parent2::Individial,
        child1::Individial, child2::Individial
        )
    if rand() < CROSS_PROB && !parent1.robot && !parent2.robot
        cross_point = rand(1:(DIM-1))
        for i in 1:cross_point
            child1.x[i] = parent2.x[i]
            child2.x[i] = parent1.x[i]
        end
    end
end

function tournament_selection(pop::Population, n = 2)::Individial
    winner = rand(pop)
    for _ in 2:n
        contestant = rand(pop)
        if (contestant.fitness > winner.fitness)
            winner = contestant
        end
    end
    return winner
end

function saveind(ind::Individial, path::String)
    open(path, "w") do file
        for i in 1:DIM
            write(file, string(ind.x[i])*"\n")
        end
    end
end

function eva(fun::BBOBFunction)
    Dfun = x -> ForwardDiff.gradient(fun, x)
    log_best = Vector{Float64}()
    #random initial population
    pop = [spawn() for _ in 1:POP_SIZE]
    #buffer for next population
    childs = [spawn() for _ in 1:POP_SIZE]
    #eval them
    for ind in pop
        eval!(ind, fun)
    end
    #then iterate:
    for gen in 1:MAX_GEN
        #generate offspring
        breed!(pop, childs)
        #mutate and calculate fitness
        Threads.@threads for ind in childs
            mutate!(ind)
            climb!(ind, Dfun)
            eval!(ind, fun)
        end

        #differential_mut!(childs)
        #share!(childs)

        #kill the parents, except elites
        Threads.@threads for n in NELITES+1:POP_SIZE
            temp = pop[n]
            pop[n] = childs[n]
            childs[n] = temp
        end
        #order by fitness
        sort!(pop, by = (ind -> ind.obj_val))
        #save results
        if (gen % GENS_PER_FRAME == 0)
            println(pop[1].obj_val, '\t', pop[1].fitness)
            println("% of robots = ", count(ind -> ind.robot, pop)/length(pop))
        end
        push!(log_best, pop[1].obj_val)
    end
    return (pop[1], log_best)
end

function main()
    bbob_funs = BlackBoxOptimizationBenchmarking.list_functions()
    for n in [2, 8]
        println()
        println("BBOB function #", n)
        println("============================================")
        log = zeros(NEXPERIMENTS, MAX_GEN)
        for experiment in 1:NEXPERIMENTS
            println("experiment #", experiment)
            @time (result, log_best) = eva(bbob_funs[n])
            for gen in 1:length(log_best)
                log[experiment, gen] = log_best[gen]
            end
            #uncomment this if you want to save the best solution
            #saveind(result, "julia-partition-results/"*EXPERIMENT_NAME*"-"*string(k)*".best")
        end
        new_path = "julia-cont_opt-results/"*EXPERIMENT_NAME*string(n)*".log"
        old_path = "julia-cont_opt-results/"*COMPARE_TO*string(n)*".log"
        fig_path = "julia-cont_opt-results/"*EXPERIMENT_NAME*string(n)*".pdf"
        plot_utils.save_log(POP_SIZE*(1:MAX_GEN), log, new_path)
        if ispath(old_path)
            p = plot_utils.plot_log(old_path, new_path)
            savefig(p, fig_path)
        else
            p = plot_utils.plot_log(new_path)
            savefig(p, fig_path)
        end
    end
end

end

ext: jl

max_stars_repo_name: OndrejKincl/evaTeaching-python

max_issues_repo_name: OndrejKincl/evaTeaching-python

max_forks_repo_name: OndrejKincl/evaTeaching-python


text:
using Adapt

"""
    Field{X, Y, Z, A, G, B} <: AbstractField{X, Y, Z, A, G}

A field defined at the location (`X`, `Y`, `Z`), each of which can be either `Cell`
or `Face`, and with data stored in a container of type `A` (typically an array).
The field is defined on a grid `G` and has field boundary conditions `B`.
"""
struct Field{X, Y, Z, A, G, B} <: AbstractField{X, Y, Z, A, G}
                   data :: A
                   grid :: G
    boundary_conditions :: B

    function Field{X, Y, Z}(data, grid, bcs) where {X, Y, Z}
        validate_field_data(X, Y, Z, data, grid)
        return new{X, Y, Z, typeof(data), typeof(grid), typeof(bcs)}(data, grid, bcs)
    end
end

"""
    Field(X, Y, Z, arch, grid, [  bcs = FieldBoundaryConditions(grid, (X, Y, Z)),
                                 data = new_data(arch, grid, (X, Y, Z)) ] )

Construct a `Field` on `grid` with `data` on architecture `arch` with
boundary conditions `bcs`. Each of `(X, Y, Z)` is either `Cell` or `Face` and determines
the field's location in `(x, y, z)`.

Example
=======

julia> ω = Field(Face, Face, Cell, CPU(), RegularCartesianmodel.grid)

"""
function Field(X, Y, Z, arch, grid,
                bcs = FieldBoundaryConditions(grid, (X, Y, Z)),
               data = new_data(eltype(grid), arch, grid, (X, Y, Z)))

    return Field{X, Y, Z}(data, grid, bcs)
end

#####
##### Convenience constructor for Field that uses a 3-tuple of locations rather than a list of locations:
#####

# Type "destantiation": convert Face() to Face and Cell() to Cell if needed.
destantiate(X) = typeof(X)
destantiate(X::DataType) = X

"""
    Field(L::Tuple, arch, grid, data, bcs)

Construct a `Field` at the location defined by the 3-tuple `L`,
whose elements are `Cell` or `Face`.
"""
Field(L::Tuple, args...) = Field(destantiate.(L)..., args...)

#####
##### Special constructors for tracers and velocity fields
#####

"""
    CellField([ FT=eltype(grid) ], arch::AbstractArchitecture, grid,
              [  bcs = TracerBoundaryConditions(grid),
                data = new_data(FT, arch, grid, (Cell, Cell, Cell) ] )

Return a `Field{Cell, Cell, Cell}` on architecture `arch` and `grid` containing `data`
with field boundary conditions `bcs`.
"""
function CellField(FT::DataType, arch, grid,
                    bcs = TracerBoundaryConditions(grid),
                   data = new_data(FT, arch, grid, (Cell, Cell, Cell)))

    return Field{Cell, Cell, Cell}(data, grid, bcs)
end

"""
    XFaceField([ FT=eltype(grid) ], arch::AbstractArchitecture, grid,
               [  bcs = UVelocityBoundaryConditions(grid),
                 data = new_data(FT, arch, grid, (Face, Cell, Cell) ] )

Return a `Field{Face, Cell, Cell}` on architecture `arch` and `grid` containing `data`
with field boundary conditions `bcs`.
"""
function XFaceField(FT::DataType, arch, grid,
                     bcs = UVelocityBoundaryConditions(grid),
                    data = new_data(FT, arch, grid, (Face, Cell, Cell)))

    return Field{Face, Cell, Cell}(data, grid, bcs)
end

"""
    YFaceField([ FT=eltype(grid) ], arch::AbstractArchitecture, grid,
               [  bcs = VVelocityBoundaryConditions(grid),
                 data = new_data(FT, arch, grid, (Cell, Face, Cell)) ] )

Return a `Field{Cell, Face, Cell}` on architecture `arch` and `grid` containing `data`
with field boundary conditions `bcs`.
"""
function YFaceField(FT::DataType, arch, grid,
                     bcs = VVelocityBoundaryConditions(grid),
                    data = new_data(FT, arch, grid, (Cell, Face, Cell)))

    return Field{Cell, Face, Cell}(data, grid, bcs)
end

"""
    ZFaceField([ FT=eltype(grid) ], arch::AbstractArchitecture, grid,
               [  bcs = WVelocityBoundaryConditions(grid),
                 data = new_data(FT, arch, grid, (Cell, Cell, Face)) ] )

Return a `Field{Cell, Cell, Face}` on architecture `arch` and `grid` containing `data`
with field boundary conditions `bcs`.
"""
function ZFaceField(FT::DataType, arch, grid,
                     bcs = WVelocityBoundaryConditions(grid),
                    data = new_data(FT, arch, grid, (Cell, Cell, Face)))

    return Field{Cell, Cell, Face}(data, grid, bcs)
end

 CellField(arch::AbstractArchitecture, grid, args...) =  CellField(eltype(grid), arch, grid, args...)
XFaceField(arch::AbstractArchitecture, grid, args...) = XFaceField(eltype(grid), arch, grid, args...)
YFaceField(arch::AbstractArchitecture, grid, args...) = YFaceField(eltype(grid), arch, grid, args...)
ZFaceField(arch::AbstractArchitecture, grid, args...) = ZFaceField(eltype(grid), arch, grid, args...)

@propagate_inbounds Base.setindex!(f::Field, v, inds...) = @inbounds setindex!(f.data, v, inds...)

Adapt.adapt_structure(to, field::Field) = Adapt.adapt(to, field.data)


text:
# ---
# title: 467. Unique Substrings in Wraparound String
# id: problem467
# author: Indigo
# date: 2021-06-22
# difficulty: Medium
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/unique-substrings-in-wraparound-string/description/>
# hidden: true
# ---
#
# Consider the string `s` to be the infinite wraparound string of
# "abcdefghijklmnopqrstuvwxyz", so `s` will look like this:
# "...zabcdefghijklmnopqrstuvwxyzabcdefghijklmnopqrstuvwxyzabcd....".
#
# Now we have another string `p`. Your job is to find out how many unique non-
# empty substrings of `p` are present in `s`. In particular, your input is the
# string `p` and you need to output the number of different non-empty substrings
# of `p` in the string `s`.
#
# **Note:** `p` consists of only lowercase English letters and the size of p
# might be over 10000.
#
# **Example 1:**
#
#
#
#     Input: "a"
#     Output: 1
#
#     Explanation: Only the substring "a" of string "a" is in the string s.
#
#
# **Example 2:**
#
#
#
#     Input: "cac"
#     Output: 2
#     Explanation: There are two substrings "a", "c" of string "cac" in the string s.
#
#
# **Example 3:**
#
#
#
#     Input: "zab"
#     Output: 6
#     Explanation: There are six substrings "z", "a", "b", "za", "ab", "zab" of string "zab" in the string s.
#
#
#
## @lc code=start
using LeetCode

function find_substring_in_wrapround_string(p::String)
    dp = fill(0, 26)
    k = 0
    iscont(prev, cur) = cur == (prev == 'z' ? 'a' : prev + 1)
    for i in 1:length(p)
        k = (i > 1 && iscont(p[i - 1], p[i])) ? k + 1 : 1
        dp[p[i] - 'a' + 1] = dp[p[i] - 'a' + 1] > k ? dp[p[i] - 'a' + 1] : k
    end
    return sum(dp)
end
## @lc code=end


text:
for N in [Float64, Float32, Rational{Int}]

    # using IA types
    X = IA.interval(N(0), N(1)) # IA
    Y = Interval(N(-1), N(2))
    @test !isdisjoint(X, Y)
    @test !isdisjoint(Y, X)

    X = IntervalBox(IA.interval(N(0), N(1)), IA.interval(N(0), N(1)))
    Y = Hyperrectangle(low=[N(-1), N(-1)], high=[N(2), N(2)])
    @test !isdisjoint(X, Y)
    @test !isdisjoint(Y, X)
end


text:
using XUnit
using Distributed

include(joinpath(@__DIR__, "subdir2/subdir2-macrodef.jl"))
import .SubModule

function fn_throws()
    throw(KeyError("Test error!"))
end

function fn_nothrows()
    return nothing
end

global_param = "global_p"
extra_param_value = "extra_p"

function before_each_fn_gen(fn_name, extra_param=extra_param_value)
    () -> println("         >>>> before_each -> $fn_name, $extra_param, $global_param")
end
function after_each_fn_gen(fn_name, extra_param=extra_param_value)
    () -> println("         >>>> after_each -> $fn_name, $extra_param, $global_param")
end

println("Running tests with $(nworkers()) processes and $(Threads.nthreads()) threads")

topname = "ABC"
bottomname = "XYZ"

function set_results_to_main(testset)
    # Set Main.XUNIT_UNITTEST_RESULTS[] to communicate the results back to the test-runner
    @info "Setting Main.XUNIT_UNITTEST_RESULTS[] to '$(XUnit.get_description(testset))'"
    Main.XUNIT_UNITTEST_RESULTS[] = testset
    # html_report(testset)
    # run(`open ./$(html_output(testset))`)
end

@testset "Top XParent $topname" metrics=SubModule.MyTestMetrics runner=DistributedTestRunner() xml_report=true success_handler=set_results_to_main failure_handler=set_results_to_main begin
    @testset "XTest 1" before_each=before_each_fn_gen("Child XTest 1") begin
        @testcase "Child XTest 1 $bottomname" begin
            @test 1 == 1
            @test 2 == 2
            @test 3 == 3
        end
        @testset "Child XTest 2" begin
            @test 22 == 22
            @test 33 == 32
            @test 44 == 44
            @testset "Grand Child XTest Suite 2-10" begin
                @testcase "Grand Grand Child XTest Case 2-10-1" begin
                    @test 221 == 221
                    @test 331 == 321
                    @test 441 == 441

                    @testset "Grand^3 Child XTest Case 2-10-1-1" begin
                        @test 2211 == 2211
                        @test 3311 == 3211
                        @test 4411 == 4411
                        @testset "Grand^4 Child XTest Case 2-10-1-1-4" begin
                            @test 22114 == 22114
                            @test 33114 == 32114
                            @test 44114 == 44114
                            @testcase "Grand^5 Child XTest Case 2-10-1-1-4-6" begin
                                @test 221146 == 221146
                                @test 331146 == 321146
                                @test 441146 == 441146
                                @testset "Grand^6 Child XTest Case 2-10-1-1-4-6-9" begin
                                    @test 2211469 == 2211469
                                    @test 3311469 == 3211469
                                    @test 4411469 == 4411469

                                    @testset "Grand^7 Child XTest Case 2-10-1-1-4-6-9-3" begin
                                        @test 22114693 == 22114693
                                        @test 33114693 == 32114693
                                        @test 44114693 == 44114693
                                    end
                                end
                            end
                        end
                    end
                end
            end
            @testcase "Grand Child XTest Case 2-11 $bottomname" begin
                @test 222 == 222
                @test 332 == 322
                @test 442 == 442
                include("sub-unittests.jl")
            end
        end
        @testcase "Child XTest 3" begin
            test_println("You should only see me as part of Child XTest 3 output!")
            @test 4 == 4
            @test 5 == 5
            @test 2 * 2 == 5
            @test fn_throws()
            @test 9 == 9
            @test 10 == 10
        end
    end
    @testset "XTest 2" begin
        @testcase "Child XTest 4" begin
            @test 6 == 6
            @test 7 == 7
            @test_throws KeyError fn_throws()
            @test_throws KeyError fn_nothrows()
            @test 8 == 8
            @test 11 == 11
        end
        @testcase "Child XTest 5" begin
            @test 11 == 11
            @test 12 == 12
            @test 13 == 13
        end
        @testcase "Child XTest 12" begin
            @test 1211 == 1211
            @test 1212 == 1212 * 2
            @test 1213 == 1213
        end
    end
    @testset "XTest 3" begin
        @testset "Child XTest 3" xxx=2 before_each=() -> begin
            before_each_fn_gen("Child XTest 3-2")()
        end begin
            @testcase "Child XTest 3-1" begin
                x = 31
                y = 52
                @testset "Grand Child XTest 3-1-1" begin
                    @test x == 31
                    @test y == 52
                    @test x + y == 31 + 52
                end
                @assert 5 == 9
                @test x == 31
                @test 32 == 32
                @test y == 33
            end
            @testcase "Child XTest 3-2" begin
                @test 36 == 36
                @test 37 == 37
                @test fn_throws()
                @test 38 == 38
            end
            @testset "Child XTest 3-3" begin
                m = 31
                n = 52
                @test m == 31
                @test n == 52
                @test m + n == 31 + 52
            end
            @testcase "Child XTest 3-4" begin
                @test 41 == 41
                @test 42 == 42
                @test 43 == 43
            end
        end
    end
    @testset "XTest 4" begin
        @testset "Child XTest 4-1" begin
            @testcase "Child XTest 4-1-1" begin
                @test 4111 == 4111
                @test fn_throws()
                @test 4113 == 4113
            end
        end
        @testcase "Child XTest 6" begin
            @test 31 == 31
            @test 32 == 32
            @test 33 == 33
        end
    end
    include("subdir1/subdir1-unittest.jl")
end


text:
function probLineSearch(func::AbstractPLSFunction, x₀::Vector{T}, α₀::T;
        search_direction::Union{Nothing, Vector{T}}=nothing,
        σ²_f::Union{Nothing, T}=nothing,
        σ²_J::Union{Nothing, T, Array{T}}=nothing,
        pars::PLSConstantParams{T}=PLSConstantParams(T),
        history::PLSHistory=PLSHistory(α₀)) where T <: AbstractFloat

    # Setup general x0 (no projection onto line yet)
    x0 = evaluate(func, x₀; history=history)

    # Get line search direction (if not specified)
    (search_direction === nothing) && (search_direction = -x0.J)

    # Project onto line
    project_1d!(x0, search_direction)  # project gradient from Jacobian, and gradient variance.

    # User specified variances
    !(σ²_f === nothing) && (x0.σ²_f = σ²_f)
    !(σ²_J === nothing) && (x0.σ²_J = σ²_J)

    # Set-up line search parameters
    search = PLSSearchDefn(T, x₀, x0.y, α₀, search_direction, x0.σ²_f, x0.σ²_J)

    # Run line search
    x, gp = probLineSearch(func, x0, search, pars=pars, history=history)
    return x, gp, search, history
end

function probLineSearch(func::AbstractPLSFunction, x0::PLSEvaluation{T}, search::PLSSearchDefn{T};
        pars::PLSConstantParams{T}=PLSConstantParams(T), history::PLSHistory{T}=PLSHistory(T)) where T <: AbstractFloat

    @unpack α₀, search_direction, denom = search
    @unpack niter, c1, c2, wolfeThresh, verbosity = pars

    # project current evaluation onto the search direction specified.
    project_1d!(x0, search.search_direction)
    x0.t = 0

    # Need to learn normalisation factor from x0, then store for rest of proc.
    norm∇_0    = abs(x0.∇y)
    if !isapprox(norm∇_0, 1.0)
        search.denom = norm∇_0   # magnitude for rescaling: |y′(0)|
        search.x₀ = x0.x
        search.f₀ = x0.y         # offset for rescaling: y(0)
        normalise_for_linesearch!(x0, search)
    else
        !isapprox(search.denom, 1.0) && (@warn "Perhaps x0 has been pre-normalised; and norm factor in searchpars≈1.0. This may cause problems.")
    end

    σ²_f, σ²_∇ = normalised_fun_var(search)
    # @show sqrt(σ²_f)
    # @show sqrt(σ²_∇)
    gp = PLSPosterior(zeros(T, 1), zeros(T, 1), [x0.∇y], σ²_f, σ²_∇) # T, Y, ∇Y, ΣY, Σ∇

    iterates = [x0]

    tt = T(1.0)   # initial step size (rescaled)

    for i = 1:niter
        #======================================================================
                            Evaluate chosen point.
         =====================================================================#
        x = evaluate(func, tt, search; history=history)
        # @show x.x
        push!(iterates, x)

        gp = update(gp, x)
        # @show gp.Ts
        # @show gp.Y
        # @show gp.∇Y
        probWolfe(gp, tt, c1, c2)
        # -- check last evaluated point for acceptance ------------------------
        if probWolfe(gp, tt, c1, c2) > wolfeThresh   # are we done?
            # display("bark1")
            push!(history.msg, "found acceptable point. (1).")
            finalise(x, gp, search, pars, history);
            return x, gp # done  => x.x, x.y, x.J, x.∇_J, search.α₀, history.ewma_α
        end
        #======================================================================
            Has updating the GP led to an existing point qualifying?
         =====================================================================#
        # ----------------- find minimum of posterior mean --------------------
        minM, minT, minj  = min_mean(gp)   # of existing points
        min∇ = d1m(gp, gp.Ts[minj])

        # -------------- ∇ ≈ 0 with high probability? -------------------------
        if abs(min∇) < 1e-5 && Vd(gp, minT) < 1e-4 # nearly deterministic
            # display("bark2")
            push!(history.msg, "found a point with almost zero gradient. Stopping, although Wolfe conditions not guaranteed.")
            x = iterates[minj]
            finalise(x, gp, search, pars, history)
            return x, gp
        end

        # -------------- Wolfe Conditions satisfied? --------------------------
        Tsort  = sort(gp.Ts);
        # ignore first point
        wolfes_ix = 1 .+ findall(probWolfe.(gp, Tsort[2:end], c1, c2) .> wolfeThresh)

        if !isempty(wolfes_ix)
            # display("bark3")
            push!(history.msg, "found acceptable point. (2).")
            smallest_wolfe = argmin(m.(gp, Tsort[wolfes_ix]))
            tt = Tsort[wolfes_ix[smallest_wolfe]]  # choose point with smallest post. mean.
            minj = findfirst(gp.Ts .== tt)
            x = iterates[minj]
            finalise(x, gp, search, pars, history)
            return x, gp;
        end

        #======================================================================
            Choose next evaluation using EI/ProbWolfe via candidates
         =====================================================================#

        # -- CANDIDATES 1: minimal means between all evaluations.--------------

        # iterate through all `cells' (O[N]), check minimum mean locations.
        Tcand = zeros(T, 0); # positions of candidate points

        for cell = 1:length(gp)-1 # loop over cells
            _cellT_ϵ   = Tsort[cell] + T(1e-6) * (Tsort[cell+1] - Tsort[cell]);
            cell_minT  = cubicMinimum(gp, _cellT_ϵ);

            # add point to candidate list if minimum lies in between T(cel) and T(cel+1)
            if Tsort[cell] < cell_minT < Tsort[cell+1];
                push!(Tcand, cell_minT);

            # ~~PATHOLOGY~~
            #= If unable to find minimum in first cell *AND* the estimated ∇ > 0
               it is possible that (assuming the search direction is ok) the scale
               is wrong, so return a point v. close to original, and reduce α₀.=#
           elseif cell==1 && d1m(gp, 0) > 0
               display("bark4")
                push!(history.msg, "function seems very steep, reevaluating close to start.")
                tt = T(0.01) * (Tsort[cell] + Tsort[cell+1]);  # => Can lead to oscillation if T[1] is too far away.
                # x = iterates[1]
                # tt = T(0)
                x = evaluate(func, tt, search; history=history)
                # @show x.x
                finalise(x, gp, search, pars, history) #; force_α= search.α₀ / 2)
                # search.α₀ /= 10
                # history.ewma_α /= 5
                return x, gp  # done
            end
        end

        # -- CANDIDATES 2: one extrapolation step -----------------------------
        push!(Tcand, maximum(gp.Ts) + search.extrap_amt);
        # display(Tcand)
        # display("before")
        # -- DISCRIMINATION: EI and ProbWolfe ---------------------------------
        # Calculate some statistics for EI.
        μs = m(gp, Tcand)         # means of candidate points
        σs = __sqrt.(V(gp, Tcand))   # s.d. of candidate points
        # display(μs)
        # display(σs)
        # display(minM)
        # display(minT)
        # display(minj)
        EIcand = EI.(μs, σs, minM); # minM = η: best (lowest) GP mean of existing points

        PPcand = probWolfe.(gp, Tcand, c1, c2)
        idx_best = argmax(EIcand .* PPcand)

        # If best point is found using extrapolation, extend extrapolation for next iter.
        (idx_best == length(Tcand)) && (search.extrap_amt *= 2)
        # display(EIcand)
        # display(PPcand)
        # display("Chosen value: $idx_best")
        # Next point chosen for evaluation (next iter!):
        tt = Tcand[idx_best];

        # makePlot();
    end

    #======================================================================
        !!! Not found an acceptable candidate in the niter limit.
     =====================================================================#
    x = evaluate(func, tt, search, history=history)
    # @show x.x
    push!(iterates, x)
    gp = update(gp, x)

    # -- check last evaluated point for acceptance ------------------------
    if probWolfe(gp, tt, c1, c2) > wolfeThresh   # are we done?
        # display("bark5")
        push!(history.msg, "found acceptable point. (3).")
        finalise(x, gp, search, pars, history);
        return x, gp # done  => x.x, x.y, x.J, x.∇_J, search.α₀, history.ewma_α
    end

    # -- return point with lowest mean ----------------------------------------
    # Code of Mahsereci and Hennig exclude t=0 from this, but I'm not sure I buy this.
    minM, minT, minj  = min_mean(gp)
    # display("bark6")
    (verbosity > 0) && @warn "reached evaluation limit. Returning 'best' known point."
    push!(history.msg, "reached evaluation limit. Returning 'best' known point.")
    x = iterates[minj]
    finalise(x, gp, search, pars, history);
    return x, gp

    # tt = Tsort[argmin(m.(gp, Tsort[wolfes_ix]))]  # choose point with smallest post. mean.
    # minj = findfirst(gp.Ts .== tt)
    # x = iterates[minj]
    # finalise(x, gp, search, pars, history)
    # makePlot();
end

function finalise(x::PLSEvaluation{T}, post::PLSPosterior{T}, searchpars::PLSSearchDefn{T},
        pars::PLSConstantParams, history::PLSHistory{T}; force_α::Union{T,Nothing}=nothing) where T <: AbstractFloat
    @unpack x₀ ,α₀, f₀, denom = searchpars
    @unpack αGrowth, verbosity, c1, c2 = pars

    # Transform `x` back to original frame of reference / undo [0,1] normalisation.
    invert_normalise_for_linesearch!(x, searchpars)

    # save step size
    α = x.t*α₀
    push!(history.αs, α)

    # next initial step size is αGrowth (def=1.3x) larger than last accepted step size
    searchpars.α₀ = something(force_α, α * αGrowth)
    check_α = (force_α === nothing)  # if new alpha not forced by input arg

    # running average for reset in case the step size becomes very small
    # this is a safeguard
    gamma = 0.95;
    history.ewma_α = gamma*history.ewma_α + (1-gamma)*something(force_α, α);

    # reset NEXT initial step size to average step size if accepted step
    # size is 100 times smaller or larger than average step size
    if check_α && !(1e-2 * history.ewma_α < α < 1e2 * history.ewma_α)
        push!(history.msg, "v large/small value of α, resetting alpha0.")
        (verbosity > 0) && @warn "v large/small value of α, resetting alpha0. To avoid cycling, ensure passing PLSHistory object through."
        searchpars.α₀ = history.ewma_α # reset step size
    end

    # reset search parameters for next iteration
    # --- any extrapolation should be baked into α₀
    searchpars.extrap_amt = 1

    # --- update noise according to x
    searchpars.σ²_f = x.σ²_f
    searchpars.σ²_J = x.σ²_J

    # Calculate acceptance statistics
    minM, _minT, _minj  = min_mean(post)   # of existing points
    push!(history.EI, EI(m(post, x.t), __sqrt(V(post, x.t)), minM))
    push!(history.Wolfe, probWolfe(post, x.t, c1, c2))
    history.recording = false
end

max_stars_repo_name: ornithos/ProbLineSearch


text:
using TableTestSets
using Test

@testset "TableTestSets.jl" begin
    results = @testset TableTestSet "outer" begin
        @testset "inner 1" begin
            @test 1 // 1 == 1
            @test false
        end

        @testset "inner 2" begin
            @test true
            @test error()
        end
    end

    io = IOBuffer()
    html_table(io, results)
    str = String(take!(io))

    @test str == "
\n
\ntestset	\npass	\nfail	\nerror	\nbroken	\ntotal	\n
\n
outer	\n2	\n1	\n1	\n0	\n4	\n

  inner 1	\n1	\n1	\n0	\n0	\n2	\n

  inner 2	\n1	\n0	\n1	\n0	\n2	\n

\n"
end


text:
"""

Run Stan's diagnose binary on a model.

$(SIGNATURES)

### Required arguments
```julia
* `model`                    : SampleModel
```
"""
function diagnose(model::SampleModel)
  #local csvfile
  n_chains = model.n_chains[1]

  samplefiles = String[]
  for i in 1:n_chains
    push!(samplefiles, "$(model.output_base)_chain_$(i).csv")
  end
  try
    pstring = joinpath("$(model.cmdstan_home)", "bin", "diagnose")
    if Sys.iswindows()
      pstring = pstring * ".exe"
    end
    cmd = `$(pstring) $(StanBase.par(samplefiles))`
    resfile = open(cmd; read=true);
    print(read(resfile, String))
  catch e
    println(e)
  end

  return
end

export
    diagnose


text:
# see NOTE in reference_ledoitwolf.jl
# this comes from appendix D in http://www.econ.uzh.ch/static/wp/econwp264.pdf
# adapted to Julia and checked for match with Octave.
# NOTE: it was modified to make the case p < n with centering work:
# - lambda[max(1, p-n+2):p] instead of lambda[max(1, p-n+1):p]
# - vcat(dtilde0 * ones(p-n+1,1), dtilde1); instead of vcat(dtilde0 * ones(p-n+1,1), dtilde1);

using LinearAlgebra
using Statistics

function matlab_ledoitwolf_analytical_shrinkage(X)
    n, p = size(X);
    @assert n ≥ 12 "there must be more than 12 samples"
    sample = cov(SimpleCovariance(), X); # note centering is applied here
    lambda, u = eigen(sample)
    perm = sortperm(lambda)
    lambda = lambda[perm]

    u = u[:, perm]
    η = ifelse(p<n, n, n-1)
    lambda = lambda[max(1, p-η+1):p]
    L = repeat(lambda, outer=(1, min(p, η)))

    h = η^(-1/3)
    H = h*L';
    x = (L-L')./H;

    ftilde = (3/4/sqrt(5))*mean(max.(1 .- x.^2 ./5,0)./H, dims=2);
    Hftemp = (-3/10/pi)*x + (3/4/sqrt(5)/pi) * (1 .- x.^2 ./5) .* log.(abs.((sqrt(5) .- x) ./ (sqrt(5) .+ x)));
    Hftemp[abs.(x).==sqrt(5)] .= (-3/10/pi)*x[abs.(x).==sqrt(5)];
    Hftilde = mean(Hftemp./H, dims=2);

    if p <= η
        dtilde = lambda ./ ((pi*(p/n)*lambda.*ftilde).^2 .+ (1 .- (p/n) .- pi*(p/n)*lambda.*Hftilde).^2);
    else
        Hftilde0 = (1/pi)*(3/10/h^2 + 3/4/sqrt(5)/h*(1-1/5/h^2) * log((1+sqrt(5)*h)/(1-sqrt(5)*h))) * mean(1 ./ lambda);
        dtilde0 = 1/(pi*(p-η)/η*Hftilde0);
        dtilde1 = lambda ./ (pi^2*lambda.^2 .* (ftilde.^2 .+ Hftilde.^2));
        dtilde = vcat(dtilde0 * ones(p-η,1), dtilde1);
    end

    u * (dtilde[:] .* u')
end


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

module Consts

import ..LibGit2: version

const HEAD_FILE  = "HEAD"
const FETCH_HEAD  = "FETCH_HEAD"
const REMOTE_ORIGIN = "origin"

# objs
@enum(OBJECT,
      OBJ_ANY    = -2,
      OBJ_BAD    = -1,
      OBJ_COMMIT = 1,
      OBJ_TREE   = 2,
      OBJ_BLOB   = 3,
      OBJ_TAG    = 4)

#revwalk
const SORT_NONE        = Cint(0)
const SORT_TOPOLOGICAL = Cint(1 << 0)
const SORT_TIME        = Cint(1 << 1)
const SORT_REVERSE     = Cint(1 << 2)

# refs
const REF_INVALID  = Cint(0)
const REF_OID      = Cint(1)
const REF_SYMBOLIC = Cint(2)
const REF_LISTALL  = REF_OID | REF_SYMBOLIC

# blame
const BLAME_NORMAL                          = Cuint(0)
const BLAME_TRACK_COPIES_SAME_FILE          = Cuint(1 << 0)
const BLAME_TRACK_COPIES_SAME_COMMIT_MOVES  = Cuint(1 << 1)
const BLAME_TRACK_COPIES_SAME_COMMIT_COPIES = Cuint(1 << 2)
const BLAME_TRACK_COPIES_ANY_COMMIT_COPIES  = Cuint(1 << 3)
const BLAME_FIRST_PARENT                    = Cuint(1 << 4)

# checkout
const CHECKOUT_NONE                    = Cuint(0)
const CHECKOUT_SAFE                    = Cuint(1 << 0)
const CHECKOUT_FORCE                   = Cuint(1 << 1)
const CHECKOUT_RECREATE_MISSING        = Cuint(1 << 2)
const CHECKOUT_ALLOW_CONFLICTS         = Cuint(1 << 4)
const CHECKOUT_REMOVE_UNTRACKED        = Cuint(1 << 5)
const CHECKOUT_REMOVE_IGNORED          = Cuint(1 << 6)
const CHECKOUT_UPDATE_ONLY             = Cuint(1 << 7)
const CHECKOUT_DONT_UPDATE_INDEX       = Cuint(1 << 8)
const CHECKOUT_NO_REFRESH              = Cuint(1 << 9)
const CHECKOUT_SKIP_UNMERGED           = Cuint(1 << 10)
const CHECKOUT_USE_OURS                = Cuint(1 << 11)
const CHECKOUT_USE_THEIRS              = Cuint(1 << 12)
const CHECKOUT_DISABLE_PATHSPEC_MATCH  = Cuint(1 << 13)
const CHECKOUT_SKIP_LOCKED_DIRECTORIES = Cuint(1 << 18)
const CHECKOUT_DONT_OVERWRITE_IGNORED  = Cuint(1 << 19)
const CHECKOUT_CONFLICT_STYLE_MERGE    = Cuint(1 << 20)
const CHECKOUT_CONFLICT_STYLE_DIFF3    = Cuint(1 << 21)
const CHECKOUT_DONT_REMOVE_EXISTING    = Cuint(1 << 22)

const CHECKOUT_UPDATE_SUBMODULES       = Cuint(1 << 16)
const CHECKOUT_UPDATE_SUBMODULES_IF_CHANGED = Cuint(1 << 17)

const CHECKOUT_NOTIFY_NONE      = Cuint(0)
const CHECKOUT_NOTIFY_CONFLICT  = Cuint(1 << 0)
const CHECKOUT_NOTIFY_DIRTY     = Cuint(1 << 1)
const CHECKOUT_NOTIFY_UPDATED   = Cuint(1 << 2)
const CHECKOUT_NOTIFY_UNTRACKED = Cuint(1 << 3)
const CHECKOUT_NOTIFY_IGNORED   = Cuint(1 << 4)
const CHECKOUT_NOTIFY_ALL       = 0x0FFFF

# diff
const DIFF_OPTIONS_VERSION = Cuint(1)

const DIFF_NORMAL                     = Cuint(0)
const DIFF_REVERSE                    = Cuint(1 << 0)
const DIFF_INCLUDE_IGNORED            = Cuint(1 << 1)
const DIFF_RECURSE_IGNORED_DIRS       = Cuint(1 << 2)
const DIFF_INCLUDE_UNTRACKED          = Cuint(1 << 3)
const DIFF_RECURSE_UNTRACKED_DIRS     = Cuint(1 << 4)
const DIFF_INCLUDE_UNMODIFIED         = Cuint(1 << 5)
const DIFF_INCLUDE_TYPECHANGE         = Cuint(1 << 6)
const DIFF_INCLUDE_TYPECHANGE_TREES   = Cuint(1 << 7)
const DIFF_IGNORE_FILEMODE            = Cuint(1 << 8)
const DIFF_IGNORE_SUBMODULES          = Cuint(1 << 9)
const DIFF_IGNORE_CASE                = Cuint(1 << 10)
const DIFF_DISABLE_PATHSPEC_MATCH     = Cuint(1 << 12)
const DIFF_SKIP_BINARY_CHECK          = Cuint(1 << 13)
const DIFF_ENABLE_FAST_UNTRACKED_DIRS = Cuint(1 << 14)

const DIFF_FORCE_TEXT               = Cuint(1 << 20)
const DIFF_FORCE_BINARY             = Cuint(1 << 21)
const DIFF_IGNORE_WHITESPACE        = Cuint(1 << 22)
const DIFF_IGNORE_WHITESPACE_CHANGE = Cuint(1 << 23)
const DIFF_IGNORE_WHITESPACE_EOL    = Cuint(1 << 24)
const DIFF_SHOW_UNTRACKED_CONTENT   = Cuint(1 << 25)
const DIFF_SHOW_UNMODIFIED          = Cuint(1 << 26)
const DIFF_PATIENCE                 = Cuint(1 << 28)
const DIFF_MINIMAL                  = Cuint(1 << 29)

const DIFF_FLAG_BINARY     = Cuint(1 << 0)
const DIFF_FLAG_NOT_BINARY = Cuint(1 << 1)
const DIFF_FLAG_VALID_OID  = Cuint(1 << 2)

const DIFF_FORMAT_PATCH        = Cuint(1)
const DIFF_FORMAT_PATCH_HEADER = Cuint(2)
const DIFF_FORMAT_RAW          = Cuint(3)
const DIFF_FORMAT_NAME_ONLY    = Cuint(4)
const DIFF_FORMAT_NAME_STATUS  = Cuint(5)

@enum(DELTA_STATUS, DELTA_UNMODIFIED = Cint(0),
                    DELTA_ADDED      = Cint(1),
                    DELTA_DELETED    = Cint(2),
                    DELTA_MODIFIED   = Cint(3),
                    DELTA_RENAMED    = Cint(4),
                    DELTA_COPIED     = Cint(5),
                    DELTA_IGNORED    = Cint(6),
                    DELTA_UNTRACKED  = Cint(7),
                    DELTA_TYPECHANGE = Cint(8))

# index
const IDXENTRY_NAMEMASK   = (0x0fff)
const IDXENTRY_STAGEMASK  = (0x3000)
const IDXENTRY_EXTENDED   = (0x4000)
const IDXENTRY_VALID      = (0x8000)
const IDXENTRY_STAGESHIFT = Cint(12)

const IDXENTRY_UPDATE            = Cint(1 << 0)
const IDXENTRY_REMOVE            = Cint(1 << 1)
const IDXENTRY_UPTODATE          = Cint(1 << 2)
const IDXENTRY_ADDED             = Cint(1 << 3)

const IDXENTRY_HASHED            = Cint(1 << 4)
const IDXENTRY_UNHASHED          = Cint(1 << 5)
const IDXENTRY_WT_REMOVE         = Cint(1 << 6)
const IDXENTRY_CONFLICTED        = Cint(1 << 7)

const IDXENTRY_UNPACKED          = Cint(1 << 8)
const IDXENTRY_NEW_SKIP_WORKTREE = Cint(1 << 9)

const INDEXCAP_IGNORE_CASE = Cuint(1)
const INDEXCAP_NO_FILEMODE = Cuint(2)
const INDEXCAP_NO_SYMLINKS = Cuint(4)
const INDEXCAP_FROM_OWNER  = ~Cuint(0)

const INDEX_ADD_DEFAULT                = Cuint(0)
const INDEX_ADD_FORCE                  = Cuint(1 << 0)
const INDEX_ADD_DISABLE_PATHSPEC_MATCH = Cuint(1 << 1)
const INDEX_ADD_CHECK_PATHSPEC         = Cuint(1 << 2)

const INDEX_STAGE_ANY = Cint(-1)

# merge
""" Option flags for git merge.
* `MERGE_FIND_RENAMES`: detect if a file has been renamed between the common
  ancestor and the "ours" or "theirs" side of the merge. Allows merges where
  a file has been renamed.
* `MERGE_FAIL_ON_CONFLICT`: exit immediately if a conflict is found rather
  than trying to resolve it.
* `MERGE_SKIP_REUC`: do not write the REUC extension on the index resulting
  from the merge.
* `MERGE_NO_RECURSIVE`: if the commits being merged have multiple merge bases,
  use the first one, rather than trying to recursively merge the bases.
"""
@enum(GIT_MERGE, MERGE_FIND_RENAMES     = 1 << 0,
                 MERGE_FAIL_ON_CONFLICT = 1 << 1,
                 MERGE_SKIP_REUC        = 1 << 2,
                 MERGE_NO_RECURSIVE     = 1 << 3)

@enum(GIT_MERGE_FILE, MERGE_FILE_DEFAULT                  = 0,       # Defaults
                      MERGE_FILE_STYLE_MERGE              = 1 << 0,  # Create standard conflicted merge files
                      MERGE_FILE_STYLE_DIFF3              = 1 << 1,  # Create diff3-style files
                      MERGE_FILE_SIMPLIFY_ALNUM           = 1 << 2,  # Condense non-alphanumeric regions for simplified diff file
                      MERGE_FILE_IGNORE_WHITESPACE        = 1 << 3,  # Ignore all whitespace
                      MERGE_FILE_IGNORE_WHITESPACE_CHANGE = 1 << 4,  # Ignore changes in amount of whitespace
                      MERGE_FILE_IGNORE_WHITESPACE_EOL    = 1 << 5,  # Ignore whitespace at end of line
                      MERGE_FILE_DIFF_PATIENCE            = 1 << 6,  # Use the "patience diff" algorithm
                      MERGE_FILE_DIFF_MINIMAL             = 1 << 7)  # Take extra time to find minimal diff
""" Option flags for git merge file favoritism.
  * `MERGE_FILE_FAVOR_NORMAL`: if both sides of the merge have changes to a section,
    make a note of the conflict in the index which `git checkout` will use to create
    a merge file, which the user can then reference to resolve the conflicts. This is
    the default.
  * `MERGE_FILE_FAVOR_OURS`: if both sides of the merge have changes to a section,
    use the version in the "ours" side of the merge in the index.
  * `MERGE_FILE_FAVOR_THEIRS`: if both sides of the merge have changes to a section,
    use the version in the "theirs" side of the merge in the index.
  * `MERGE_FILE_FAVOR_UNION`: if both sides of the merge have changes to a section,
    include each unique line from both sides in the file which is put into the index.
"""
@enum(GIT_MERGE_FILE_FAVOR, MERGE_FILE_FAVOR_NORMAL = 0,
                            MERGE_FILE_FAVOR_OURS   = 1,
                            MERGE_FILE_FAVOR_THEIRS = 2,
                            MERGE_FILE_FAVOR_UNION  = 3)
""" The user's instructions for how to perform a possible merge.
* `MERGE_PREFERENCE_NONE`: the user has no preference.
* `MERGE_PREFERENCE_NO_FASTFORWARD`: do not allow any fast-forward merges.
* `MERGE_PREFERENCE_FASTFORWARD_ONLY`: allow only fast-forward merges and no
  other type (which may introduce conflicts).
"""
@enum(GIT_MERGE_PREFERENCE, MERGE_PREFERENCE_NONE             = 0,
                            MERGE_PREFERENCE_NO_FASTFORWARD   = 1,
                            MERGE_PREFERENCE_FASTFORWARD_ONLY = 2)
""" Result of analysis on merge possibilities.
* `MERGE_ANALYSIS_NONE`: it is not possible to merge the elements of the input commits.
* `MERGE_ANALYSIS_NORMAL`: a regular merge, when HEAD and the commits that the
  user wishes to merge have all diverged from a common ancestor. In this case the
  changes have to be resolved and conflicts may occur.
* `MERGE_ANALYSIS_UP_TO_DATE`: all the input commits the user wishes to merge can
  be reached from HEAD, so no merge needs to be performed.
* `MERGE_ANALYSIS_FASTFORWARD`: the input commit is a descendant of HEAD and so no
  merge needs to be performed - instead, the user can simply checkout the
  input commit(s).
* `MERGE_ANALYSIS_UNBORN`: the HEAD of the repository refers to a commit which does not
  exist. It is not possible to merge, but it may be possible to checkout the input
  commits.
"""
@enum(GIT_MERGE_ANALYSIS, MERGE_ANALYSIS_NONE        = 0,
                          MERGE_ANALYSIS_NORMAL      = 1 << 0,
                          MERGE_ANALYSIS_UP_TO_DATE  = 1 << 1,
                          MERGE_ANALYSIS_FASTFORWARD = 1 << 2,
                          MERGE_ANALYSIS_UNBORN      = 1 << 3)

# reset
const RESET_SOFT  = Cint(1) # Move the head to the given commit
const RESET_MIXED = Cint(2) # SOFT plus reset index to the commit
const RESET_HARD  = Cint(3) # MIXED plus changes in working tree discarded

# rebase
""" Options for what rebase operation is currently being performed on a commit.
* `REBASE_OPERATION_PICK`: cherry-pick the commit in question.
* `REBASE_OPERATION_REWORD`: cherry-pick the commit in question, but rewrite its
  message using the prompt.
* `REBASE_OPERATION_EDIT`: cherry-pick the commit in question, but allow the user
  to edit the commit's contents and its message.
* `REBASE_OPERATION_SQUASH`: squash the commit in question into the previous commit.
  The commit messages of the two commits will be merged.
* `REBASE_OPERATION_FIXUP`: squash the commit in question into the previous commit.
  Only the commit message of the previous commit will be used.
* `REBASE_OPERATION_EXEC`: do not cherry-pick a commit. Run a command and continue if
  the command exits successfully.
"""
@enum(GIT_REBASE_OPERATION, REBASE_OPERATION_PICK   = Cint(0),
                            REBASE_OPERATION_REWORD = Cint(1),
                            REBASE_OPERATION_EDIT   = Cint(2),
                            REBASE_OPERATION_SQUASH = Cint(3),
                            REBASE_OPERATION_FIXUP  = Cint(4),
                            REBASE_OPERATION_EXEC   = Cint(5))

# fetch_prune
const FETCH_PRUNE_UNSPECIFIED = Cint(0)
const FETCH_PRUNE             = Cint(1)
const FETCH_NO_PRUNE          = Cint(2)

# remote_autotag
const REMOTE_DOWNLOAD_TAGS_UNSPECIFIED = Cint(0)
const REMOTE_DOWNLOAD_TAGS_AUTO        = Cint(1)
const REMOTE_DOWNLOAD_TAGS_NONE        = Cint(2)
const REMOTE_DOWNLOAD_TAGS_ALL         = Cint(3)

# clone
const CLONE_LOCAL_AUTO     = Cint(0)
const CLONE_LOCAL          = Cint(1)
const CLONE_NO_LOCAL       = Cint(2)
const CLONE_LOCAL_NO_LINKS = Cint(3)

# describe
const DESCRIBE_DEFAULT = Cuint(0)
const DESCRIBE_TAGS    = Cuint(1 << 0)
const DESCRIBE_ALL     = Cuint(1 << 1)

# status
const STATUS_CURRENT          = Cuint(0)
const STATUS_INDEX_NEW        = Cuint(1 << 0)
const STATUS_INDEX_MODIFIED   = Cuint(1 << 1)
const STATUS_INDEX_DELETED    = Cuint(1 << 2)
const STATUS_INDEX_RENAMED    = Cuint(1 << 3)
const STATUS_INDEX_TYPECHANGE = Cuint(1 << 4)
const STATUS_WT_NEW           = Cuint(1 << 7)
const STATUS_WT_MODIFIED      = Cuint(1 << 8)
const STATUS_WT_DELETED       = Cuint(1 << 9)
const STATUS_WT_TYPECHANGE    = Cuint(1 << 10)
const STATUS_WT_RENAMED       = Cuint(1 << 11)
const STATUS_WT_UNREADABLE    = Cuint(1 << 12)
const STATUS_IGNORED          = Cuint(1 << 14)
const STATUS_CONFLICTED       = Cuint(1 << 15)

# status show
const STATUS_SHOW_INDEX_AND_WORKDIR = Cint(0)
const STATUS_SHOW_INDEX_ONLY        = Cint(1)
const STATUS_SHOW_WORKDIR_ONLY      = Cint(2)

# status options
const STATUS_OPT_INCLUDE_UNTRACKED                = Cuint(1 << 0)
const STATUS_OPT_INCLUDE_IGNORED                  = Cuint(1 << 1)
const STATUS_OPT_INCLUDE_UNMODIFIED               = Cuint(1 << 2)
const STATUS_OPT_EXCLUDE_SUBMODULES               = Cuint(1 << 3)
const STATUS_OPT_RECURSE_UNTRACKED_DIRS           = Cuint(1 << 4)
const STATUS_OPT_DISABLE_PATHSPEC_MATCH           = Cuint(1 << 5)
const STATUS_OPT_RECURSE_IGNORED_DIRS             = Cuint(1 << 6)
const STATUS_OPT_RENAMES_HEAD_TO_INDEX            = Cuint(1 << 7)
const STATUS_OPT_RENAMES_INDEX_TO_WORKDIR         = Cuint(1 << 8)
const STATUS_OPT_SORT_CASE_SENSITIVELY            = Cuint(1 << 9)
const STATUS_OPT_SORT_CASE_INSENSITIVELY          = Cuint(1 << 10)
const STATUS_OPT_RENAMES_FROM_REWRITES            = Cuint(1 << 11)
const STATUS_OPT_NO_REFRESH                       = Cuint(1 << 12)
const STATUS_OPT_UPDATE_INDEX                     = Cuint(1 << 13)
const STATUS_OPT_INCLUDE_UNREADABLE               = Cuint(1 << 14)
const STATUS_OPT_INCLUDE_UNREADABLE_AS_UNTRACKED  = Cuint(1 << 15)

@enum(GIT_SUBMODULE_IGNORE, SUBMODULE_IGNORE_UNSPECIFIED  = -1, # use the submodule's configuration
                            SUBMODULE_IGNORE_NONE         = 1,  # any change or untracked == dirty
                            SUBMODULE_IGNORE_UNTRACKED    = 2,  # dirty if tracked files change
                            SUBMODULE_IGNORE_DIRTY        = 3,  # only dirty if HEAD moved
                            SUBMODULE_IGNORE_ALL          = 4)  # never dirty

"""
Option flags for `GitRepo`.

* `REPOSITORY_OPEN_NO_SEARCH` - Only open the repository if it can be immediately found in the `path`.  Do not walk up from the `path` looking at parent directories.
* `REPOSITORY_OPEN_CROSS_FS` - Unless this flag is set, open will not continue searching across filesystem boundaries. (E.g. Searching in a user's home directory `/home/user/source/` will not return `/.git/` as the found repo if `/` is a different filesystem than `/home`.)
* `REPOSITORY_OPEN_BARE` - Open repository as a bare repo regardless of core.bare config, and defer loading config file for faster setup.
"""
@enum(GIT_REPOSITORY_OPEN, REPOSITORY_OPEN_DEFAULT   = 0,
                           REPOSITORY_OPEN_NO_SEARCH = 1<<0,
                           REPOSITORY_OPEN_CROSS_FS  = 1<<1,
                           REPOSITORY_OPEN_BARE      = 1<<2)

@enum(GIT_BRANCH, BRANCH_LOCAL = 1, BRANCH_REMOTE = 2)

@enum(GIT_FILEMODE, FILEMODE_UNREADABLE          = 0o000000,
                    FILEMODE_TREE                = 0o040000,
                    FILEMODE_BLOB                = 0o100644,
                    FILEMODE_BLOB_EXECUTABLE     = 0o100755,
                    FILEMODE_LINK                = 0o120000,
                    FILEMODE_COMMIT              = 0o160000)

@enum(GIT_CREDTYPE, CREDTYPE_USERPASS_PLAINTEXT = Cuint(1 << 0),
                    CREDTYPE_SSH_KEY            = Cuint(1 << 1),
                    CREDTYPE_SSH_CUSTOM         = Cuint(1 << 2),
                    CREDTYPE_DEFAULT            = Cuint(1 << 3),
                    CREDTYPE_SSH_INTERACTIVE    = Cuint(1 << 4),
                    CREDTYPE_USERNAME           = Cuint(1 << 5),
                    CREDTYPE_SSH_MEMORY         = Cuint(1 << 6))

@enum(GIT_FEATURE, FEATURE_THREADS = Cuint(1 << 0),
                   FEATURE_HTTPS   = Cuint(1 << 1),
                   FEATURE_SSH     = Cuint(1 << 2),
                   FEATURE_NSEC    = Cuint(1 << 3))

if version() >= v"0.24.0"
    @doc """
    Priority level of a config file.

    These priority levels correspond to the natural escalation logic (from higher to lower) when searching for config entries in git.

    * `CONFIG_LEVEL_DEFAULT` - Open the global, XDG and system configuration files if any available.
    * `CONFIG_LEVEL_PROGRAMDATA` - System-wide on Windows, for compatibility with portable git
    * `CONFIG_LEVEL_SYSTEM` - System-wide configuration file; `/etc/gitconfig` on Linux systems
    * `CONFIG_LEVEL_XDG` - XDG compatible configuration file; typically `~/.config/git/config`
    * `CONFIG_LEVEL_GLOBAL` - User-specific configuration file (also called Global configuration file); typically `~/.gitconfig`
    * `CONFIG_LEVEL_LOCAL` - Repository specific configuration file; `\$WORK_DIR/.git/config` on non-bare repos
    * `CONFIG_LEVEL_APP` - Application specific configuration file; freely defined by applications
    * `CONFIG_HIGHEST_LEVEL` - Represents the highest level available config file (i.e. the most specific config file available that actually is loaded)
    """ ->
    @enum(GIT_CONFIG, CONFIG_LEVEL_DEFAULT     = 0,
                      CONFIG_LEVEL_PROGRAMDATA = 1,
                      CONFIG_LEVEL_SYSTEM      = 2,
                      CONFIG_LEVEL_XDG         = 3,
                      CONFIG_LEVEL_GLOBAL      = 4,
                      CONFIG_LEVEL_LOCAL       = 5,
                      CONFIG_LEVEL_APP         = 6,
                      CONFIG_HIGHEST_LEVEL     =-1)
else
    @doc """
    Priority level of a config file.

    These priority levels correspond to the natural escalation logic (from higher to lower) when searching for config entries in git.

    * `CONFIG_LEVEL_DEFAULT` - Open the global, XDG and system configuration files if any available.
    * `CONFIG_LEVEL_SYSTEM` - System-wide configuration file; `/etc/gitconfig` on Linux systems
    * `CONFIG_LEVEL_XDG` - XDG compatible configuration file; typically `~/.config/git/config`
    * `CONFIG_LEVEL_GLOBAL` - User-specific configuration file (also called Global configuration file); typically `~/.gitconfig`
    * `CONFIG_LEVEL_LOCAL` - Repository specific configuration file; `\$WORK_DIR/.git/config` on non-bare repos
    * `CONFIG_LEVEL_APP` - Application specific configuration file; freely defined by applications
    * `CONFIG_HIGHEST_LEVEL` - Represents the highest level available config file (i.e. the most specific config file available that actually is loaded)
    """ ->
    @enum(GIT_CONFIG, CONFIG_LEVEL_DEFAULT     = 0,
                      CONFIG_LEVEL_SYSTEM      = 1,
                      CONFIG_LEVEL_XDG         = 2,
                      CONFIG_LEVEL_GLOBAL      = 3,
                      CONFIG_LEVEL_LOCAL       = 4,
                      CONFIG_LEVEL_APP         = 5,
                      CONFIG_HIGHEST_LEVEL     =-1)
end

"""
Global library options.

These are used to select which global option to set or get and are used in `git_libgit2_opts()`.
"""
@enum(GIT_OPT, GET_MWINDOW_SIZE         = 0,
               SET_MWINDOW_SIZE         = 1,
               GET_MWINDOW_MAPPED_LIMIT = 2,
               SET_MWINDOW_MAPPED_LIMIT = 3,
               GET_SEARCH_PATH          = 4,
               SET_SEARCH_PATH          = 5,
               SET_CACHE_OBJECT_LIMIT   = 6,
               SET_CACHE_MAX_SIZE       = 7,
               ENABLE_CACHING           = 8,
               GET_CACHED_MEMORY        = 9,
               GET_TEMPLATE_PATH        = 10,
               SET_TEMPLATE_PATH        = 11,
               SET_SSL_CERT_LOCATIONS   = 12)

"""
Option flags for `GitProxy`.

* `PROXY_NONE`: do not attempt the connection through a proxy.
* `PROXY_AUTO`: attempt to figure out the proxy configuration from the git configuration.
* `PROXY_SPECIFIED`: connect using the URL given in the `url` field of this struct.
"""
@enum(GIT_PROXY, PROXY_NONE,
                 PROXY_AUTO,
                 PROXY_SPECIFIED)

end

max_stars_repo_name: BenChung/Julia-BTI

max_issues_repo_name: BenChung/Julia-BTI

max_forks_repo_name: BenChung/Julia-BTI


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.6.2
#     language: julia
#     name: julia-1.6
# ---

# %% [markdown]
# https://discourse.julialang.org/t/concrete-example-0805-how-to-prevent-using-global-variable-in-the-module/65917

# %%
module Mixture

export musigma

mutable struct Mean_covar
    mu::Array{Float64,2}
    sigma::Array{Float64,2}
    w::Float64
end

global const musigma = Array{Mean_covar,1}()

function mean_covar_init(kmix::Int64,dim_p::Int64
                        ,weight::Array{Float64,1}
                        ,sigma::Array{Float64,3}
                        ,mu::Array{Float64,2})
    @assert length(weight) == kmix
    @assert size(sigma) == (kmix,dim_p,dim_p)
    @assert size(mu) == (kmix,dim_p)
    resize!(musigma, kmix)
    for k in 1:kmix
        musigma[k] = Mean_covar(zeros(dim_p,1),zeros(dim_p,dim_p),0.0)
        musigma[k].mu[1,1] = mu[k,1]
        musigma[k].mu[2,1] = mu[k,2]
        musigma[k].sigma[1,1] = sigma[k,1,1]
        musigma[k].sigma[2,2] = sigma[k,2,2]
        musigma[k].w = weight[k]
    end
    return nothing
end

end

kmix = 5
dim_p = 3
weight = rand(kmix)
sigma = rand(kmix, dim_p, dim_p)
mu = rand(kmix, dim_p)

using .Mixture
Mixture.mean_covar_init(kmix, dim_p, weight, sigma, mu)
@show typeof(musigma) size(musigma)
@show musigma[1].mu musigma[1].sigma musigma[1].w;

# %%
module O

export musigma

struct Mean_covar{T}
    mu::Vector{T}
    sigma::Matrix{T}
    w::Array{T,0}
end

function mean_covar_init(kmix, dim_p,
        weight::AbstractVector{T},
        sigma::AbstractArray{T,3},
        mu::AbstractMatrix{T}) where T
    @assert size(weight) == (kmix,)
    @assert size(sigma) == (kmix, dim_p, dim_p)
    @assert size(mu) == (kmix, dim_p)
    musigma = Vector{Mean_covar{T}}(undef, kmix)
    for k in 1:kmix
        musigma[k] = Mean_covar(zeros(T, dim_p), zeros(T, dim_p, dim_p), fill(zero(T)))
        musigma[k].mu[1,1] = mu[k,1]
        musigma[k].mu[2,1] = mu[k,2]
        musigma[k].sigma[1,1] = sigma[k,1,1]
        musigma[k].sigma[2,2] = sigma[k,2,2]
        musigma[k].w[] = weight[k]
    end
    musigma
end

end

musigma64 = O.mean_covar_init(kmix, dim_p, weight, sigma, mu)
@show typeof(musigma64) size(musigma64)
@show musigma64[1].mu musigma64[1].sigma musigma64[1].w;

# %%
weight32 = rand(Float32, kmix)
sigma32 = rand(Float32, kmix, dim_p, dim_p)
mu32 = rand(Float32, kmix, dim_p)

musigma32 = O.mean_covar_init(kmix, dim_p, weight32, sigma32, mu32)
@show typeof(musigma32) size(musigma32)
@show musigma32[1].mu musigma32[1].sigma musigma32[1].w;

# %%
module P

export musigma

struct Mean_covar{M,S,W}
    mu::M
    sigma::S
    w::W
end

function mean_covar_init(kmix, dim_p,
        weight::AbstractVector{T},
        sigma::AbstractArray{T,3},
        mu::AbstractMatrix{T}) where T
    @assert size(weight) == (kmix,)
    @assert size(sigma) == (kmix, dim_p, dim_p)
    @assert size(mu) == (kmix, dim_p)
    M = typeof(similar(mu, (0,)))
    S = typeof(similar(sigma, (0, 0)))
    W = typeof(similar(weight, ()))
    musigma = Vector{Mean_covar{M,S,W}}(undef, kmix)
    for k in 1:kmix
        musigma[k] = Mean_covar(zeros(T, dim_p), zeros(T, dim_p, dim_p), fill(zero(T)))
        musigma[k].mu[1,1] = mu[k,1]
        musigma[k].mu[2,1] = mu[k,2]
        musigma[k].sigma[1,1] = sigma[k,1,1]
        musigma[k].sigma[2,2] = sigma[k,2,2]
        musigma[k].w[] = weight[k]
    end
    musigma
end

end

musigma64 = P.mean_covar_init(kmix, dim_p, weight, sigma, mu)
@show typeof(musigma64) size(musigma64)
@show musigma64[1].mu musigma64[1].sigma musigma64[1].w;

# %%
musigma32 = O.mean_covar_init(kmix, dim_p, weight32, sigma32, mu32)
@show typeof(musigma32) size(musigma32)
@show musigma32[1].mu musigma32[1].sigma musigma32[1].w;

# %%

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
using OpenQuantumSystems
import QuantumOpticsBase, LinearAlgebra, OrdinaryDiffEq, QuadGK, DelayDiffEq

function master_ansatz(
    rho0::T,
    tspan::Array,
    p;
    reltol::Float64 = 1.0e-12,
    abstol::Float64 = 1.0e-12,
    alg::Any = DelayDiffEq.MethodOfSteps(DelayDiffEq.Vern6()),
    fout::Union{Function,Nothing} = nothing,
    kwargs...,
) where {B<:Basis,T<:Operator{B,B}}
    history_fun(p, t) = T(rho0.basis_l, rho0.basis_r, zeros(ComplexF64, size(rho0.data)))
    # (du,u,h,p,t)
    Ham_S, Ham_int, H_lambda, H_S, H_Sinv, Ham_B, W0, W0_bath, agg, FCProd, aggIndices, vibindices, elementtype = p
    tmp1 = copy(W0.data)
    tmp2 = copy(W0.data)
    tmp3 = copy(W0.data)
    dmaster_(t, rho::T, drho::T, history_fun, p) =
        dmaster_ansatz(t, rho, drho, history_fun, tmp1, tmp2, tmp3, p, history_fun)
    # tspan_ = convert(Vector{float(eltype(tspan))}, tspan)
    tspan_ = deepcopy(tspan)
    x0 = rho0.data
    state = T(rho0.basis_l, rho0.basis_r, rho0.data)
    dstate = T(rho0.basis_l, rho0.basis_r, rho0.data)
    OpenQuantumSystems.integrate_delayed(
        tspan_,
        dmaster_,
        history_fun,
        x0,
        state,
        dstate,
        fout;
        p = p,
        reltol = reltol,
        abstol = abstol,
        alg = alg,
        kwargs...,
    )
end

function dmaster_ansatz(
    t::AbstractFloat,
    rho::T,
    drho::T,
    history_fun,
    tmp1::Array,
    tmp2::Array,
    tmp3::Array,
    p,
    h
) where {B<:Basis,T<:Operator{B,B}}
    Ham_S, Ham_int, H_lambda, H_S, H_Sinv, Ham_B, W0, W0_bath, agg, FCProd, aggIndices, vibindices, elementtype = p
    # Ham_II_t = getInteractionHamIPicture(Ham_S, Ham_int, t)
    Ham_II_t = getInteractionHamIPictureA(Ham_int.data, H_lambda, H_S, H_Sinv, t)
    kernel_integrated, err = QuadGK.quadgk(
        s -> MemoryKernel(t, s, tmp1, tmp2, tmp3, h, p, Ham_II_t),
        0,
        t,
        rtol = 1e-3,
        atol = 1e-3
    )

    LinearAlgebra.mul!(tmp1, Ham_II_t, W0_bath.data, -elementtype(im), zero(elementtype))
    LinearAlgebra.mul!(tmp1, W0_bath.data, Ham_II_t, elementtype(im), one(elementtype))
    K_traced = trace_bath(tmp1, agg, FCProd, aggIndices, vibindices)
    # rho0 = trace_bath(W0, agg, FCProd, aggIndices, vibindices)
    # println(t)
    # println(K_traced)
    # println(K_traced .* rho0.data)

    # println(data)
    LinearAlgebra.mul!(drho.data, one(elementtype), K_traced, -elementtype(im), zero(elementtype))
    data = kernel_integrated .* rho.data
    LinearAlgebra.mul!(drho.data, one(elementtype), data, -elementtype(1), one(elementtype))
    # println(-1im * (K_traced .* rho.data))
    # println(kernel_integrated .* rho.data)

    return drho
end

function MemoryKernel(t, s, tmp1, tmp2, tmp3, h, p, Ham_II_t)
    Ham_S, Ham_int, H_lambda, H_S, H_Sinv, Ham_B, W0, W0_bath, agg, FCProd, aggIndices, vibindices, elementtype = p
    # sprintln(t)
    # Ham_II_s = getInteractionHamIPictureA(Ham_S, Ham_int, s)
    tmp3 .= getInteractionHamIPictureA(Ham_int.data, H_lambda, H_S, H_Sinv, s)
    # LinearAlgebra.mul!(tmp1, H_S, exp(-1im * H_lambda * s), one(elementtype), zero(elementtype))
    # LinearAlgebra.mul!(tmp3, tmp1, H_Sinv, one(elementtype), zero(elementtype))

    # LinearAlgebra.mul!(tmp1, adjoint(tmp3), Ham_int.data, one(elementtype), zero(elementtype))
    # LinearAlgebra.mul!(tmp2, tmp1, tmp3, one(elementtype), zero(elementtype))
    # println(tmp2)
    rho = h(p, s)

    if (typeof(rho) <: Operator)
        rho = rho.data
    end
    # println("rho", rho)
    U_b_s = evolutionOperator(Ham_B, s)
    W0_bath_int_s = U_b_s' * W0_bath * U_b_s
    tmp2 .= ad(rho, W0_bath_int_s.data, agg, FCProd, aggIndices, vibindices)

    QuantumOpticsBase.mul!(tmp1, tmp3, tmp2, elementtype(1), zero(elementtype))
    QuantumOpticsBase.mul!(tmp1, tmp2, tmp3, -elementtype(1), one(elementtype))

    QuantumOpticsBase.mul!(tmp2, Ham_II_t, tmp1, elementtype(1), zero(elementtype))
    QuantumOpticsBase.mul!(tmp2, tmp1, Ham_II_t, -elementtype(1), one(elementtype))
    MK_traced = trace_bath(tmp2, agg, FCProd, aggIndices, vibindices)
    # println(MK_traced)
    return MK_traced
end


text:
using SumOfSquares
using DynamicPolynomials
using MosekTools
using Test

include("sparse_polynomials.jl")

function test(test_function, factory)
    sparse_model = sos_lower_bound(test_function, factory, VariableSparsity())
    dense_model  = sos_lower_bound(test_function, factory, NoSparsity())
    return abs(objective_value(sparse_model) - objective_value(dense_model))
end

@testset "chordal_sos" begin
    factory = Mosek.Optimizer

    val = test(chained_singular(12), factory)
    @test abs(val) < 1e-6

    val = test(broyden_banded(6), factory)
    @test abs(val) < 1e-6

    val = test(broyden_tridiagonal(12), factory)
    @test abs(val) < 1e-6

    val = test(chained_wood(12), factory)
    @test abs(val) < 1e-6

    val = test(generalized_rosenbrock(12), factory)
    @test abs(val) < 1e-6

end


text:
@with_kw struct Pars
    # Primitives for any model
    σ²ₐ    = 1.0       # variance of fundamental θ
    σ²ᵤ    = 1.0       # variance of noise traders
    ca     = 0.01      # Parameter for cost function
    cb     = 2.0       # Parameter for cost function
    aₜ     = 0.0       # Mean of distiribution for a
    εₜ     = 0.0       # Mean of distribution for epsilon

    # Switches
    type::Symbol       = :CES   # Payoff function
end

@with_kw struct BPars

    ## Economy
    β   # information precision economy
    σ²  # uncertainty economy

    # sliced by signal
    w_z  # weight on z market
    w_zε # weight on ε in the z part
    w_x  # weight on x market
    w_xε # weight on ε in the x=z part

    # sliced by fund/noise
    w_p  # weight on prior in the price
    w_ε  # weight on ε in the price
    w_a  # weight on a in the price

    ## Market parameters
    βᵢ   # information precision market
    σ²ᵢ  # uncertainty market

    # sliced by signal
    wᵢ_z  # weight on z market
    wᵢ_zε # weight on ε in the z part
    wᵢ_x  # weight on x market
    wᵢ_xε # weight on ε in the x=z part

    # sliced by fund/noise
    wᵢ_p  # weight on prior in the price
    wᵢ_ε  # weight on ε in the price
    wᵢ_a  # weight on a in the price

    ## Objective parameters
    σ²ᵢ_obj  # uncertainty objective observer

    # sliced by signal
    wᵢ_obj_z  # weight on z objective observer
    wᵢ_obj_zε # weight on ε in the z part
    #wᵢ_obj_x  # weight on x objective observer
    #wᵢ_obj_xε # weight on ε in the x=z part

    # sliced by fund/noise
    wᵢ_obj_p  # weight on prior
    wᵢ_obj_ε  # total weight on ε
    wᵢ_obj_a  # total weight on a

    ## Individual parameters
    βᵢⱼ    # information precision i
    σ²ᵢⱼ   # uncertainty trader i

    # sliced by signal
    wᵢⱼ_z  # weight on z trader i
    wᵢⱼ_zε  # weight on ε in z trader i
    wᵢⱼ_x  # weight on x trader i
    wᵢⱼ_xε # weight on ε in the x=z part

    # sliced by fund/noise
    wᵢⱼ_p  # weight on prior trader i
    wᵢⱼ_a  # weight on private and public signal combined trader i
    wᵢⱼ_ε # weight on ε in the x=z part

    # Additional parameters for convenience
    Ωᵢⱼ_x # scaling factor for private signal
    Ωᵢⱼ_z # scaling factor for public signal
    Ωᵢⱼ_p # scaling factor for prior
    Ωᵢⱼ_ε # scaling factor for noise

    # From the expected utility with linear utility
    σ²ₐi_tilde
    σ²_i_ε_tilde
end


text:
# Copyright 2021 The FirstOrderLp Authors
#
# 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.

@enum OptimalityNorm L_INF L2

# Relevant readings on infeasibility certificates:
# (1) https://docs.mosek.com/modeling-cookbook/qcqo.html provides references
# explaining why the primal rays imply dual infeasibility and dual rays imply
# primal infeasibility.
# (2) The termination criteria for Mosek's linear programming optimizer
# https://docs.mosek.com/9.0/pythonfusion/solving-linear.html.
# (3) The termination criteria for OSQP is in section 3.3 of
# https://web.stanford.edu/~boyd/papers/pdf/osqp.pdf.
# (4) The termination criteria for SCS is in section 3.5 of
# https://arxiv.org/pdf/1312.3039.pdf.

"A description of solver termination criteria."
mutable struct TerminationCriteria
  "The norm that we are measuring the optimality criteria in."
  optimality_norm::OptimalityNorm

  # Let p correspond to the norm we are using as specified by optimality_norm.
  # If the algorithm terminates with termination_reason =
  # TERMINATION_REASON_OPTIMAL then the following hold:
  # | primal_objective - dual_objective | <= eps_optimal_absolute +
  #  eps_optimal_relative * ( | primal_objective | + | dual_objective | )
  # norm(primal_residual, p) <= eps_optimal_absolute + eps_optimal_relative *
  #  norm(right_hand_side, p)
  # norm(dual_residual, p) <= eps_optimal_absolute + eps_optimal_relative *
  #   norm(objective_vector, p)
  # It is possible to prove that a solution satisfying the above conditions
  # also satisfies SCS's optimality conditions (see link above) with ϵ_pri =
  # ϵ_dual = ϵ_gap = eps_optimal_absolute = eps_optimal_relative. (ϵ_pri,
  # ϵ_dual, and ϵ_gap are SCS's parameters).

  """
  Absolute tolerance on the duality gap, primal feasibility, and dual
  feasibility.
  """
  eps_optimal_absolute::Float64

  """
  Relative tolerance on the duality gap, primal feasibility, and dual
  feasibility.
  """
  eps_optimal_relative::Float64

  """
  If the following two conditions hold we say that we have obtained an
  approximate dual ray, which is an approximate certificate of primal
  infeasibility.
    (1) dual_ray_objective > 0.0,
    (2) max_dual_ray_infeasibility / dual_ray_objective <=
        eps_primal_infeasible.
  """
  eps_primal_infeasible::Float64

  """
  If the following three conditions hold we say we have obtained an
  approximate primal ray, which is an approximate certificate of dual
  infeasibility.
    (1) primal_ray_linear_objective < 0.0,
    (2) max_primal_ray_infeasibility / (-primal_ray_linear_objective) <=
        eps_dual_infeasible,
    (3) primal_ray_quadratic_norm / (-primal_ray_linear_objective) <=
        eps_dual_infeasible.
  """
  eps_dual_infeasible::Float64

  """
  If termination_reason = TERMINATION_REASON_TIME_LIMIT then the solver has
  taken at least time_sec_limit time.
  """
  time_sec_limit::Float64

  """
  If termination_reason = TERMINATION_REASON_ITERATION_LIMIT then the solver has
  taken at least iterations_limit iterations.
  """
  iteration_limit::Int32

  """
  If termination_reason = TERMINATION_REASON_KKT_MATRIX_PASS_LIMIT then
  cumulative_kkt_matrix_passes is at least kkt_pass_limit.
  """
  kkt_matrix_pass_limit::Float64
end

function construct_termination_criteria(;
  optimality_norm = L2,
  eps_optimal_absolute = 1.0e-6,
  eps_optimal_relative = 1.0e-6,
  eps_primal_infeasible = 1.0e-8,
  eps_dual_infeasible = 1.0e-8,
  time_sec_limit = Inf,
  iteration_limit = typemax(Int32),
  kkt_matrix_pass_limit = Inf,
)
  return TerminationCriteria(
    optimality_norm,
    eps_optimal_absolute,
    eps_optimal_relative,
    eps_primal_infeasible,
    eps_dual_infeasible,
    time_sec_limit,
    iteration_limit,
    kkt_matrix_pass_limit,
  )
end

function validate_termination_criteria(criteria::TerminationCriteria)
  if criteria.eps_primal_infeasible < 0
    error("eps_primal_infeasible must be nonnegative")
  end
  if criteria.eps_dual_infeasible < 0
    error("eps_dual_infeasible must be nonnegative")
  end
  if criteria.time_sec_limit <= 0
    error("time_sec_limit must be positive")
  end
  if criteria.iteration_limit <= 0
    error("iteration_limit must be positive")
  end
  if criteria.kkt_matrix_pass_limit <= 0
    error("kkt_matrix_pass_limit must be positive")
  end
end

"""
Information about the quadratic program that is used in the termination
criteria. We store it in this struct so we don't have to recompute it.
"""
struct CachedQuadraticProgramInfo
  l_inf_norm_primal_linear_objective::Float64
  l_inf_norm_primal_right_hand_side::Float64
  l2_norm_primal_linear_objective::Float64
  l2_norm_primal_right_hand_side::Float64
end

function cached_quadratic_program_info(qp::QuadraticProgrammingProblem)
  return CachedQuadraticProgramInfo(
    norm(qp.objective_vector, Inf),
    norm(qp.right_hand_side, Inf),
    norm(qp.objective_vector, 2),
    norm(qp.right_hand_side, 2),
  )
end

"""
Check if the algorithm should terminate declaring the optimal solution is found.
"""
function optimality_criteria_met(
  optimality_norm::OptimalityNorm,
  abs_tol::Float64,
  rel_tol::Float64,
  convergence_information::ConvergenceInformation,
  qp_cache::CachedQuadraticProgramInfo,
)

  ci = convergence_information
  abs_obj = abs(ci.primal_objective) + abs(ci.dual_objective)
  gap = abs(ci.primal_objective - ci.dual_objective)

  if optimality_norm == L_INF
    primal_err = ci.l_inf_primal_residual
    primal_err_baseline = qp_cache.l_inf_norm_primal_right_hand_side
    dual_err = ci.l_inf_dual_residual
    dual_err_baseline = qp_cache.l_inf_norm_primal_linear_objective
  elseif optimality_norm == L2
    primal_err = ci.l2_primal_residual
    primal_err_baseline = qp_cache.l2_norm_primal_right_hand_side
    dual_err = ci.l2_dual_residual
    dual_err_baseline = qp_cache.l2_norm_primal_linear_objective
  else
    error("Unknown optimality_norm")
  end

  return dual_err < abs_tol + rel_tol * dual_err_baseline &&
         primal_err < abs_tol + rel_tol * primal_err_baseline &&
         gap < abs_tol + rel_tol * abs_obj

end

# TODO: Allow the termination norm to change for infeasibility
# detection.
"""
Check if the algorithm should terminate declaring the primal is infeasible.
"""
function primal_infeasibility_criteria_met(
  eps_primal_infeasible::Float64,
  infeasibility_information::InfeasibilityInformation,
)
  ii = infeasibility_information
  if ii.dual_ray_objective <= 0.0
    return false
  end
  return ii.max_dual_ray_infeasibility / ii.dual_ray_objective <=
         eps_primal_infeasible
end

"""
Check if the algorithm should terminate declaring the dual is infeasible.
"""
function dual_infeasibility_criteria_met(
  eps_dual_infeasible::Float64,
  infeasibility_information::InfeasibilityInformation,
)
  ii = infeasibility_information
  if ii.primal_ray_linear_objective >= 0.0
    return false
  end
  return ii.max_primal_ray_infeasibility / (-ii.primal_ray_linear_objective) <=
         eps_dual_infeasible &&
         ii.primal_ray_quadratic_norm / (-ii.primal_ray_linear_objective) <=
         eps_dual_infeasible
end

"""
Checks if the given iteration_stats satisfy the termination criteria. Returns
a TerminationReason if so, and false otherwise.
"""
function check_termination_criteria(
  criteria::TerminationCriteria,
  qp_cache::CachedQuadraticProgramInfo,
  iteration_stats::IterationStats,
)
  for convergence_information in iteration_stats.convergence_information
    if optimality_criteria_met(
      criteria.optimality_norm,
      criteria.eps_optimal_absolute,
      criteria.eps_optimal_relative,
      convergence_information,
      qp_cache,
    )
      return TERMINATION_REASON_OPTIMAL
    end
  end
  for infeasibility_information in iteration_stats.infeasibility_information
    if primal_infeasibility_criteria_met(
      criteria.eps_primal_infeasible,
      infeasibility_information,
    )
      return TERMINATION_REASON_PRIMAL_INFEASIBLE
    end
    if dual_infeasibility_criteria_met(
      criteria.eps_dual_infeasible,
      infeasibility_information,
    )
      return TERMINATION_REASON_DUAL_INFEASIBLE
    end
  end
  if iteration_stats.iteration_number >= criteria.iteration_limit
    return TERMINATION_REASON_ITERATION_LIMIT
  elseif iteration_stats.cumulative_kkt_matrix_passes >=
         criteria.kkt_matrix_pass_limit
    return TERMINATION_REASON_KKT_MATRIX_PASS_LIMIT
  elseif iteration_stats.cumulative_time_sec >= criteria.time_sec_limit
    return TERMINATION_REASON_TIME_LIMIT
  else
    return false # Don't terminate.
  end
end

function termination_reason_to_string(termination_reason::TerminationReason)
  # Strip TERMINATION_REASON_ prefix.
  return string(termination_reason)[20:end]
end


text:
module Query

import IterableTables
using DataValues
using MacroTools: postwalk
using QueryOperators

export @from, @query, @count, Grouping, key

export @map, @filter, @groupby, @orderby, @orderby_descending, @unique,
	@thenby, @thenby_descending, @groupjoin, @join, @mapmany, @take, @drop

export @select, @rename, @mutate, @dissallowna, @dropna, @replacena

export isna, NA

include("query_utils.jl")
include("query_translation.jl")
include("standalone_query_macros.jl")
include("table_query_macros.jl")

macro from(range::Expr, body::Expr)
	if range.head!=:call || (range.args[1]!=:in && range.args[1]!=in)
		error()
	end

	if body.head!=:block
		error()
	end

	body.args = filter(i->!isa(i, LineNumberNode),body.args)

	insert!(body.args,1,:( @from $(range.args[2]) in $(range.args[3]) ))

	translate_query(body)

	return body.args[1]
end

macro query(range::Symbol, body::Expr)
	if body.head!=:block
		error()
	end

	f_arg = gensym()
	x = x = esc(:($f_arg -> $Query.@from $in($range, $f_arg) $body))
	return x
end

end # module


text:
#   Unit examples of the FourierAnalysis Package for julia language
#
#   MIT License
#   Copyright (c) 2019-2022,
#   Marco Congedo, CNRS, Grenobe, France:
#   https://sites.google.com/site/marcocongedo/home

# ? CONTENTS :
#   This example shows hot to create and plot tapering window for FFT.

#   ~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~  #
#                                                                             #
#   ~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~  #

############ plot tapering windows #############
using FourierAnalysis, Plots

# using the standard plot function
tapers=[TaperKind(i) for i=1:8]
X=zeros(t, 8)
for i=1:8 X[:, i] = taper(tapers[i], t).y end
mylabels=Array{String}(undef, 1, 8)
for i=1:8 mylabels[1, i]=string(tapers[i]) end
plot(X; labels=mylabels)

# discrete prolate spheroid sequences
H=taper(slepian, t, α=4, n=7)
# This below is a plot recipe declared in the `recipes.jl unit`
plot(H)

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


text:
using YStockData
using Base.Test
using Base.Dates
using DataFrames

# write your own tests here
@test typeof(YStockData.get_previous_close("IBM").value) == Float64
@test typeof(YStockData.get_company_name("IBM").value) <: AbstractString
typeof(YStockData.get("IBM")["previous_close"].value) == Float64

@test size(YStockData.get_historical("IBM", now()-Year(1), now()), 1) > 200
@test size(YStockData.get_historical(["IBM", "YHOO"], now()-Year(1), now()), 1) > 400