text:
" Outputs the unstructured mesh with the solution fields "
type VTKOutput <: Output
end

"""
    Create the lists of Cells and Points for Outputting VTK

    Returns a tuple of `cells, points`
"""
function createVTKGeometry(mesh::Mesh)
    # Create the cells
    elems = mesh.elements
    n_elem = length(elems)

    cells = Array{MeshCell}(n_elem)

    for i in 1:n_elem
        elem = mesh.elements[i]
        cells[i] = MeshCell(VTKCellTypes.VTK_QUAD, [ node.id for node in elem.nodes ])
    end

    # Create the Points array
    nodes = mesh.nodes
    n_nodes = length(nodes)

    points = Matrix{Float64}((2, n_nodes))

    for i in 1:n_nodes
        node = mesh.nodes[i]
        points[1,i] = node.coords[1]
        points[2,i] = node.coords[2]
    end

    return cells, points
end

"""
    Creates the point data suitable for outputing to VTK

    returns a Dict{VariableName, Array{Float64}}
"""
function createPointData(solver::Solver, solution::Array)
    # One entry per variable
    point_data = Dict{String, Array{Float64}}()

    sys = solver.system
    mesh = sys.mesh
    nodes = mesh.nodes
    n_nodes = length(nodes)

    for var in sys.variables
        id = var.id
        data = Array{Float64}(n_nodes)

        for i in 1:n_nodes
            data[i] = solution[mesh.nodes[i].dofs[id]]
        end

        point_data[var.name] = data
    end

    return point_data
end

" Write out a VTK file "
function output(out::VTKOutput, solver::Solver, filebase::String)
    startLog(main_perf_log, "output(VTKOutput)")

    serialized_solution = serializeToZero(solver.solution)

    if MPI.Comm_rank(MPI.COMM_WORLD) == 0
        cells, points = createVTKGeometry(solver.system.mesh)

        vtkfile = vtk_grid(filebase, points, cells)

        point_data = createPointData(solver, serialized_solution)

        for name_data in point_data
            vtk_point_data(vtkfile, name_data[2], name_data[1])
        end

        outfiles = vtk_save(vtkfile)
    end

    stopLog(main_perf_log, "output(VTKOutput)")
end


text:
using StructuralIdentifiability
using ModelingToolkit

@parameters b c g
@variables t S(t) I(t) R(t) N(t)
N = S + I + R
D = Differential(t);

sir_eqs1 = [
  D(S) ~ -b*S*I,
  D(I) ~ b*S*I-g*I,
  D(R) ~ g*I
]
sir_ode1 = ODESystem(sir_eqs1, t, name=:SIR1)

sir_eqs2 = [
  D(S) ~ -b*c*S*I,
  D(I) ~ b*c*S*I-g*I,
  D(R) ~ g*I
]
sir_ode2 = ODESystem(sir_eqs2, t, name=:SIR2)

sir_eqs3 = [
  D(S) ~ -b*c*S*I/N,
  D(I) ~ b*c*S*I/N-g*I,
  D(R) ~ g*I
]
sir_ode3 = ODESystem(sir_eqs3, t, name=:SIR3)

@variables y(t)
measured_quantities1 = [y ~ b*S*I];

funcs_to_check1 = [b,g,b/g]
ident1 = assess_identifiability(sir_ode1; measured_quantities= measured_quantities1, funcs_to_check = funcs_to_check1)

funcs_to_check2 = [b,c,g,b*c]
measured_quantities2 = [y ~ b*c*S*I]
ident2 = assess_identifiability(sir_ode2; measured_quantities= measured_quantities2, funcs_to_check = funcs_to_check2)

funcs_to_check3 = [b,c,g,b*c]
measured_quantities3 = [y ~ b*c*S*I/N]
ident3 = assess_identifiability(sir_ode3; measured_quantities= measured_quantities3, funcs_to_check = funcs_to_check3)

@variables n(t)
measured_quantities3a = [y ~ b*c*S*I/N, n ~ N]
ident3a = assess_identifiability(sir_ode3; measured_quantities= measured_quantities3a, funcs_to_check = funcs_to_check3)

exp_eqs = [
  D(I) ~ b*S*I-g*I,
  D(R) ~ g*I
]
exp_ode = ODESystem(exp_eqs, t, name=:EXP)
assess_identifiability(exp_ode;
                       measured_quantities= [y ~ b*S*I],
                       funcs_to_check = [b,g])

ext: jl


text:
using Documenter, Metrics

makedocs(;
    modules=[Metrics],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
        "Regression Metrics" => "regression.md",
        "Classification Metrics" => "classification.md",
        "NLP Metrics" => "nlp.md",
        "CV Metrics" => "cv.md",
        "Ranking Metrics" => "rank.md",
        "Utils" => "utils.md"
    ],
    repo="https://github.com/AdarshKumar712/Metrics.jl/blob/{commit}{path}#L{line}",
    sitename="Metrics.jl",
    authors="Adarsh Kumar <Adarshkumar712.ak@gmail.com>",
    assets=String[],
)

deploydocs(;
    repo="github.com/AdarshKumar712/Metrics.jl",
    target = "build",
    push_preview = true
)


text:
using .MPI
import Libdl

const depsfile = joinpath(@__DIR__, "..", "deps", "deps.jl")
if isfile(depsfile)
    include(depsfile)
else
    error("HDF5 is not properly installed. Please run Pkg.build(\"HDF5\") ",
          "and restart Julia.")
end

# Check whether the HDF5 libraries were compiled with parallel support.
HAS_PARALLEL[] = Libdl.dlopen(libhdf5) do lib
    Libdl.dlsym(lib, :H5Pget_fapl_mpio, throw_error=false) !== nothing
end

# Low-level MPI handles.
const MPIHandle = Union{MPI.MPI_Comm, MPI.MPI_Info}

# MPI.jl wrapper types.
const MPIHandleWrapper = Union{MPI.Comm, MPI.Info}

const H5MPIHandle = let csize = sizeof(MPI.MPI_Comm)
    @assert csize in (4, 8)
    csize == 4 ? API.Hmpih32 : API.Hmpih64
end

h5_to_mpi(handle::H5MPIHandle) = reinterpret(MPI.MPI_Comm, handle)

mpi_to_h5(handle::MPIHandle) = reinterpret(H5MPIHandle, handle)
mpi_to_h5(mpiobj::MPIHandleWrapper) = mpi_to_h5(mpiobj.val)

# Set MPIO properties in HDF5.
# Note: HDF5 creates a COPY of the comm and info objects.
function set_fapl_mpio(fapl_id, comm::MPI.Comm, info::MPI.Info)
    h5comm = mpi_to_h5(comm)
    h5info = mpi_to_h5(info)
    API.h5p_set_fapl_mpio(fapl_id, h5comm, h5info)
end

# Retrieves the copies of the comm and info MPIO objects from the HDF5 property list.
function get_fapl_mpio(fapl_id)
    h5comm, h5info = API.h5p_get_fapl_mpio(fapl_id, H5MPIHandle)
    comm = MPI.Comm(h5_to_mpi(h5comm))
    info = MPI.Info(h5_to_mpi(h5info))
    comm, info
end

function check_hdf5_parallel()
    has_parallel() && return
    error(
        "HDF5.jl has no parallel support." *
        " Make sure that you're using MPI-enabled HDF5 libraries, and that" *
        " MPI was loaded before HDF5." *
        " See HDF5.jl docs for details."
    )
end

"""
    h5open(filename, [mode="r"], comm::MPI.Comm, [info::MPI.Info]; pv...)

Open or create a parallel HDF5 file using the MPI-IO driver.

Equivalent to `h5open(filename, mode; fapl_mpio=(comm, info), pv...)`.
Throws an informative error if the loaded HDF5 libraries do not include parallel
support.

See the [HDF5 docs](https://portal.hdfgroup.org/display/HDF5/H5P_SET_FAPL_MPIO)
for details on the `comm` and `info` arguments.
"""
function h5open(
        filename::AbstractString, mode::AbstractString,
        comm::MPI.Comm, info::MPI.Info = MPI.Info(); pv...
    )
    check_hdf5_parallel()
    h5open(filename, mode; fapl_mpio=(comm, info), pv...)
end

h5open(filename::AbstractString, comm::MPI.Comm, args...; pv...) =
    h5open(filename, "r", comm, args...; pv...)


text:
module LogicProcessing

#################################################################################
# Creador: Juan David Zorrilla Herrergit config --global user.email "EMAIL"a y Cristhian Silvita Cartagena            #
# Fecha: 17/01/2018                                                             #
# Empresa: KNAR S.A.                                                            #
# Lenguaje: Julia 0.6.2                                                         #
# Titulo: Librerias.jl                                                          #
# Versión Actual: 2.0                                                           #
# Descripción: Este programa es la creacion de unas librerias necesarias para la#
# generación de arreglos de numeros aleatorios utilizando la CPU en donde ya te-#
# nemos la generacion de matrices aleatorias y una funcion logica (||) implemen-#
# tada                                                                          #
# Estado: TRABAJANDO EN LA SIGUIENTE git config --global user.email "EMAIL"VERSIÓN                                    #
#################################################################################
## Variables Globales
global Mat=Array[], num=0, ID=Array[],Res=Array[]
################################################################################
export Evento,LogicOpMat,DispMat
using Distributions
using PyPlot

# FUNCIÓN QUE ME RETORNA UN NUMERO ALEATORIO VARIANDO SEGÚN LA FUNCIÓN DE DISTRIBUCIÓN
function probability_function(type_fc;μ=NaN,σ=NaN,α=NaN,θ=NaN)
  # Función con distribución de tipo Normal
  if type_fc=="normal"        # Tipo de función
    if isequal(μ,NaN)||isequal(σ,NaN) # Verifica que se hayan ingresado los parametros correctos
      println("Se requiere μ(Media) y σ(Desviación Estandar)")
      return NaN
    end
    return rand(Normal(μ,σ))  #Retorna el valor aleatorio
  end
  # Función con distribución de tipo Weibull
  if type_fc=="weibull" # Tipo de función
    if isequal(α,NaN)||isequal(θ,NaN) # Verifica que se hayan ingresado los parametros correctos
      println("Se requiere α y θ")
      return NaN
    end
    return rand(Weibull(α,θ)) #Retorna el valor aleatorio
  end
  # Función con distribución de tipo Exponencial
  if type_fc=="exponential" # Tipo de función
    if isequal(θ,NaN)       # Verifica que se hayan ingresado los parametros correctos
      println("Se requiere θ")
      return NaN
    end
    return rand(Exponential(θ)) #Retorna el valor aleatorio
  end
  # Error al ingresar los parametros
  println("Ingreso mal la función de probabilidad");
  return NaN
end
#FUNCION QUE ME RETORNA UNA MATRIZ DE VALORES QUE ME INDICAN LAS FALLAS DE UN EQUIPO
function VecAlea(t_final,fix,type_fc;miui=NaN,sigmai=NaN,alphai=NaN,thetai=NaN)
  timenow=0;  # Tiempo actual
  fail=[0.0]  #Vector FLOAT que indica las FALLAS
  while true #Continue si el tiempo actual es menor al tiempo final de la simulacion
    Numb=probability_function(type_fc;μ=miui,σ=sigmai,α=alphai,θ=thetai); #Hallamos el numero con las caracteristicas inficadas
    aux=+(timenow,Numb,fix) #Hallamos el valor del tiempo que estariamos si lo repararamos
    if aux>t_final
      break;  #Si se salio del tiempo establecido, salir del ciclo
    end
    push!(fail,timenow+Numb)  #Agregamos el momento en que se creo la falla
    timenow=aux;
    push!(fail,timenow)git config --global user.email "EMAIL"
  end
  shift!(fail)  #Eliminamos el 0.0 del vector de fallas ya que solo era para declaracion
  return fail#Retornamos el vector de fallas
end

# FUNCION QUE NOS RETORNA UNA MATRIZ DE EVENTOS TENIENDO EN CUENTA LAS CARACTERISTICAS DE CADA EVENTO
# Recibe como parametros el Numero de iteraciones, el tiempo de simulacion, el tipo de funcion de
# distribucion con sus respectivos parametros despues de ;
function Evento(NumIter,t_final,fix,type_fc,equipo::String;miu=NaN,sigma=NaN,alpha=NaN,theta=NaN)
  global Mat, ID, num
  num+=1
  push!(ID,[num,equipo])
  Aux=[VecAlea(t_final,fix,type_fc;miui=miu,sigmai=sigma,alphai=alpha,thetai=theta)] # Declaramos el vector aleatorio con el primer caso
  for i = 2:NumIter # Realizamos el mismo procedimiento anterior pero Numero de iteraciones - 1
    push!(Aux,VecAlea(t_final,fix,type_fc;miui=miu,sigmai=sigma,alphai=alpha,thetai=theta)) # Agregamos el resultado al vector aleatorio
  end
  push!(Mat,Aux)
  return "Exito: Evento Equipo $equipo Id $num" # Retornamos un vector de vectores de diferentes tamaños
end

# FUNCION QUE NOS PERMITE OPERAR UN PAR DE VECTORES CON UNA FUNCION LOGICA EN CONTINUO
# Recibe como parametros el primer vector "VecA", el segundo "VecB" y la operacion logica "AND" "OR"
function LogicOpVec(VecA,VecB,LogOp)
  eventF=[0.0]  # Declaramos el primer valor de la salida para que posea el formato deseado en la salida
  if LogOp=="OR"  # Seccion de operacion "OR"
    iai=1 # Indice de inicio del Vector A
    iaf=2 # Indice de final del Vector A
    ibi=1 # Indice de inicio del Vector B
    ibf=2 # Indice de final del Vector B
    flagA,flagB=1,1 # Bandera de los vectores A y B que indica el final de los vectores
    findend=1; # Bandera que permite que se agregue el comienzo del vector A o B en la salida
    while true # Ciclo que finalizaa cuando ya se hayan analizado ambos vectoes

      if iaf>endof(VecA[1]) # Si el indice final de A sobrepasa del Vector A regrese
        iaf-=2
        iai-=2
        flagA=0
      end
      if ibf>endof(VecB[1]) # Si el indice final de B sobrepasa del Vector B regrese
        ibf-=2
        ibi-=2
        flagB=0;
      end
      # Si las banderas son 0 y algun indice es mayor al tamaño del vector salga del while
      if (flagA==0 && ibf>endof(VecB[1])) || (flagB==0 && iaf>endof(VecA[1])) || (flagA==0 && flagB==0)
        break;
      end
      # Se agrega al vector de salida el menor de los inicios de los vectores
      if VecA[1][iai]<VecB[1][ibi] && flagA==1
        if findend==1 # Siempre y cuagit config --global user.email "EMAIL"ndo el anterior dato sea un final
          push!(eventF,VecA[1][iai])
        end
        # Se analiza el caso en el que el inicio del vector B este antes del final del vector A
        if VecA[1][iaf]>=VecB[1][ibi] && VecA[1][iaf]<=VecB[1][ibf]
          findend=0;  #  Se activa la bandera indicando que no se ha hallado el final
          iaf+=2  # Pasamos al siguiente caso
          iai+=2
        else
          findend=1 # Significa que podemos agregar el final del evento
          push!(eventF,VecA[1][iaf])  # Agregamos a la salida

          if VecB[1][ibf]<VecA[1][iaf]  # Analizamos el caso que el evento de la vector B este dentro del evento del vector A
            ibf+=2  # Se es asi pasamos al siguente caso de B
            ibi+=2
          end
          # Pasamos al siguente caso de A
          iaf+=2
          iai+=2
        end
        # Miramos el mismo caso del anterior teniendo en cuenta que el menor sera el vector B
      elseif VecA[1][iai]>=VecB[1][ibi] && flagB==1

        if findend==1
          push!(eventF,VecB[1][ibi])
        end

        if VecB[1][ibf]>=VecA[1][iai] && VecB[1][ibf]<=VecA[1][iaf]
          findend=0;
          ibf+=2
          ibi+=2
        else
          findend=1
          push!(eventF,VecB[1][ibf])

          if VecA[1][iaf]<VecB[1][ibf]
            iaf+=2
            iai+=2
          end
          ibf+=2
          ibi+=2
        end
        # Analizamos el caso cuando ya se llego al final del vector A
      elseif VecA[1][iai]<VecB[1][ibi] && flagA==0
        if findend==1 # Si esta disponible agregamos el inicio de B
          push!(eventF,VecB[1][ibi])
        end
        push!(eventF,VecB[1][ibf])  # y el final del vector B
        ibf+=2
        ibi+=2
        # Analizamos el mismo caso anterior pero cuando se llego al final del vector B
      elseif VecA[1][iai]>=VecB[1][ibi] && flagB==0
        if findend==1
          push!(eventF,VecA[1][iai])
        end
        push!(eventF,VecA[1][iaf])
        iaf+=2
        iai+=2
      end
    end
    shift!(eventF)  # Eliminamos el primer valor ya que este fue de declaracion
  end
  return(eventF) # Retornamos el vector de resultado
end

# FUNCION QUE REALIZA LA OPERACION LOGICA DE MATICES
function LogicOpMat(Opera,ID1,ID2) # Recibe La Matriz A, Matriz B y la Operacion
  global Mat,Res,num,ID
  num+=1
  MatA=Mat[ID1]
  MatB=Mat[ID2]

  Numiter=length(MatA)  # Hallamos el numero de iteciones por medio del tamaño de la matriz A
  Res=[LogicOpVec([MatA[1]],[MatB[1]],"OR")] # Declaramos la salida "OUT" haciendo el primer caso
  for i=2:Numiter # Comenzamos desde i=2 porque i=1 se desarrollo en la anterior
    push!(Res,LogicOpVec([MatA[i]],[MatB[i]],"OR")) # Realizamos la operacion por linea
  end

  push!(ID,[num,"Resultado"])
  push!(Mat,Res)
  return "Exito: Operacion Logica $Opera - ID $num"  # Retornamos el valor
end

function DispVect(VecA,t_sim)
  sVecA=length(VecA[1])
  aux=0.0
  for i=1:2:sVecA-2
    aux=+(-(VecA[1][+(i,1)],VecA[1][i]),aux)
  end
  aux=(1-aux/t_sim)*100
  return aux
end

function DispMat(t_sim,Steps)
  global Res
  sRes=length(Res)
  aux=[0.0]
  for i=1:sRes
    push!(aux,DispVect([Res[i]],t_sim))
  end
  shift!(aux)

  plt[:hist](aux,Steps)
  return "Exito: Calculo de disponibilidad"
end

function Tabla()
  global ID,Mat
  return ID,Mat
end

end # module


text:
show(io::IO, c::Colorant)              = _show(io, c)
show(io::IO, c::ColorantUFixed)        = show_ufixed(io, c)
showcompact(io::IO, c::Colorant)       = _showcompact(io, c)
showcompact(io::IO, c::ColorantUFixed) = show_ufixed(io, c)

for N = 1:4
    component = N >= 3 ? (:comp1, :comp2, :comp3, :alpha) : (:comp1, :alpha)
    printargs = Array(Any, 2, N)
    for i = 1:N
        printargs[1,i] = :(show(io, $(component[i])(c)))
        chr = i < N ? ',' : ')'
        printargs[2,i] = :(print(io, $chr))
    end
    @eval begin
        function _show{T}(io::IO, c::Colorant{T,$N})
            print(io, colorant_string(typeof(c)), "{", T, "}(")
            $(printargs[:]...)
        end
    end
    for i = 1:N
        printargs[1,i] = :(showcompact(io, $(component[i])(c)))
    end
    @eval begin
        function _showcompact{T}(io::IO, c::Colorant{T,$N})
            print(io, colorant_string(typeof(c)), "{", T, "}(")
            $(printargs[:]...)
        end
        # Special handling for UFixed types: don't print the giant type name
        function show_ufixed{T,f}(io::IO, c::Colorant{FixedPointNumbers.UFixed{T,f},$N})
            print(io, colorant_string(typeof(c)), "{UFixed", f, "}(")
            $(printargs[:]...)
        end
        function show_ufixed(io::IO, c::Colorant{U8,$N})
            print(io, colorant_string(typeof(c)), "{U8}(")
            $(printargs[:]...)
        end
    end
end


text:
@testset "Decreasing2LP: $(fct_type), dimension $(dim), $(T)" for fct_type in ["vector of variables", "vector affine function"], dim in [2, 3], T in [Int, Float64]
    mock = MOIU.MockOptimizer(MILPModel{T}())
    model = COIB.Decreasing2LP{T}(mock)

    if T == Int
        @test MOI.supports_constraint(model, MOI.VariableIndex, MOI.Integer)
    end
    @test MOI.supports_constraint(
        model,
        MOI.ScalarAffineFunction{T},
        MOI.EqualTo{T},
    )
    @test MOIB.supports_bridging_constraint(
        model,
        MOI.VectorAffineFunction{T},
        CP.Decreasing,
    )

    if T == Int
        x, _ = MOI.add_constrained_variables(model, [MOI.Integer() for _ in 1:dim])
    elseif T == Float64
        x = MOI.add_variables(model, dim)
    end

    fct = if fct_type == "vector of variables"
        MOI.VectorOfVariables(x)
    elseif fct_type == "vector affine function"
        MOIU.vectorize(x)
    else
        @assert false
    end
    c = MOI.add_constraint(model, fct, CP.Decreasing(dim))

    for i in 1:dim
        @test MOI.is_valid(model, x[i])
    end
    @test MOI.is_valid(model, c)

    bridge = MOIBC.bridges(model)[MOI.ConstraintIndex{MOI.VectorOfVariables, CP.Decreasing}(-1)]

    @testset "Bridge properties" begin
        @test MOIBC.concrete_bridge_type(typeof(bridge), MOI.VectorOfVariables, CP.Decreasing) == typeof(bridge)
        @test MOIB.added_constrained_variable_types(typeof(bridge)) == Tuple{Type}[]
        @test MOIB.added_constraint_types(typeof(bridge)) == [(MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T})]

        @test MOI.get(bridge, MOI.NumberOfVariables()) == 0
        @test MOI.get(bridge, MOI.NumberOfConstraints{MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T}}()) == dim - 1

        @test Set(MOI.get(bridge, MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T}}())) == Set(collect(values(bridge.cons)))
    end

    @testset "Set of constraints" begin
        @test length(bridge.cons) == dim - 1
        for i in 1:(dim - 1)
            @test MOI.is_valid(model, bridge.cons[i])
            f = MOI.get(model, MOI.ConstraintFunction(), bridge.cons[i])
            @test length(f.terms) == 2
            @test MOI.get(model, MOI.ConstraintSet(), bridge.cons[i]) == MOI.GreaterThan(zero(T))

            t1 = f.terms[1]
            @test t1.coefficient === one(T)
            @test t1.variable == x[i]

            t2 = f.terms[2]
            @test t2.coefficient === -one(T)
            @test t2.variable == x[i + 1]
        end
    end
end


text:
@testset "poisson" begin
    n = 20
    h = 1/n
    node = zeros((n+1)^2, 2)
    node_mid = zeros(2n^2,2)
    elem = zeros(Int32, 2n^2, 3)
    k = 1
    for i = 1:n
        for j = 1:n
            node[i+(j-1)*(n+1), :] = [(i-1)*h;(j-1)*h]
            node[i+1+(j-1)*(n+1), :] = [i*h;(j-1)*h]
            node[i+(j)*(n+1), :] = [(i-1)*h;(j)*h]
            node[i+1+(j)*(n+1), :] = [(i)*h;(j)*h]
            elem[k,:] = [i+(j-1)*(n+1) i+1+(j-1)*(n+1) i+(j)*(n+1)];
            elem[k+1,:] = [i+1+(j)*(n+1) i+1+(j-1)*(n+1) i+(j)*(n+1)];
            node_mid[k, :] = mean(node[elem[k,:], :], dims=1)
            node_mid[k+1, :] = mean(node[elem[k+1,:], :], dims=1)
            k += 2
        end
    end
    bdnode = Int32[]
    for i = 1:n+1
        for j = 1:n+1
            if i==1 || j == 1 || i==n+1 || j==n+1
                push!(bdnode, i+(j-1)*(n+1))
            end
        end
    end
    x, y = node_mid[:,1], node_mid[:,2]
    f = constant(@. 2y*(1-y)+2x*(1-x))
    a = constant(ones(2n^2))
    u = pde_poisson(node,elem,bdnode,f,a)
    uval = run(sess, u)
    x, y = node[:,1], node[:,2]
    uexact = @. x*(1-x)*y*(1-y)
    @test maximum(abs.(uval - uexact))<1e-3
end


text:
using DataInterpolations, Test

@testset "Interface" begin include("interface.jl") end
@testset "Interpolation Tests" begin include("interpolation_tests.jl") end
@testset "Derivative Tests" begin include("derivative_tests.jl") end


text:
# Dantzig-Wolfe Optim
# Functions for Master Problem
# Edward J. Xu
# 2019.5.28
########################################################################################################################

function setModelMas(mat_e, vec_b, num_sub, gurobi_env)
    ## X1: initial matrix of extreme points
    modMas = Model(solver = GurobiSolver(OutputFlag = 0, gurobi_env))
    (m_e, n_e) = size(mat_e)
    vec_p = vec_b[1: m_e, 1]
    K = 1
    # In this case we do not use a starting set of extreme points.
    # we just use a dummy starting column
    @variable(modMas, vecLambda[1: K] >= 0 )
    # Remember to consider if we need to maximize or minimize
    @objective(modMas, Max, sum(- 1e6 * vecLambda[j] for j = 1: K) )
    # remember to change "==" if your master-problem uses a different type of constraints!
    @constraint(modMas, vecConsRef[i = 1: m_e], sum(vec_p[i] * vecLambda[j] for j = 1: K) <= vec_p[i])
    @constraint(modMas, vecConsConvex[k = 1: num_sub], sum(vecLambda[j] for j = 1: K) == 1)
    return (modMas, vecConsRef, vecConsConvex, vecLambda)
end

function solveMas(modMas, vecConsRef, vecConsConvex)
    status = solve(modMas)
    if status != :Optimal
        throw("Error: Non-optimal master-problem status")
    end
    vec_pi = getdual(vecConsRef)
    vec_kappa = getdual(vecConsConvex)
    obj_master = getobjectivevalue(modMas)
    ## ensure that vec_pi and vec_kappa are row vectors
    vec_pi = reshape(vec_pi, 1, length(vec_pi))
    vec_kappa = reshape(vec_kappa, 1, length(vec_kappa))
    return (vec_pi, vec_kappa, obj_master)
end

function addColToMaster(modMas::JuMP.Model, modSub::ModelSub, vec_x_result, vecConsRef, consConvex)
    (m, n) = size(modSub.mat_e)
    A0x = modSub.mat_e * vec_x_result
    vec_consTouched = ConstraintRef[]
    vals = Float64[]
    for i = 1: m
        # only insert non-zero elements (this saves memory and may make the master problem easier to solve)
        if A0x[i] != 0
            push!(vec_consTouched, vecConsRef[i])
            push!(vals, A0x[i])
        end
    end
    # add variable to convexity constraint.
    push!(vec_consTouched, consConvex)
    push!(vals, 1)
    objCoef = sum(modSub.vec_l[j] * vec_x_result[j] for j = 1: n)
    # println("objCoef = $objCoef.")
    @variable(
        modMas,
        lambdaNew >= 0,                     # New variable to be added
        objective = objCoef,                # cost coefficient of new varaible in the objective
        inconstraints = vec_consTouched,    # constraints to be modified
        coefficients = vals                 # the coefficients of the variable in those constraints
        )
    return lambdaNew
end


text:
module LP

using JuMP
using Clp

function example()
    # inspired by Linear Algebra and its applications Ch. 9

    m = Model(Clp.Optimizer)

    # unknowns
    @variable(m, 0 <= x1)
    @variable(m, 0 <= x2)
    @variable(m, 0 <= x3)

    # dot(C, x)
    @objective(m, Max, 2x1 + 3x2 + 4x3)

    # A ∈ R(mxn), with m constraints and n unknowns
    @constraint(m, x1 + x2 + x3 <= 50)
    @constraint(m, x1 + 2x2 + 4x3 <= 80)

    optimize!(m)

    println("Optimal X: x1: $(value(x1)) x2: $(value(x2)) x3: $(value(x3))")
    println("Optimal Objective value: $(objective_value(m))")

end # example

function practice_problem_one()
    # inspired by Linear Algebra and its applications Ch. 9

    m = Model(Clp.Optimizer)

    # unknowns
    @variable(m, 0 <= x1)
    @variable(m, 0 <= x2)

    # dot(C, x)
    @objective(m, Max, 2x1 + x2)

    @constraint(m, -x1 + x2 <= 8)
    @constraint(m, 3x1 + 2x2  <= 24)

    optimize!(m)

    println("Optimal X: x1: $(value(x1)) x2: $(value(x2))")
    println("Optimal Objective value: $(objective_value(m))")

end # practice problem one

function exercise_seven()
    # inspired by Linear Algebra and its applications Ch. 9

    m = Model(Clp.Optimizer)

    # unknowns
    @variable(m, 0 <= x1)
    @variable(m, 0 <= x2)

    # dot(C, x)
    @objective(m, Max, 80x1 + 65x2)

    @constraint(m, 2x1 + x2 <= 32)
    @constraint(m, x1 + x2  <= 18)
    @constraint(m, x1 + 3x2  <= 24)

    optimize!(m)

    println("Optimal X: x1: $(value(x1)) x2: $(value(x2))")
    println("Optimal Objective value: $(objective_value(m))")

end # exercise seven

end # module


text:
function _ilog2(n, i)
    if n <= (one(n) << i)
        i
    else
        _ilog2(n, i+1)
    end
end
function ilog2(n::Integer)
    @assert n > zero(n)
    _ilog2(n, zero(n))
end

"""
    GeoMeanBridge{T, F, G, H}

The `GeometricMeanCone` is `SecondOrderCone` representable; see [1, p. 105].
The reformulation is best described in an example.
Consider the cone of dimension 4
```math
t \\le \\sqrt[3]{x_1 x_2 x_3}
```
This can be rewritten as ``\\exists x_{21} \\ge 0`` such that
```math
\\begin{align*}
  t & \\le x_{21},\\\\
  x_{21}^4 & \\le x_1 x_2 x_3 x_{21}.
\\end{align*}
```
Note that we need to create ``x_{21}`` and not use ``t^4`` directly as ``t`` is allowed to be negative.
Now, this is equivalent to
```math
\\begin{align*}
  t & \\le x_{21}/\\sqrt{4},\\\\
  x_{21}^2 & \\le 2x_{11} x_{12},\\\\
  x_{11}^2 & \\le 2x_1 x_2, & x_{12}^2 & \\le 2x_3(x_{21}/\\sqrt{4}).
\\end{align*}
```

[1] Ben-Tal, Aharon, and Arkadi Nemirovski. *Lectures on modern convex optimization: analysis, algorithms, and engineering applications*. Society for Industrial and Applied Mathematics, 2001.
"""
struct GeoMeanBridge{T, F, G, H} <: AbstractBridge
    # Initially, (t, x) is of dimension d so x is dimension (d-1)
    # TODO the sentence below is a little confusing, the n isn't used anywhere
    # We create n new variables so that there are 2^l = d-1+n variables x_i
    # We then need to create 2^l-1 new variables (1+2+...+2^{l-1})
    # So the total number of variables is d+2^l-1
    d::Int
    xij::Vector{MOI.VariableIndex}
    tubc::CI{F, MOI.LessThan{T}}
    socrc::Vector{CI{G, MOI.RotatedSecondOrderCone}}
    nonneg::Union{Nothing, CI{H, MOI.Nonnegatives}}
end
function bridge_constraint(::Type{GeoMeanBridge{T, F, G, H}}, model,
                           f::MOI.AbstractVectorFunction,
                           s::MOI.GeometricMeanCone) where {T, F, G, H}
    d = s.dimension
    if d <= 1
        # TODO change to a standard error: https://github.com/JuliaOpt/MathOptInterface.jl/issues/967
        error("Dimension of GeometricMeanCone must be greater than 1.")
    end
    n = d - 1
    l = ilog2(n)
    N = 1 << l
    socrc = Vector{CI{G, MOI.RotatedSecondOrderCone}}(undef, N - 1)
    f_scalars = MOIU.eachscalar(f)
    t = f_scalars[1]
    SG = MOIU.scalar_type(G)

    if d == 2
        xij = MOI.VariableIndex[]
        tubc = MOIU.normalize_and_add_constraint(
            model, MOIU.operate!(-, T, t, f_scalars[2]), MOI.LessThan(zero(T)),
            allow_modify_function=true)
        nonneg = MOI.add_constraint(model, f_scalars[2:end], MOI.Nonnegatives(1))
        return GeoMeanBridge{T, F, G, H}(d, xij, tubc, socrc, nonneg)
    end

    xij = MOI.add_variables(model, N - 1)
    xl1 = MOI.SingleVariable(xij[1])
    sN = one(T) / √N
    function _getx(i)::SG
        if i > n
            return sN * xl1
        else
            return f_scalars[1 + i]
        end
    end

    # With sqrt(2)^l*t - xl1, we should scale both the ConstraintPrimal and ConstraintDual
    tubc = MOIU.normalize_and_add_constraint(
        model, MOIU.operate!(+, T, t, -sN * xl1), MOI.LessThan(zero(T)),
        allow_modify_function=true)

    offset = offset_next = 0
    for i in 1:l
        num_lvars = 1 << (i - 1)
        offset_next = offset + num_lvars
        for j in 1:num_lvars
            if i == l
                a = _getx(2j - 1)
                b = _getx(2j)
            else
                a = convert(SG, MOI.SingleVariable(xij[offset_next + 2j - 1]))
                b = convert(SG, MOI.SingleVariable(xij[offset_next + 2j]))
            end
            c = MOI.SingleVariable(xij[offset+j])
            socrc[offset + j] = MOI.add_constraint(
                model, MOIU.operate(vcat, T, a, b, c),
                MOI.RotatedSecondOrderCone(3))
        end
        offset = offset_next
    end
    @assert offset == length(socrc)
    return GeoMeanBridge{T, F, G, H}(d, xij, tubc, socrc, nothing)
end

function MOI.supports_constraint(::Type{GeoMeanBridge{T}},
                                ::Type{<:MOI.AbstractVectorFunction},
                                ::Type{MOI.GeometricMeanCone}) where T
    return true
end
MOIB.added_constrained_variable_types(::Type{<:GeoMeanBridge}) = Tuple{DataType}[]
function MOIB.added_constraint_types(::Type{<:GeoMeanBridge{T, F, G}}) where {T, F, G}
    return [(F, MOI.LessThan{T}), (G, MOI.RotatedSecondOrderCone), (G, MOI.Nonnegatives)]
end
function concrete_bridge_type(::Type{<:GeoMeanBridge{T}},
                              H::Type{<:MOI.AbstractVectorFunction},
                              ::Type{MOI.GeometricMeanCone}) where T
    S = MOIU.scalar_type(H)
    A = MOIU.promote_operation(*, T, T, MOI.SingleVariable)
    F = MOIU.promote_operation(+, T, S, A)
    G = MOIU.promote_operation(vcat, T, A, S, MOI.SingleVariable)
    return GeoMeanBridge{T, F, G, H}
end

# Attributes, Bridge acting as a model
MOI.get(b::GeoMeanBridge, ::MOI.NumberOfVariables) = length(b.xij)
MOI.get(b::GeoMeanBridge, ::MOI.ListOfVariableIndices) = b.xij
function MOI.get(b::GeoMeanBridge{T, F},
                 ::MOI.NumberOfConstraints{F, MOI.LessThan{T}}) where {T, F}
    return 1 # t ≤ x_{l1}/sqrt(N)
end
function MOI.get(b::GeoMeanBridge{T, F, G},
                 ::MOI.NumberOfConstraints{G, MOI.RotatedSecondOrderCone}) where {T, F, G}
    return length(b.socrc)
end
function MOI.get(b::GeoMeanBridge{T, F, G},
                 ::MOI.NumberOfConstraints{G, MOI.Nonnegatives}) where {T, F, G}
    return (b.d > 2 ? 0 : 1)
end
function MOI.get(b::GeoMeanBridge{T, F},
                 ::MOI.ListOfConstraintIndices{F, MOI.LessThan{T}}) where {T, F}
    return [b.tubc]
end
function MOI.get(b::GeoMeanBridge{T, F, G},
                 ::MOI.ListOfConstraintIndices{G, MOI.RotatedSecondOrderCone}) where {T, F, G}
    return b.socrc
end
function MOI.get(b::GeoMeanBridge{T, F, G, H},
                 ::MOI.ListOfConstraintIndices{H, MOI.Nonnegatives}) where {T, F, G, H}
    return (b.d > 2 ? MOI.ConstraintIndex{H, MOI.Nonnegatives}[] : [b.nonneg])
end

# References
function MOI.delete(model::MOI.ModelLike, bridge::GeoMeanBridge)
    MOI.delete(model, bridge.xij)
    MOI.delete(model, bridge.tubc)
    MOI.delete(model, bridge.socrc)
    if bridge.d == 2
        MOI.delete(model, bridge.nonneg)
    end
end

# Attributes, Bridge acting as a constraint
function MOI.get(::MOI.ModelLike, ::MOI.ConstraintSet, bridge::GeoMeanBridge)
    return MOI.GeometricMeanCone(bridge.d)
end
function MOI.get(model::MOI.ModelLike, attr::MOI.ConstraintFunction,
                 bridge::GeoMeanBridge{T, F, G, H}) where {T, F, G, H}
    d = bridge.d
    f_scalars = Vector{MOIU.scalar_type(H)}(undef, bridge.d)
    tub = MOI.get(model, attr, bridge.tubc)
    rhs = MOI.constant(MOI.get(model, MOI.ConstraintSet(), bridge.tubc))
    tub = MOIU.operate(-, T, tub, rhs)
    if d == 2
        x = MOI.get(model, attr, bridge.nonneg)
        f_scalars[2] = MOIU.eachscalar(x)[1]
        f_scalars[1] = MOIU.operate(+, T, tub, f_scalars[2])
    else
        t = MOIU.remove_variable(tub, bridge.xij[1])
        f_scalars[1] = t
        n = d - 1
        l = ilog2(n)
        num_lvars = 1 << (l - 1)
        offset = num_lvars - 1
        for j in 1:num_lvars
            if 2j <= bridge.d
                func = MOI.get(model, attr, bridge.socrc[offset + j])
                func_scalars = MOIU.eachscalar(func)
                f_scalars[2j] = func_scalars[1]
                if 2j + 1 <= bridge.d
                    f_scalars[2j + 1] = func_scalars[2]
                end
            end
        end
    end
    return MOIU.vectorize(f_scalars)
end
function _getconstrattr(model, attr, bridge::GeoMeanBridge{T}) where T
    output = Vector{T}(undef, bridge.d)
    output[1] = MOI.get(model, attr, bridge.tubc)
    if bridge.d == 2
        output[2] = MOI.get(model, attr, bridge.nonneg)[1]
    else
        N = length(bridge.xij) + 1
        # div(N, 2) gets layer before original problem variables are involved
        offset = div(N, 2) - 1 # 1 + 2 + ... + n/4
        for i in 1:(bridge.d - 1)
            j = ((i - 1) >> 1) + 1
            k = i - 2(j - 1)
            output[1+i] = MOI.get(model, attr, bridge.socrc[offset + j])[k]
        end
    end
    output
end
function MOI.get(model::MOI.ModelLike, attr::MOI.ConstraintPrimal, bridge::GeoMeanBridge)
    output = _getconstrattr(model, attr, bridge)
    N = length(bridge.xij) + 1
    # the constraint is t - x_l1/sqrt(2^l) ≤ 0, we need to add the value of x_l1
    if bridge.d == 2
        output[1] += MOI.get(model, attr, bridge.nonneg)[1]
    else
        output[1] += MOI.get(model, MOI.VariablePrimal(), bridge.xij[1]) / sqrt(N)
    end
    return output
end
function MOI.get(model::MOI.ModelLike, attr::MOI.ConstraintDual, bridge::GeoMeanBridge)
    output = _getconstrattr(model, attr, bridge)
    if bridge.d == 2
        # The transformation is:
        # [t]                [1 -1] [t]
        # [x] in GeoMean <=> [0  1] [x] in R- x R+
        # Hence the dual transformation is
        # [ 1 0] [u]                 [u]
        # [-1 1] [v] in GeoMean* <=> [v] in R- x R+
        output[2] -= MOI.get(model, attr, bridge.tubc)
    end
    # Otherwise, the transformation is:
    # [t]                                  [t]
    # [x] in GeoMean <=> exists xij s.t. A [x] + B xij in R- x RSOC^(N-1)
    # Hence the dual transformation is
    # A' y in GeoMean* <=>
    # B' y = 0         <=> y in R- x RSOC^(N-1)
    # Here A has at most one nonzero entry by line and column and this entry is 1.
    # Therefore we can A' = inv(A) hence we can use the same `getconstrattr` function
    # than for getting the `ConstraintPrimal`.
    # We don't need to use B' y = 0.
    return output
end


text:
function test_content(f::Base.Callable, fname)
  try
    open(fname) do io
      f(read(io, String))
    end
  finally
    rm(fname, force = true)
  end
end

@testset "Send" begin
  logfile = tempname()
  server = "smtp://127.0.0.1:1025"
  addr = "<julian@julialang.org>"
  mock = joinpath(dirname(@__FILE__), "mock.py")

  cmd = `python3.7 $mock $logfile`
  smtpsink = run(pipeline(cmd, stderr=stdout), wait = false)
  sleep(.5)  # wait for fake smtp server ready

  try
    let  # send with body::IOBuffer
      body = IOBuffer("body::IOBuffer test")
      send(server, [addr], addr, body)
      test_content(logfile) do s
        @test occursin("body::IOBuffer test", s)
      end
    end

    let  # send with body::IOStream
      mktemp() do path, io
        write(io, "body::IOStream test")
        seekstart(io)

        send(server, [addr], addr, io)
        test_content(logfile) do s
          @test occursin("body::IOStream test", s)
        end
      end
    end

    let  # AUTH PLAIN
      opts = SendOptions(username = "foo@example.org", passwd = "bar")
      body = IOBuffer("AUTH PLAIN test")
      send(server, [addr], addr, body, opts)
      test_content(logfile) do s
        @test occursin("AUTH PLAIN test", s)
      end
    end

    let
      opts = SendOptions(username = "foo@example.org", passwd = "invalid")
      body = IOBuffer("invalid password")
      @test_throws Exception send(server, [addr], addr, body, opts)
    end

    let  # multiple RCPT TO
      body = IOBuffer("multiple rcpt")
      rcpts = ["<foo@example.org>", "<bar@example.org>", "<baz@example.org>"]
      send(server, rcpts, addr, body)

      test_content(logfile) do s
        @test occursin("multiple rcpt", s)
        @test occursin("X-RCPT: foo@example.org", s)
        @test occursin("X-RCPT: bar@example.org", s)
        @test occursin("X-RCPT: baz@example.org", s)
      end
    end

    let  # non-blocking send
      body = IOBuffer("non-blocking send")
      task = @async send(server, [addr], addr, body)
      wait(task)
      test_content(logfile) do s
        @test occursin("non-blocking send", s)
      end
    end

    let  # SendOptions.verbose no error
      opts = SendOptions(verbose = true, username = "foo@example.org", passwd = "bar")
      body = IOBuffer("SendOptions.verbose")
      send(server, [addr], addr, body, opts)
      test_content(logfile) do s
        @test occursin("SendOptions.verbose", s)
      end
    end

    let  # send using get_body
      message = "body mime message"
      subject = "test message"

      mime_message = get_mime_msg(message, Val(:usascii))
      body = get_body([addr], addr, subject, mime_message)

      send(server, [addr], addr, body)

      test_content(logfile) do s
        @test occursin("From: $addr", s)
        @test occursin("To: $addr", s)
        @test occursin("Subject: $subject", s)
        @test occursin(message, s)
      end
    end

    let  # send using get_body with UTF-8 encoded message
      message =
        "ABCDEFGHIJKLMNOPQRSTUVWXYZ /0123456789\r\n" *
        "abcdefghijklmnopqrstuvwxyz £©µÀÆÖÞßéöÿ\r\n" *
        "–—‘“”„†•…‰™œŠŸž€ ΑΒΓΔΩαβγδω АБВГДабвгд\r\n" *
        "∀∂∈ℝ∧∪≡∞ ↑↗↨↻⇣ ┐┼╔╘░►☺♀ ﬁ�⑀₂ἠḂӥẄɐː⍎אԱა\r\n"
      subject = "test message in UTF-8"

      mime_message = get_mime_msg(message, Val(:utf8))
      body = get_body([addr], addr, subject, mime_message)

      send(server, [addr], addr, body)

      test_content(logfile) do s
        @test occursin("From: $addr", s)
        @test occursin("To: $addr", s)
        @test occursin("Subject: $subject", s)
        @test occursin("ABCDEFGHIJKLMNOPQRSTUVWXYZ /0123456789", s)
        @test occursin("abcdefghijklmnopqrstuvwxyz £©µÀÆÖÞßéöÿ", s)
        @test occursin("–—‘“”„†•…‰™œŠŸž€ ΑΒΓΔΩαβγδω АБВГДабвгд", s)
        @test occursin("∀∂∈ℝ∧∪≡∞ ↑↗↨↻⇣ ┐┼╔╘░►☺♀ ﬁ�⑀₂ἠḂӥẄɐː⍎אԱა", s)
      end
    end

    let  # send using get_body with extra fields
        message = "body mime message with extra fields"
        subject = "test message with extra fields"

        mime_message = get_mime_msg(message)
        from = addr
        to = ["<foo@example.org>", "<bar@example.org>"]
        cc = ["<baz@example.org>", "<qux@example.org>"]
        bcc = ["<quux@example.org>"]
        replyto = addr
        body = get_body(to, from, subject, mime_message;
            cc = cc, replyto = replyto)
        rcpts = vcat(to, cc, bcc)

        send(server, rcpts, addr, body)

        test_content(logfile) do s
          @test occursin("From: $addr", s)
          @test occursin("Subject: $subject", s)
          @test occursin("Cc: <baz@example.org>, <qux@example.org>", s)
          @test occursin("Reply-To: $addr", s)
          @test occursin("To: <foo@example.org>, <bar@example.org>", s)
          @test occursin(message, s)
          @test occursin("X-RCPT: foo@example.org", s)
          @test occursin("X-RCPT: bar@example.org", s)
          @test occursin("X-RCPT: baz@example.org", s)
          @test occursin("X-RCPT: qux@example.org", s)
          @test occursin("X-RCPT: quux@example.org", s)
      end
    end

    let  # send with attachment
      message = "body mime message with attachment"
      subject = "test message with attachment"
      svg_str = """<?xml version="1.0" encoding="UTF-8"?>
        <svg xmlns="http://www.w3.org/2000/svg" xmlns:xlink="http://www.w3.org/1999/xlink" width="320pt" height="200pt" viewBox="0 0 320 200" version="1.1">
        <g id="surface61">
        <path style=" stroke:none;fill-rule:nonzero;fill:rgb(0%,0%,0%);fill-opacity:1;" d="M 67.871094 164.3125 C 67.871094 171.847656 67.023438 177.933594 65.328125 182.566406 C 63.632812 187.203125 61.222656 190.800781 58.09375 193.363281 C 54.96875 195.925781 51.21875 197.640625 46.847656 198.507812 C 42.476562 199.371094 37.613281 199.804688 32.265625 199.804688 C 25.027344 199.804688 19.488281 198.675781 15.648438 196.414062 C 11.804688 194.152344 9.882812 191.441406 9.882812 188.273438 C 9.882812 185.636719 10.953125 183.414062 13.101562 181.605469 C 15.25 179.796875 18.132812 178.894531 21.75 178.894531 C 24.464844 178.894531 26.632812 179.628906 28.25 181.097656 C 29.871094 182.566406 31.210938 184.019531 32.265625 185.449219 C 33.46875 187.03125 34.488281 188.085938 35.316406 188.613281 C 36.144531 189.140625 36.898438 189.40625 37.578125 189.40625 C 39.007812 189.40625 40.101562 188.558594 40.855469 186.863281 C 41.609375 185.167969 41.984375 181.871094 41.984375 176.972656 L 41.984375 84.050781 L 67.871094 76.929688 L 67.871094 164.3125 M 104.738281 79.414062 L 104.738281 139.214844 C 104.738281 140.875 105.058594 142.4375 105.699219 143.90625 C 106.339844 145.375 107.226562 146.640625 108.355469 147.695312 C 109.488281 148.75 110.804688 149.597656 112.3125 150.238281 C 113.820312 150.878906 115.441406 151.199219 117.175781 151.199219 C 119.132812 151.199219 121.359375 150.101562 124.070312 148.203125 C 128.363281 145.195312 130.964844 143.128906 130.964844 140.683594 C 130.964844 140.097656 130.964844 79.414062 130.964844 79.414062 L 156.738281 79.414062 L 156.738281 164.3125 L 130.964844 164.3125 L 130.964844 156.398438 C 127.574219 159.261719 123.957031 161.558594 120.113281 163.292969 C 116.269531 165.027344 112.539062 165.894531 108.921875 165.894531 C 104.703125 165.894531 100.78125 165.195312 97.164062 163.800781 C 93.546875 162.40625 90.382812 160.503906 87.671875 158.09375 C 84.957031 155.683594 82.828125 152.855469 81.28125 149.613281 C 79.738281 146.375 78.964844 142.90625 78.964844 139.214844 L 78.964844 79.414062 L 104.738281 79.414062 M 192.882812 164.3125 L 167.222656 164.3125 L 167.222656 45.277344 L 192.882812 38.15625 L 192.882812 164.3125 M 203.601562 84.050781 L 229.375 76.929688 L 229.375 164.3125 L 203.601562 164.3125 L 203.601562 84.050781 M 283.226562 120.449219 C 280.738281 121.507812 278.230469 122.730469 275.707031 124.125 C 273.183594 125.519531 270.882812 127.046875 268.8125 128.703125 C 266.738281 130.359375 265.0625 132.132812 263.78125 134.015625 C 262.5 135.898438 261.859375 137.859375 261.859375 139.894531 C 261.859375 141.476562 262.066406 143.003906 262.480469 144.472656 C 262.894531 145.941406 263.480469 147.203125 264.234375 148.257812 C 264.988281 149.3125 265.816406 150.160156 266.722656 150.800781 C 267.625 151.441406 268.605469 151.761719 269.660156 151.761719 C 271.769531 151.761719 273.898438 151.121094 276.046875 149.839844 C 278.195312 148.558594 280.585938 146.941406 283.226562 144.980469 L 283.226562 120.449219 M 309.109375 164.3125 L 283.226562 164.3125 L 283.226562 157.527344 C 281.792969 158.734375 280.398438 159.847656 279.042969 160.863281 C 277.6875 161.878906 276.160156 162.765625 274.464844 163.519531 C 272.769531 164.273438 270.867188 164.855469 268.753906 165.273438 C 266.644531 165.6875 264.15625 165.894531 261.296875 165.894531 C 257.375 165.894531 253.851562 165.328125 250.726562 164.199219 C 247.597656 163.066406 244.941406 161.523438 242.757812 159.5625 C 240.570312 157.605469 238.894531 155.285156 237.726562 152.609375 C 236.558594 149.9375 235.972656 147.015625 235.972656 143.851562 C 235.972656 140.609375 236.59375 137.671875 237.839844 135.03125 C 239.082031 132.394531 240.777344 130.023438 242.925781 127.910156 C 245.074219 125.800781 247.578125 123.917969 250.441406 122.257812 C 253.304688 120.601562 256.378906 119.074219 259.65625 117.679688 C 262.933594 116.285156 266.34375 115.007812 269.886719 113.839844 C 273.425781 112.671875 276.933594 111.558594 280.398438 110.503906 L 283.226562 109.824219 L 283.226562 101.460938 C 283.226562 96.035156 282.1875 92.191406 280.117188 89.929688 C 278.042969 87.667969 275.273438 86.539062 271.808594 86.539062 C 267.738281 86.539062 264.910156 87.519531 263.328125 89.476562 C 261.746094 91.4375 260.953125 93.808594 260.953125 96.597656 C 260.953125 98.179688 260.785156 99.726562 260.445312 101.234375 C 260.109375 102.742188 259.523438 104.058594 258.695312 105.191406 C 257.867188 106.320312 256.679688 107.226562 255.132812 107.902344 C 253.589844 108.582031 251.648438 108.921875 249.3125 108.921875 C 245.695312 108.921875 242.757812 107.882812 240.496094 105.8125 C 238.234375 103.738281 237.105469 101.121094 237.105469 97.953125 C 237.105469 95.015625 238.101562 92.285156 240.097656 89.761719 C 242.097656 87.234375 244.789062 85.066406 248.183594 83.261719 C 251.574219 81.449219 255.492188 80.019531 259.9375 78.964844 C 264.382812 77.910156 269.09375 77.382812 274.066406 77.382812 C 280.171875 77.382812 285.429688 77.929688 289.839844 79.019531 C 294.246094 80.113281 297.882812 81.675781 300.746094 83.710938 C 303.609375 85.746094 305.71875 88.195312 307.074219 91.058594 C 308.433594 93.921875 309.109375 97.128906 309.109375 100.667969 L 309.109375 164.3125 "/>
        <path style=" stroke:none;fill-rule:nonzero;fill:rgb(79.6%,23.5%,20%);fill-opacity:1;" d="M 235.273438 55.089844 C 235.273438 64.757812 227.4375 72.589844 217.773438 72.589844 C 208.105469 72.589844 200.273438 64.757812 200.273438 55.089844 C 200.273438 45.425781 208.105469 37.589844 217.773438 37.589844 C 227.4375 37.589844 235.273438 45.425781 235.273438 55.089844 "/>
        <path style=" stroke:none;fill-rule:nonzero;fill:rgb(25.1%,38.8%,84.7%);fill-opacity:1;" d="M 72.953125 55.089844 C 72.953125 64.757812 65.117188 72.589844 55.453125 72.589844 C 45.789062 72.589844 37.953125 64.757812 37.953125 55.089844 C 37.953125 45.425781 45.789062 37.589844 55.453125 37.589844 C 65.117188 37.589844 72.953125 45.425781 72.953125 55.089844 "/>
        <path style=" stroke:none;fill-rule:nonzero;fill:rgb(58.4%,34.5%,69.8%);fill-opacity:1;" d="M 277.320312 55.089844 C 277.320312 64.757812 269.484375 72.589844 259.820312 72.589844 C 250.15625 72.589844 242.320312 64.757812 242.320312 55.089844 C 242.320312 45.425781 250.15625 37.589844 259.820312 37.589844 C 269.484375 37.589844 277.320312 45.425781 277.320312 55.089844 "/>
        <path style=" stroke:none;fill-rule:nonzero;fill:rgb(22%,59.6%,14.9%);fill-opacity:1;" d="M 256.300781 18.671875 C 256.300781 28.335938 248.464844 36.171875 238.800781 36.171875 C 229.132812 36.171875 221.300781 28.335938 221.300781 18.671875 C 221.300781 9.007812 229.132812 1.171875 238.800781 1.171875 C 248.464844 1.171875 256.300781 9.007812 256.300781 18.671875 "/>
        </g>
        </svg>
      """
      filename = joinpath(tempdir(), "julia_logo_color.svg")
      open(filename, "w") do f
        write(f, svg_str)
      end
      readme = open(f->read(f, String), joinpath("..", "README.md"))

      mime_message = get_mime_msg(message, Val(:utf8))
      attachments = [joinpath("..", "README.md"), filename]
      body = get_body([addr], addr, subject, mime_message, attachments = attachments)

      send(server, [addr], addr, body)

      test_content(logfile) do s
        m = match(r"Content-Type:\s*multipart\/mixed;\s*boundary=\"(.+)\"\n", s)
        @test m !== nothing
        boundary = m.captures[1]
        @test occursin("To: $addr", s)
        @test occursin("Subject: $subject", s)
        @test occursin(message, s)
        splt = split(s)
        ind = findall(v -> occursin("--$boundary", v), splt)
        @test length(ind) == 6
        @test String(base64decode(splt[ind[4]-1])) == readme
        @test String(base64decode(splt[ind[6]-1])) == svg_str
      end
      rm(filename)
    end

  finally
    kill(smtpsink)
    rm(logfile, force = true)
  end
end  # @testset "Send"


text:
module DifferencesBase

using SparseArrays

export AbstractDifference,
    VectorDifference,
    MatrixDifference,
    DictDifference,
    NamedTupleDifference,
    SetDifference,
    added_indices,
    modified_indices,
    removed_indices,
    added,
    common,
    modified,
    removed

include("types.jl")
include("differences.jl")
include("math.jl")
include("show.jl")

end # module


text:
# ===========================================================================
# Spacecraft Rendezvous model
# See https://easychair.org/publications/paper/S2V for a description and
# https://gitlab.com/goranf/ARCH-COMP/blob/master/2018/NLN/C2E2/ARCH18_NLN/rendezvous/c2e2_nonlinear_pass_4d.hyxml
# for a reference model
# ===========================================================================
using SparseArrays, HybridSystems, MathematicalSystems, LazySets, Reachability, TaylorIntegration
using Reachability: solve

include("eqs_mymul.jl")

function spacecraft_rendezvous(;T=200.0, orderT=10, orderQ=2, abs_tol=1e-10,
                                max_steps=500, plot_vars=[1, 2], project_reachset=false)
    # variables
    x = 1   # x position (negative!)
    y = 2   # y position (negative!)
    vx = 3  # x velocity
    vy = 4  # y velocity
    t = 5   # time

    # number of variables
    n = 4 + 1

    # time of abort
    t_abort = 120.0

    # discrete structure (graph)
    automaton = LightAutomaton(3)

    # mode 1 ("approaching")
    𝐹 = spacecraft_approaching!
    invariant = HPolyhedron([
        HalfSpace(sparsevec([x], [1.], n), -100.),   # x <= -100
        HalfSpace(sparsevec([t], [1.], n), t_abort)  # t <= t_abort
       ])
    m₁ = CBBCS(𝐹, 5, invariant)

    # mode 2 ("rendezvous attempt")
    𝐹 = spacecraft_rendezvous_attempt!
    invariant = HPolyhedron([
        HalfSpace(sparsevec([x], [-1.], n), 100.),           # x >= -100
        HalfSpace(sparsevec([x], [1.], n), 100.),            # x <= 100
        HalfSpace(sparsevec([y], [-1.], n), 100.),           # y >= -100
        HalfSpace(sparsevec([y], [1.], n), 100.),            # y <= 100
        HalfSpace(sparsevec([x, y], [-1., -1.], n), 141.1),  # x + y >= -141.1
        HalfSpace(sparsevec([x, y], [1., 1.], n), 141.1),    # x + y <= 141.1
        HalfSpace(sparsevec([x, y], [1., -1.], n), 141.1),   # -x + y >= -141.1
        HalfSpace(sparsevec([x, y], [-1., 1.], n), 141.1),   # -x + y <= 141.1
        HalfSpace(sparsevec([t], [1.], n), t_abort)          # t <= t_abort
       ])
    m₂ = CBBCS(𝐹, 5, invariant)

    # mode 3 ("aborting")
    𝐹 = spacecraft_aborting!
    invariant = Universe(n)
    m₃ = CBBCS(𝐹, 5, invariant)

    # modes
    modes = [m₁, m₂, m₃]

    # transition 1 -> 2
    add_transition!(automaton, 1, 2, 1)
    guard = HPolyhedron([
        HalfSpace(sparsevec([x], [-1.], n), 100.),           # x >= -100
        HalfSpace(sparsevec([x], [1.], n), 100.),            # x <= 100
        HalfSpace(sparsevec([y], [-1.], n), 100.),           # y >= -100
        HalfSpace(sparsevec([y], [1.], n), 100.),            # y <= 100
        HalfSpace(sparsevec([x, y], [-1., -1.], n), 141.1),  # x + y >= -141.1
        HalfSpace(sparsevec([x, y], [1., 1.], n), 141.1),    # x + y <= 141.1
        HalfSpace(sparsevec([x, y], [1., -1.], n), 141.1),   # -x + y >= -141.1
        HalfSpace(sparsevec([x, y], [-1., 1.], n), 141.1)    # -x + y <= 141.1
       ])
    t₁ = ConstrainedIdentityMap(n, guard)

    # transition 1 -> 3
    add_transition!(automaton, 1, 3, 2)
    guard = HalfSpace(sparsevec([t], [-1.], n), -t_abort)  # t >= t_abort
    t₂ = ConstrainedIdentityMap(n, guard)

    # transition 2 -> 3
    add_transition!(automaton, 2, 3, 3)
    guard = HalfSpace(sparsevec([t], [-1.], n), -t_abort)  # t >= t_abort
    t₃ = ConstrainedIdentityMap(n, guard)

    # transition annotations
    resetmaps = [t₁, t₂, t₃]

    # switching
    switchings = [AutonomousSwitching()]

    ℋ = HybridSystem(automaton, modes, resetmaps, switchings)

    X0 = Hyperrectangle([-900., -400., 0., 0., 0.], [25., 25., 0., 0., 0.])
    initial_condition = [(1, X0)]

    𝑃 = InitialValueProblem(ℋ, initial_condition)

    # safety property in "Mode 2"
    velocity = 0.055 * 60.  # meters per minute
    cx = velocity * cos(π / 8)  # x-coordinate of the octagon's first (ENE) corner
    cy = velocity * sin(π / 8)  # y-coordinate of the octagon's first (ENE) corner

    tan30 = tan(π/6)
    #=
    octagon = [
        HalfSpace(sparsevec([vx], [1.], n), cx),                 # vx <= cx
        HalfSpace(sparsevec([vx, vy], [1., 1.], n), cy + cx),    # vx + vy <= cy + cx
        HalfSpace(sparsevec([vy], [1.], n), cx),                 # vy <= cx
        HalfSpace(sparsevec([vx, vy], [-1., 1.], n), cy + cx),   # -vx + vy <= cy + cx
        HalfSpace(sparsevec([vx], [-1.], n), cx),                # vx >= -cx
        HalfSpace(sparsevec([vx, vy], [-1., -1.], n), cy + cx),  # -vx - vy <= cy + cx
        HalfSpace(sparsevec([vy], [-1.], n), cx),                # vy >= -cx
        HalfSpace(sparsevec([vx, vy], [1., -1.], n), cy + cx)    # vx - vy <= cy + cx
       ]
    cone = [
        HalfSpace(sparsevec([x], [-1.], n), 100.),          # x >= -100
        HalfSpace(sparsevec([x, y], [tan30, -1.], n), 0.),  # -x tan(30°) + y >= 0
        HalfSpace(sparsevec([x, y], [tan30, 1.], n), 0.),   # -x tan(30°) - y >= 0
       ]
    property_rendezvous = SafeStatesProperty(HPolytope([octagon; cone]))
    =#

    property_rendezvous = (t, x) -> (x[3] <= cx) && (x[3]+x[4]<=cx+cy) && (x[4]<=cx) && (-x[3]+x[4] <= cx+cy) && (x[3] >= -cx) &&
                                    (-x[3]-x[4] <= cx+cy) && (x[4] >= -cx) && (x[3]-x[4] <= cx+cy) &&
                                    (x[1] >= -100) && (-x[1]*tan30+x[2]>=0) && (-x[1]*tan30-x[2]>=0)

    # safety property in "Passive"
    #=
    target = HPolytope([
        HalfSpace(sparsevec([x], [1.], n), 0.2),   # x <= 0.2
        HalfSpace(sparsevec([x], [-1.], n), 0.2),  # x >= -0.2
        HalfSpace(sparsevec([y], [1.], n), 0.2),   # y <= 0.2
        HalfSpace(sparsevec([y], [-1.], n), 0.2),  # y >= -0.2
       ])
    property_aborting = BadStatesProperty(target)
    =#

    property_aborting = (t, x) -> !( (x[1] <= 0.2) && (x[1] >= -0.2) && (x[2] <= 0.2) && (x[2] >= -0.2))

    # safety properties
    property = Dict{Int, Function}(1 => (t, x) -> true,
                                   2 => property_rendezvous,
                                   3 => property_aborting)

    # global options
    𝑂 = Options(:T=>T, :property=>property, :plot_vars=>plot_vars,
                :project_reachset=>project_reachset, :mode=>"reach")

    # algorithm-specific options
    #𝑂jets = Options(:orderT=>orderT, :orderQ=>orderQ, :abs_tol=>abs_tol, :max_steps=>max_steps)

    return 𝑃, 𝑂
end

max_stars_repo_name: Aadesha/ReachabilityBenchmarks

max_issues_repo_name: Aadesha/ReachabilityBenchmarks

max_forks_repo_name: Aadesha/ReachabilityBenchmarks


text:
module JGraph

using Distributions
using Random
using DataStructures

# export
export  AbstractGraph,
        Graph,
        nv,
        ne,
        has_edge,
        add_edge!,
        add_vertex!,
        remove_edge!,
        remove_vertex!,
        degree,
        neighbors,
        vertices,
        edges

export  Edge,
        dst,
        src,
        EdgeIter

export  erdos_renyi

export  bfs,
        connected_components

# Include
include("edge.jl")
include("edgeiter.jl")
include("core.jl")
include("generator.jl")
include("traversal.jl")
end # module


text:
################################################################################
# Fundamental "neural" data types
################################################################################
struct Spike
	n::UInt64
	t::Float64
end
Base.show(io::IO, s::Spike) = print(io, "(Neuron $(s.n): $(s.t) s)")

struct Synapse
	strength::Float64
	source::UInt64
	dest::UInt64
	delay::UInt16 	# in number of samples
end
Base.show(io::IO, syn::Synapse) = print(io, "(Source: $(syn.source), Destination: $(syn.dest), Delay: $(syn.delay), Strength: $(syn.strength))")

abstract type Neuron end

################################################################################
# Izhikevich neuron
################################################################################
struct SimpleNeuron <: Neuron
	a::Float64
	b::Float64
	c::Float64
	d::Float64
end

struct SimpleNeuronType{T} end
SimpleNeuronType(s::Symbol) = SimpleNeuronType{s}()
SimpleNeuronType(s::AbstractString) = SimpleNeuronType{Symbol(s)}()

SimpleNeuron(s::Symbol) = SimpleNeuron(SimpleNeuronType{Symbol(s)}())
SimpleNeuron(s::AbstractString) = SimpleNeuron(SimpleNeuronType(s))

# Regular Spiking (Exc)
SimpleNeuron(::SimpleNeuronType{:rs}) = SimpleNeuron(0.02, 0.2, -65.0, 8.0)
SimpleNeuron(::SimpleNeuronType{:RS}) = SimpleNeuron(0.02, 0.2, -65.0, 8.0)
# Intrinsically Bursting (Exc)
SimpleNeuron(::SimpleNeuronType{:ib}) = SimpleNeuron(0.02, 0.2, -55.0, 4.0)
SimpleNeuron(::SimpleNeuronType{:IB}) = SimpleNeuron(0.02, 0.2, -55.0, 4.0)
# Chattering (Exc)
SimpleNeuron(::SimpleNeuronType{:ch}) = SimpleNeuron(0.02, 0.2, -50.0, 2.0)
SimpleNeuron(::SimpleNeuronType{:CH}) = SimpleNeuron(0.02, 0.2, -50.0, 2.0)
# Fast Spiking (Inh)
SimpleNeuron(::SimpleNeuronType{:fs}) = SimpleNeuron(0.1, 0.2, -65.0, 2.0)
SimpleNeuron(::SimpleNeuronType{:FS}) = SimpleNeuron(0.1, 0.2, -65.0, 2.0)
# Low-Threshold Spiking (Inh)
SimpleNeuron(::SimpleNeuronType{:lts}) = SimpleNeuron(0.02, 0.25, -65.0, 2.0)
SimpleNeuron(::SimpleNeuronType{:LTS}) = SimpleNeuron(0.02, 0.25, -65.0, 2.0)
# Thalamo-Cortical
SimpleNeuron(::SimpleNeuronType{:tc}) = SimpleNeuron(0.02, 0.25, -65.0, 0.05)
SimpleNeuron(::SimpleNeuronType{:TC}) = SimpleNeuron(0.02, 0.25, -65.0, 0.05)
# Resonator
SimpleNeuron(::SimpleNeuronType{:rz}) = SimpleNeuron(0.1, 0.26, -65.0, 2.0)
SimpleNeuron(::SimpleNeuronType{:RZ}) = SimpleNeuron(0.1, 0.26, -65.0, 2.0)

# Random Excitatory
SimpleNeuron(::SimpleNeuronType{:exc}) = begin
	c = rand()^2.0
	SimpleNeuron(0.02, 0.2, -65.0+15.0c, 8.0-6.0c)
end
SimpleNeuron(::SimpleNeuronType{:EXC}) = SimpleNeuron(SimpleNeuronType(:exc))
# Random Inhibitory
SimpleNeuron(::SimpleNeuronType{:inh}) = begin
	r = rand()
	SimpleNeuron(0.02+0.08r, 0.25-0.05r, -65.0, 2.0)
end
SimpleNeuron(::SimpleNeuronType{:INH}) = SimpleNeuron(SimpleNeuronType(:inh))

################################################################################
# Gerstner neuron
################################################################################
mutable struct SRMNeuron <: Neuron
	u::Float64
	u_max::Float64
	u_rest::Float64
	Θ::Float64
	τ_refr::Float64
	τ_memb::Float64
	t_lastevent::Float64
	t_lastspike::Float64
	ϵ_last::Float64
end

SRMNeuron() = SRMNeuron(-65.0, 8.0, -65.0, -50.0, 7.0/1000.0, 3.0/1000.0,
						0.0, 0.0, 0.0)

SRMNeuron(τ_membrane) = SRMNeuron(-65.0, 8.0, -65.0, -50.0, 7.0/1000.0, τ_membrane,
								  0.0, 0.0, 0.0)

######################################################################
# Delay lines
######################################################################

mutable struct DelayLine
	numstored::Int64
	delaylen::Int64
	counts::Vector{Int64}
end

DelayLine(delaylen::Int64, maxevents::Int64) =
	DelayLine(0, delaylen, zeros(Int64, maxevents))

DelayLine(delaylen::Int64) = DelayLine(0, delaylen, zeros(Int64, 64))

DelayLine(numstored::Int64, delaylen::Int64, maxevents::Int64) =
	DelayLine(numstored, delaylen, zeros(Int64, maxevents))

################################################################################
# Types for managing experiments
################################################################################
struct ModelParams
	fs::Float64

	# Graph (revise later) -- now specified directly, but still used here
	num_neurons::UInt64
	p_contact::Float64
	p_exc::Float64
	maxdelay::Float64 # ms

	# STDP
	synmax::Float64
	tau_pre::Float64
	tau_post::Float64
	a_pre::Float64
	a_post::Float64
end
ModelParams() = ModelParams(2000.0, 1000, 0.1, 0.8, 20.0, 10.0,
							0.02, 0.02, 1.20, 1.0)

struct TrialParams
	dur::Float64
	lambda::Float64
	randspikesize::Float64
	randinput::Bool
	inhibition::Bool
    stdp::Bool
	inputmode::UInt64
	multiinputmode::UInt64
	inputweight::Float64
	recordstart::Float64
	recordstop::Float64
	lambdainput::Float64
	inputrefractorytime::Float64
end

TrialParams() = TrialParams(1.0, 3.0, 20.0, 1, 1, 1, 1, 1, 20.0, 0.0, 1.0, 0.5, 0.0)

function Base.show(io::IO, tp::TrialParams)
	println(io, "\tInput mode: $(tp.inputmode == 1 ? "periodic" : "poissonian")")
	println(io, "\tInput timing density (if random): $(tp.lambdainput) (1/s)")
	println(io, "\tInput size: $(tp.inputweight)")
	println(io, "\tDuration: $(tp.dur) (s)")
	println(io, "\tRandom input: $(tp.randinput)")
	println(io, "\tRand spike density (λ): $(tp.lambda) (1/s)")
	println(io, "\tRandom spike size: $(tp.randspikesize)")
	println(io, "\tInhibition: $(tp.inhibition)")
	println(io, "\tSTDP: $(tp.stdp)")
	println(io, "\tRecord start: $(tp.recordstart)")
	println(io, "\tRecord stop: $(tp.recordstop)")
end

struct Experiment
	name::String
	mp::ModelParams
	tp::TrialParams
	input::Array{Spike,1}
end

struct Model
	params::ModelParams
	delgraph::Array{Int64,2}
	syngraph::Array{Float64,2}
	neurons::Array{SimpleNeuron,1}
end

struct ExperimentOutput
	modelname::AbstractString
	trialname::AbstractString
	tp::TrialParams
	input::Array{Array{Spike,1},1}
	output::Array{Spike,1}
	inputtimes::Array{Float64,1}
	inputids::Array{Int64,1}
	synapses::Array{Synapse,1}
end

function Base.show(io::IO, exp::ExperimentOutput)
	println(IO, "Model Name: $(exp.modelname)")
	println(IO, "Trial Name: $(exp.trialname)")
	println(IO, "Trial Parameters")
	Base.show(io, exp.tp)
end

struct ExperimentType{T} end

ExperimentType(s::AbstractString) = ExperimentType{Symbol(s)}()

################################################################################
# Delay learning types
################################################################################
struct CategoricalSpiker
	dls::Array{DelayLine,1}
	neuron::SRMNeuron
	weight::Float64
end

CategoricalSpiker(inputs::T, possibledelays::T, weight;
				  τ_membrane=3.0/1000.0,
				  maxevents=64) where T <: Union{Array{Int64,1}, UnitRange{Int64}} = begin
	dls = [DelayLine(rand(possibledelays), maxevents) for _ ∈ inputs]
	neuron = SRMNeuron(τ_membrane)
	CategoricalSpiker(dls, neuron, weight)
end

max_stars_repo_name: analogouscircuit/DelayNetTools

max_issues_repo_name: analogouscircuit/DelayNetTools

max_forks_repo_name: analogouscircuit/DelayNetTools


text:
# tests for homotopy.jl

const test_homotopy_inputfile = "input_vals_channel.jl"
const test_homotopy_moddict = Dict{String, Any}(
  "Flux_name" => "RoeFlux",
  "use_DG" => true,
  "IC_name" => "ICFreeStream",
  "BC1_name" => "FreeStreamBC",
  "jac_type" => 1,
  "jac_method" => 2,
  "new_fname" => "input_vals_channel_dg_homotopy")

function test_homotopy(mesh, sbp, eqn, opts)

  @testset "----- Testing Homotopy operators ------" begin
  # the initial condition is uniform flow, so the residual of the homotopy
  # should be zero

  res = zeros(eqn.res)
  fill!(eqn.res, 42)  # make sure eqn.res is not modified

  EulerEquationMod.calcHomotopyDiss(mesh, sbp, eqn, opts, res)

  for i=1:mesh.numEl
    @test isapprox( norm(res[:, :, i]), 0.0) atol=1e-12
  end

  for i=1:length(eqn.res)
    @test ( eqn.res[i] )== 42
  end

  # make homotopy look like a physics
  function homotopy_physics_test(mesh, sbp, eqn, opts, t=0.0)
    evalHomotopy(mesh, sbp, eqn, opts, eqn.res)
  end

  ctx_residual = (homotopy_physics_test,)

  # test jacobian
  t = 0.0
  EulerEquationMod.ICExp(mesh, sbp, eqn, opts, eqn.q_vec)
  array1DTo3D(mesh, sbp, eqn, opts, eqn.q_vec, eqn.q)

  pc, lo = NonlinearSolvers.getHomotopyPCandLO(mesh, sbp, eqn, opts)
  ls_dense = StandardLinearSolver(pc, lo, eqn.comm, opts)
#  newton_data_dense, rhs_vec = NonlinearSolvers.setupNewton(mesh, mesh, sbp, eqn, opts, ls_dense)

  opts2 = copy(opts)
  opts2["jac_type"] = 2  # SparseMatrixCSC
  pc, lo = NonlinearSolvers.getHomotopyPCandLO(mesh, sbp, eqn, opts2)
  ls_sparse = StandardLinearSolver(pc, lo, eqn.comm, opts)

#  newton_data_sparse, rhs_vec_sparse = NonlinearSolvers.setupNewton(mesh, mesh, sbp, eqn, opts2, ls_sparse)

  lo2_dense = getBaseLO(ls_dense.lo)
  lo2_sparse = getBaseLO(ls_sparse.lo)
  @test ( typeof(lo2_dense.A) <: Array )== true
  @test ( typeof(lo2_sparse.A) <: SparseMatrixCSC )==  true

  calcLinearOperator(ls_dense.lo, mesh, sbp, eqn, opts, ctx_residual, t)
  calcLinearOperator(ls_sparse.lo, mesh, sbp, eqn, opts2, ctx_residual, t)
#  NonlinearSolvers.physicsJac(newton_data_dense, mesh, sbp, eqn, opts, jac, ctx_residual)
#  NonlinearSolvers.physicsJac(newton_data_sparse, mesh, sbp, eqn, opts2, jac_sparse, ctx_residual)

  jac_dense = lo2_dense.A
  jac_dense2 = full(lo2_sparse.A)

  for i=1:mesh.numDof
    for j=1:mesh.numDof
      @test isapprox( jac_dense2[j, i], jac_dense[j, i]) atol=1e-12
    end
  end

  #test matrix-free products
  opts2["jac_type"] = 4
  println("constructing mat-free linear operator")
  pc_free, lo_free = NonlinearSolvers.getHomotopyPCandLO(mesh, sbp, eqn, opts2)
  println("typeof(lo_free) = ", typeof(lo_free))
  println("calculating linear operaotr")
  calcLinearOperator(lo_free, mesh, sbp, eqn, opts2, ctx_residual, t)
  x = rand(mesh.numDof)
  b = zeros(x)
  b2 = zeros(x)

  println("applying linear operator")
  applyLinearOperator(ls_dense.lo, mesh, sbp, eqn, opts, ctx_residual, t, x, b)
  applyLinearOperator(lo_free, mesh, sbp, eqn, opts2, ctx_residual, t, x, b2)

  @test isapprox( norm(b - b2), 0.0) atol=1e-13

#=
  # explicit jacobian computation
  println("testing explicit jacobian calculation")
  opts2["jac_type"] = 2
  pc, lo = NonlinearSolvers.getHomotopyPCandLO(mesh, sbp, eqn, opts)
  pc, lo2 = NonlinearSolvers.getHomotopyPCandLO(mesh, sbp, eqn, opts)

  calcLinearOperator(lo, mesh, sbp, eqn, opts, ctx_residual, t)
  opts2["calc_jac_explicit"] = true
  calcLinearOperator(lo2, mesh, sbp, eqn, opts2, ctx_residual, t)

  A = getBaseLO(lo).A
  A2 = getBaseLO(lo2).A
  @test isapprox( norm(full(A) - full(A2)), 0.0) atol=1e-13
=#

end  # end do
  return nothing
end

add_func3!(EulerTests, test_homotopy, test_homotopy_inputfile, test_homotopy_moddict, [TAG_HOMOTOPY, TAG_SHORTTEST])

function test_homotopy_convergence()

  println("\ntesting homotopy jacobians")
  # run for 2 iterations, with explicit, coloring jacobians, make sure
  # residuals are the same
  opts_tmp = read_input_file("input_vals_homotopy.jl")
  opts_tmp["calc_jac_explicit"] = true
  opts_tmp["itermax"] = 2
  make_input(opts_tmp, "input_vals_homotopy_tmp.jl")
  opts_tmp["calc_jac_explicit"] = false
  make_input(opts_tmp, "input_vals_homotopy_tmp2.jl")

  mesh, sbp, eqn, opts = run_solver("input_vals_homotopy_tmp.jl")
  mesh2, sbp2, eqn2, opts2 = run_solver("input_vals_homotopy_tmp2.jl")

  @test isapprox( calcNorm(eqn, eqn.q), calcNorm(eqn2, eqn2.q)) atol=1e-13

  # now run a full case
  println("\nTesting homotopy")
  rmfile("convergence.dat")
  mesh, sbp, eqn, opts = run_solver("input_vals_homotopy.jl")

  @test  calcNorm(eqn, eqn.res_vec, strongres=true)  < opts["res_abstol"]
  data = readdlm("convergence.dat")

  @test  size(data, 1)  < 30# 16 iterations to convergence

  return nothing
end

add_func1!(EulerTests, test_homotopy_convergence, [TAG_HOMOTOPY, TAG_LONGTEST])


text:
export
    # From libim.jl
    imFileOpen,
    imFileOpenAs,
    imFileNew,
    imFileClose,
    imFileHandle,
    imFileGetInfo,
    imFileSetInfo,
    imFileSetAttribute,
    imFileGetAttribute,
    imFileGetAttributeList,
    imFileGetPalette,
    imFileSetPalette,
    imFileReadImageInfo,
    imFileWriteImageInfo,
    imFileReadImageData,
    imFileWriteImageData,
    imFormatRegisterInternal,
    imFormatRemoveAll,
    imFormatList,
    imFormatInfo,
    imFormatInfoExtra,
    imFormatCompressions,
    imFormatCanWriteImage,
    imEncodeColor,
    imDecodeColor,
    imFileFormat,
    imImageInfo,
    imLoadRGB,
    imSaveRGB,
    imLoadMap,
    imSaveMap,
    imRGB2Map,
    imMap2RGB,
    imRGB2Gray,
    imMap2Gray,
    imResize,
    imStretch,
    imRegisterCallback,
    #
    # From im_image.jl        -- A submodule perhaps?
    #
    imImageCreate,
    imImageInit,
    imImageCreateBased,
    imImageDestroy,
    imImageAddAlpha,
    imImageSetAlpha,
    imImageRemoveAlpha,
    imImageReshape,
    imImageCopy,
    imImageCopyData,
    imImageCopyAttributes,
    imImageMergeAttributes,
    imImageCopyPlane,
    imImageDuplicate,
    imImageClone,
    imImageSetAttribute,
    imImageGetAttribute,
    imImageGetAttributeList,
    imImageClear,
    imImageIsBitmap,
    imImageSetPalette,
    imImageMatchSize,
    imImageMatchColor,
    imImageMatchDataType,
    imImageMatchColorSpace,
    imImageMatch,
    imImageSetMap,
    imImageSetBinary,
    imImageSetGray,
    imImageMakeBinary,
    imImageMakeGray,
    imFileLoadImage,
    imFileLoadImageFrame,
    imFileLoadBitmap,
    imFileLoadImageRegion,
    imFileLoadBitmapFrame,
    imFileSaveImage,
    imFileImageLoad,
    imFileImageLoadBitmap,
    imFileImageLoadRegion,
    imFileImageSave,
    #
    # From libim_h.jl
    #
    IM_INTERRUPTED,
    IM_ALL,
    IM_COUNTER_CB,
    IM_RESOLUTION_CB,
    IM_GIF_TRANSPARENT_COLOR_CB,
    IM_TIF_IMAGE_DESCRIPTION_CB,
    IM_BYTE,
    IM_SHORT,
    IM_USHORT,
    IM_INT,
    IM_FLOAT,
    IM_CFLOAT,
    IM_RGB,
    IM_MAP,
    IM_GRAY,
    IM_BINARY,
    IM_CMYK,
    IM_YCBCR,
    IM_LAB,
    IM_LUV,
    IM_XYZ,
    IM_ALPHA,
    IM_PACKED,
    IM_TOPDOWN,
    IM_ERR_NONE,
    IM_ERR_OPEN,
    IM_ERR_ACCESS,
    IM_ERR_FORMAT,
    IM_ERR_DATA,
    IM_ERR_COMPRESS,
    IM_ERR_MEM,
    IM_ERR_COUNTER,
    IM_ERR_READ,
    IM_ERR_WRITE,
    IM_ERR_TYPE,
    IM_ERR_COMP,
    IM_BMP,
    IM_PCX,
    IM_GIF,
    IM_TIF,
    IM_RAS,
    IM_SGI,
    IM_JPG,
    IM_LED,
    IM_TGA,
    IM_NONE,
    IM_DEFAULT,
    IM_COMPRESSED,
    IM_RES_NONE,
    IM_RES_DPI,
    IM_RES_DPC,
    # From iupim.jl
    IupLoadImage, IupSaveImage, IupGetNativeHandleImage, IupGetImageNativeHandle, IupImageFromImImage,
    # From im_image_h.jl
    imImage    # A composite type

include("libim.jl")
include("im_image.jl")
include("iupim.jl")
include("im_convert.jl")
include("im_process.jl")


text:
using XKCD
using Test

@testset "XKCD.jl" begin
    c = XKCD.Comic(999)
    @test c.json.num == 999
    @test c.json.title == "Cougars"
    @test c.json.img == "https://imgs.xkcd.com/comics/cougars.png"
end


text:
## Stationary problems
function stationary_gmp(MP::MarkovProcess, order::Int, solver)
    test_fxns = polynomial.(monomials(MP.x, 1:order))
    gmp = GMPModel(solver)
    set_approximation_mode(gmp, PRIMAL_RELAXATION_MODE())
    @variable(gmp, P∞, Meas(MP.x, support = MP.X))
    @constraint(gmp, dynamics[ϕ in test_fxns], Mom(inf_generator(MP, ϕ), P∞) == 0.0)
    @constraint(gmp, normalization, Mom(1, P∞) == 1.0)
    return gmp
end

stationary_gmp(MP::ReactionProcess, order::Int, solver) = stationary_gmp(MP.JumpProcess, order, solver)

function stationary_moment(gmp, x)
    P = gmp[:P∞]
    @objective(gmp, Min, Mom(polynomial(x), P))
    optimize!(gmp)
    lb, stat_lb, time_lb = extract_solution(gmp)
    @objective(gmp, Max, Mom(polynomial(x), P))
    optimize!(gmp)
    ub, stat_ub, time_ub = extract_solution(gmp)
    return [lb, ub], [stat_lb, stat_ub], [time_lb, time_ub]
end

function prepare_reaction_system(rn::ReactionSystem, x0::Dict, solver, scales::Dict, auto_scaling)
    P = ReactionProcess(rn)
    project_into_subspace!(P, [x0[s] for s in species(rn)])
    if auto_scaling
        scales = stoich_bounds(rn, x0, solver)
    end
    x_scale = [scales[P.state_to_species[x]] for x in P.JumpProcess.x]
    x0 = [x0[P.state_to_species[x]]/scales[P.state_to_species[x]] for x in P.JumpProcess.x]
    rescale_state!(P, x_scale)
    return P, x0
end

function stationary_mean(MP::MarkovProcess, states::AbstractVector, order::Int, solver; G = [], f = [])
    if !isempty(G)
        project_into_subspace!(MP, G, f)
    end
    gmp = stationary_gmp(MP, order, solver)
    res = Dict(s => stationary_moment(gmp, s) for s in states)
    sol = Bounds((:mean, Inf), order, Dict(s => res[s][1] for s in states), InfoData(Dict(s => res[s][2] for s in states), Dict(s => res[s][3] for s in states)))
    return sol
end

function stationary_mean(rn::ReactionSystem, specs::AbstractVector, x0::Dict, order::Int, solver;
                         scales = Dict(spec => 1.0 for spec in species(rn)), auto_scaling = false)
    P, x0 = prepare_reaction_system(rn, x0, solver, scales, auto_scaling)
    sol = stationary_mean(P.JumpProcess, [P.species_to_state[s] for s in specs], order, solver)
    remap_states!(P, sol)
    return sol
end

function transient_gmp(MP::MarkovProcess, x0::AbstractVector, order::Int, trange::AbstractVector, solver)
    nt = length(trange)
    Δt = [trange[1], [trange[i] - trange[i-1] for i in 2:nt]...]
    @polyvar(t)
    z = [MP.x..., t]
    test_fxns = polynomial.(monomials(z, 0:order))
    gmp = GMPModel(solver)
    set_approximation_mode(gmp, PRIMAL_RELAXATION_MODE())
    @variable(gmp, Pₜ[i in 1:nt], Meas(MP.x, support = MP.X))
    @variable(gmp, ξₜ[i in 1:nt], Meas(z, support = intersect(MP.X, @set(t >= 0), @set(Δt[i] - t >= 0))))
    @constraint(gmp, dynamics[ϕ in test_fxns, i in 1:nt],
                     louiville(MP,Pₜ[i],ξₜ[i],ϕ,t,Δt[i]) == (i == 1 ? split_poly(ϕ,MP.x,0)(x0...) : Mom(split_poly(ϕ, MP.x, 0), Pₜ[i-1])))
    return gmp
end

transient_gmp(MP::ReactionProcess, x0::AbstractVector, order::Int, trange::AbstractVector, solver) = transient_gmp(MP.JumpProcess, x0, order, trange, solver)

function transient_moment(gmp, x, idx)
    P = gmp[:Pₜ][idx]
    @objective(gmp, Min, Mom(polynomial(x), P))
    optimize!(gmp)
    lb, stat_lb, time_lb = extract_solution(gmp)
    @objective(gmp, Max, Mom(polynomial(x), P))
    optimize!(gmp)
    ub, stat_ub, time_ub = extract_solution(gmp)
    return [lb, ub], [stat_lb, stat_ub], [time_lb, time_ub]
end

function transient_mean(MP::MarkovProcess, states::AbstractVector, x0::AbstractVector, order::Int, trange::AbstractVector, solver; G = [], f = [], t = 0)
    if !isempty(G)
        project_into_subspace!(MP, G, f)
    end
    if t != 0
        @assert t in trange "t not an element of trange"
        idx = findfirst(isequal(t), trange)
    else
        idx = length(trange)
    end
    gmp = transient_gmp(MP, x0, order, trange, solver)
    res = Dict(s => transient_moment(gmp, s, idx) for s in states)
    sol = Bounds((:mean, trange[idx]), order, Dict(s => res[s][1] for s in states), InfoData(Dict(s => res[s][2] for s in states), Dict(s => res[s][3] for s in states)))
    return sol
end

function transient_mean(rn::ReactionSystem, specs::AbstractVector, x0::Dict, order::Int, trange::AbstractVector, solver;
                        scales = Dict(spec => 1.0 for spec in species(rn)), auto_scaling = false, t=0)
    P, x0 = prepare_reaction_system(rn, x0, solver, scales, auto_scaling)
    sol = transient_mean(P.JumpProcess, [P.species_to_state[s] for s in specs], x0, order, trange, solver, t=t)
    remap_states!(P, sol)
    return sol
end

# variance
function stationary_variance(MP::MarkovProcess, states::AbstractVector, order::Int, solver; G = [], f = [])
    if !isempty(G)
        project_into_subspace!(MP, G, f)
    end
    gmp = stationary_gmp(MP, order, solver)
    model, measures = build_approximate_model(gmp, approximation_mode(gmp))
    set_optimizer(model, solver)
    P = measures[1]
    res = Dict()
    for s in states
        slack = @variable(model)
        @constraint(model, [MM.expectation(s^2,P)-slack MM.expectation(s,P);
                            MM.expectation(s,P) 1] in PSDCone())
        @objective(model, Max, slack)
        optimize!(model)
        res[s] = extract_solution(model)
    end
    sol = Bounds((:variance, Inf), order, Dict(s => res[s][1] for s in states), InfoData(Dict(s => res[s][2] for s in states), Dict(s => res[s][3] for s in states)))
    return sol
end

function stationary_variance(rn::ReactionSystem, specs::AbstractVector, x0::Dict, order::Int, solver;
                             scales = Dict(spec => 1.0 for spec in species(rn)), auto_scaling = false)
    P, x0 = prepare_reaction_system(rn, x0, solver, scales, auto_scaling)
    sol = stationary_variance(P.JumpProcess, [P.species_to_state[s] for s in specs], order, solver)
    remap_states!(P, sol)
    return sol
end

function transient_variance(MP::MarkovProcess, states::AbstractVector, x0::AbstractVector, order::Int, trange::AbstractVector, solver; G = [], f = [], t=0)
    if !isempty(G)
        project_into_subspace!(MP, G, f)
    end
    if t != 0
        @assert t in trange "t not an element of trange"
        idx = findfirst(isequal(t), trange)
    else
        idx = length(trange)
    end
    gmp = transient_gmp(MP, x0, order, trange, solver)
    model, measures = build_approximate_model(gmp, approximation_mode(gmp))
    set_optimizer(model, solver)
    P = measures[idx]
    res = Dict()
    for s in states
        slack = @variable(model)
        @constraint(model, [MM.expectation(s^2,P)-slack MM.expectation(s,P);
                            MM.expectation(s,P) 1] in PSDCone())
        @objective(model, Max, slack)
        optimize!(model)
        res[s] = extract_solution(model)
    end
    sol = Bounds((:variance, trange[idx]), order, Dict(s => res[s][1] for s in states), InfoData(Dict(s => res[s][2] for s in states), Dict(s => res[s][3] for s in states)))
    return sol
end

function transient_variance(rn::ReactionSystem, specs::AbstractVector, x0::Dict, order::Int, trange::AbstractVector, solver;
                            scales = Dict(spec => 1.0 for spec in species(rn)), auto_scaling = false, t=0)
    P, x0 = prepare_reaction_system(rn, x0, solver, scales, auto_scaling)
    sol = transient_variance(P.JumpProcess, [P.species_to_state[s] for s in specs], x0, order, trange, solver, t=t)
    remap_states!(P, sol)
    return sol
end

# Probability mass
function stationary_probability_mass(event, P::MarkovProcess, order::Int, solver)
    test_fxns = polynomial.(monomials(MP.x, 1:order))
    gmp = GMPModel(solver)
    set_approximation_mode(gmp, PRIMAL_RELAXATION_MODE())
    @variable(gmp, P∞, Meas(P.x, support = P.X))
    @variable(gmp, Q∞, Meas(P.x, support = event))
    @constraint(gmp, dynamics[ϕ in test_fxns], Mom(inf_generator(MP, ϕ), P∞ + Q∞) == 0.0)
    @constraint(gmp, normalization, Mom(1, P∞ + Q∞) == 1.0)
    @objective(gmp, Max, Mom(1, Q∞))
    optimize!(m)
    ub, stat_ub, time_ub = extract_solution(gmp)
    @objective(gmp, Min, Mom(1,ν∞))
    optimize!(m)
    lb, stat_lb, time_lb = extract_solution(gmp)
    return Bounds(:probability_mass, order, Dict(event => [lb,ub]), InfoData([stat_lb,stat_ub], [time_lb,time_ub]))
end

# volume of stationary ellipsoid
function stationary_confidence_ellipsoid(MP::MarkovProcess, states, order::Int, solver)
    # vol of confidence ellipsoide ~ det(𝔼[xxᵀ])
    gmp = stationary_gmp(MP, order, solver)
    m, measures = build_approximate_model(gmp, approximation_mode(gmp))
    set_optimizer(m, solver)
    P = measures[1]
    n = length(states)
    @variable(m, z[1:div(n*(n+1),2)])
    @variable(m, G[1:n,1:n])
    Z = UpperDiag(z, n)
    U = MM.expectation.(states*states', P)
    v = MM.expectation.(states, P)
    vᵀ = transpose(v)
    d = Diagonal(diag(Z))
    Zᵀ = transpose(Z)
    @constraint(m, Schur_Cov, [U .- G v;
                               vᵀ 1] in PSDCone())
    @constraint(m, Schur_logdet, [G  Z;
                                  Zᵀ d] in PSDCone())
    @variable(m, t[1:n])
    @constraint(m, logsum[i in 1:n], [t[i], 1, Z[i,i]] in MOI.ExponentialCone())
    @objective(m, Max, sum(t[i] for i in 1:n))
    optimize!(m)
    ub, stat, time = extract_solution(m)
    return Bounds(:volume_confidence_ellipsoid, order, Dict(states => exp(ub)), InfoData(stat, time))
end

function stationary_confidence_ellipsoid(rn::ReactionSystem, x0::Dict, species::AbstractVector, order::Int, solver;
                                         scales = Dict(spec => 1.0 for spec in species(rn)), auto_scaling = false)
    P, x0 = prepare_reaction_system(rn, x0, solver, scales, auto_scaling)
    states = [P.species_to_state[s] for s in species]
    bnds = stationary_confidence_ellipsoid(P.JumpProcess, states, order, solver)
    return Bounds(bnds.type, order, Dict(species => bnds.bounds[states]), bnds.info)
end

## Control Processes

function add_PathChanceConstraint!(gmp, cc, test_fxns, Δt, CP)
    MP = CP.MarkovProcess
    nt = length(Δt)
    ∂X = ∂(cc.X)
    Qₜ = @variable(gmp, [i in 1:nt], Meas(MP.x, support = cc.X)) # residence measure
    Rₜ = @variable(gmp, [k in 1:length(∂X), i in 1:nt],
                        Meas([MP.x..., CP.u..., CP.t], support = intersect(∂X[k], @set(CP.t >= 0), @set(1 - CP.t >= 0)))) # exit measure
    Sₜ = @variable(gmp, [i in 1:nt], Meas([MP.x..., CP.u..., CP.t], support = intersect(CP.ChanceConstraint.X, CP.U, @set(CP.t >= 0), @set(1 - CP.t >= 0))))
    @constraint(gmp, [ϕ in test_fxns, i in 1:nt],
                     louiville(MP,Qₜ[i],Sₜ[i],ϕ,CP.t,1,scale=Δt[i]) + sum(Mom(ϕ, Rₜ[k,i]) for k in length(∂X)) ==
                     (i == 1 ? split_poly(ϕ,MP.x,0)(x0...) : Mom(split_poly(ϕ, MP.x, 0), Qₜ[i-1])))
    @constraint(gmp, [1], Mom(1, Qₜ[nt]) >= 1 - cc.α)
    Tₜ = @variable(gmp, [i in 1:nt], Meas(MP.x, support = MP.X)) # slack measure
    @constraint(gmp, [i in 1:nt], Tₜ[i] + Qₜ[i] == gmp[:Pₜ][i])
end

function add_TerminalChanceConstraint!(gmp, cc, nt, CP)
    MP = CP.MarkovProcess
    Q = @variable(gmp, [nt], Meas(MP.x, support = cc.X))
    S = @variable(gmp, [nt], Meas(MP.x, support = MP.X))
    @constraint(gmp, [nt], Q[nt] + S[nt] == gmp[:Pₜ][nt])
    @constraint(gmp, [nt], Mom(1, Q[nt]) >= 1-cc.α)
end

function control_gmp(CP::ControlProcess, x0::AbstractVector, order::Int, trange::AbstractVector, solver)
    MP = CP.MP
    nt = length(trange)
    Δt = [trange[1], [trange[i] - trange[i-1] for i in 2:nt]...]
    gmp = GMPModel(solver)
    set_approximation_mode(gmp, PRIMAL_RELAXATION_MODE())
    test_fxns = polynomial.(monomials([MP.x..., CP.t], 0:order))
    if typeof(CP.Objective) == LagrangeMayer || typeof(CP.Objective) == TerminalSetProbability
        @variable(gmp, Pₜ[i in 1:nt], Meas(MP.x, support = MP.X))
        @variable(gmp, ξₜ[i in 1:nt], Meas([MP.x..., CP.u..., CP.t], support = intersect(MP.X, CP.U, @set(CP.t >= 0), @set(1 - CP.t >= 0))))
        @constraint(gmp, dynamics[ϕ in test_fxns, i in 1:nt],
                         louiville(MP,Pₜ[i],ξₜ[i],ϕ,CP.t,1,scale=Δt[i]) == (i == 1 ? split_poly(ϕ,MP.x,0)(x0...) : Mom(split_poly(ϕ, MP.x, 0), Pₜ[i-1])))

        if !isnothing(CP.PathChanceConstraints)
            for cc in CP.PathChanceConstraints
                add_PathChanceConstraint(gmp, cc, test_fxns, Δt, CP)
            end
        end

        if !isnothing(CP.TerminalChanceConstraints)
            for cc in CP.TerminalChanceConstraints
                add_TerminalChanceConstraint!(gmp, cc, nt, CP)
            end
        end

        if typeof(CP.Objective) == LagrangeMayer
            @objective(gmp, Min, sum(Mom(Δt[i]*CP.Objective.l, gmp[:ξₜ][i]) for i in 1:nt) + Mom(CP.Objective.m, gmp[:Pₜ][end]))
        else
            Q = @variable(gmp, [nt], Meas(CP.Objective.l, support = CP.Objective.X))
            R = @variable(gmp, [nt], Meas(CP.MP, support = MP.X))
            @constraint(gmp, R[nt] + Q[nt] == Pₜ[nt])
            @objective(gmp, Max, Mom(1, Q[nt]))
        end

    elseif typeof(CP.Objective) == ExitProbability
        ∂X = ∂(CP.Objective.X)
        Qₜ = @variable(gmp, [i in 1:nt], Meas(MP.x, support = CP.Objective.X))
        Rₜ = @variable(gmp, [k in 1:length(∂X), i in 1:nt],
                            Meas([MP.x..., CP.t], support = intersect(∂X[k], @set(CP.t >= 0), @set(1 - CP.t >= 0))))
        Sₜ = @variable(gmp, [i in 1:nt], Meas([MP.x..., CP.u..., CP.t], support = intersect(CP.Objective.X, CP.U, @set(CP.t >= 0), @set(1 - CP.t >= 0))))
        @constraint(gmp, dynamics[ϕ in test_fxns, i in 1:nt],
                         louiville(MP,Qₜ[i],Sₜ[i],ϕ,CP.t,1,scale=Δt[i]) + sum(Mom(ϕ, Rₜ[k,i]) for k in length(∂X)) ==
                         (i == 1 ? split_poly(ϕ,MP.x,0)(x0...) : Mom(split_poly(ϕ, MP.x, 0), Qₜ[i-1])))
        @objective(gmp, Max, Mom(1, Qₜ[nt]))
    end
    return gmp
end

function value_function_approximation(gmp, t, trange)
    V_pieces = Dict(i => polynomial(0.0*t) for i in 1:length(trange))
    for idx in keys(gmp[:dynamics])
        i, ϕ = idx[2], idx[1]
        t0 = (i > 1 ? trange[i-1] : 0)
        V_pieces[i] += dual(gmp[:dynamics][idx])[1] * subs(ϕ, t => (t-t0)/(trange[i]-t0))
    end
    V = function (x,t)
            if t > trange[end] || t < 0
                @warn string("t = ", t, " outside the domain of V")
                i = t < 0 ? 0 : length(trange) - 1
            else
                i = findlast(ti -> ti <= t, trange[1:end-1])
            end
            return isnothing(i) ? (V_pieces[1], V_pieces[1](x...,t)) : (V_pieces[i+1], V_pieces[i+1](x...,t))
        end
    return V
end

function optimal_control(CP::ControlProcess, x0::AbstractVector, order::Int, trange::AbstractVector, solver; value_function = false)
    gmp = control_gmp(CP, x0, order, trange, solver)
    optimize!(gmp)
    lb, stat, time = extract_solution(gmp)
    if value_function
        return Bounds(:control_problem, order, lb, InfoData(stat, time)), value_function_approximation(gmp, CP.t, trange)
    else
        return Bounds(:control_problem, order, lb, InfoData(stat, time))
    end
end


text:
#using ClustForOpt_priv
using JuMP
using Clp
using DataFrames
 # parameters

acf = 0.1 # annualized cost of capital factor

c_cap_s = (11.94e6 + acf*2.58e9)*20/100 # $/GW   # 20/100 ~ 0.5$/W, 40/100~ 1$/W
c_cap_w = 30.98e6 + acf*1.74e9 # $/GW
c_cap_gas = 11.96e6 + acf*0.968e9 # $/GW
c_cap_nuc = 92.04e6 + acf*3.82e9 # $/GW
c_cap_coal = 28.15e6 + acf*2.22e9 # $/GW
c_cap = [c_cap_coal,c_cap_nuc,c_cap_gas,c_cap_w,c_cap_s]

c_var_s = 0.5 # $/GWh
c_var_w = 0.5 # $/GWh
c_var_gas = 41500 # $/GWh
c_var_nuc = 7025 # $/GWh
c_var_coal = 23000 # $/GWh
c_var = [c_var_s,c_var_w,c_var_gas,c_var_nuc,c_var_coal]
c_var = [c_var_coal,c_var_nuc,c_var_gas,c_var_w,c_var_s]

c_lost = 1000e6 # $/GW # VOLL

n_years = 10  # depreciation period for capital
r = 0.07 # real discount rate

 # load input data

de= readtable("/home/hteich/.julia/v0.6/ClustForOpt_priv/data/texas_merrick/demand.txt",separator=' ')[:DEM] # MW
d=reshape(de,(size(de)[1],1))
 # load growth
d=1.486*d
d=d/1000 # GW

 #d = rand(8760,1)  # demand
 # load growth *1.5 (footnote 6 merrick)
a_s= readtable("/home/hteich/.julia/v0.6/ClustForOpt_priv/data/texas_merrick/TexInsolationFactorV1.txt",separator=' ')[:solar_61]
a_s=reshape(a_s,(size(a_s)[1],1))
  a_s = a_s/1000
 #a_s = rand(8760,1) # availability solar
 #a_w = rand(8760,1) # availability wind
a_w= readtable("/home/hteich/.julia/v0.6/ClustForOpt_priv/data/texas_merrick/windfactor2.txt",separator=' ')[:Wind_61]
a_w=reshape(a_w,(size(a_w)[1],1))
a_gas = ones(8760,1)
a_nuc = ones(8760,1)
a_coal =ones(8760,1)
a= [a_coal,a_nuc,a_gas,a_w,a_s]

N=ones(8760,1)

del_t = 1 # hour

  ### SETS ####
K = size(d)[1]
T = size(d)[2]
G_names = ["coal","nuclear","gas","wind","solar"] # generators
n_G = length(G_names)

  ##### MODEL #####
m = Model(optimizer=ClpSolver())

  ### VARIABLES ##
  ### DESIGN VARIABLES ###
# generation capacity build
@variable(m,G_max[1:n_G] >= 0) # GWh

 ### OPERATIONAL VARIABLES ####
# generation of individual generator per time step
@variable(m,G[1:n_G,1:K,1:T] >=0)

 ##### OBJECTIVE FUNCTION #####
 @objective(m,Min,sum(c_cap[g]*G_max[g] for g=1:n_G)  +
     sum(sum(sum( c_var[g]*G[g,k,t]
           for k=1:K )
         for t=1:T )
       for g=1:n_G )
     )

 #### CONSTRAINTS ######
 # energy balance constraint
@constraint(m,[k=1:K,t=1:T], sum(G[g,k,t] for g=1:n_G) >= d[k,t])
 # maximum power production constraint
@constraint(m,[g=1:n_G,k=1:K,t=1:T], G[g,k,t] <= G_max[g] * a[g][k,t] )

  #### SOLVE ####
  status = solve(m)

 #### Output results ####
 G_max = getvalue(G_max)
 G = getvalue(G) # installed capacity
 G_gen =sum(G[:,:,1],2) # generated power

 if true
    using PyPlot

    figure()
    bt = zeros(length(G_max))
    for i=2:length(bt)
      bt[i]=cumsum(G_max)[i-1]
    end
    for i=1:length(G_max)
      bar(1, G_max[i],bottom=bt[i],color=cols[i],label=G_names[i])
    end
    ylabel("GW")
    legend()
    title("Capacity built")

    figure()
    bt = zeros(length(G_gen))
    for i=2:length(bt)
      bt[i]=cumsum(G_gen)[i-1]
    end
    for i=1:length(G_gen)
      bar(1, G_gen[i],bottom=bt[i],color=cols[i],label=G_names[i])
    end
    ylabel("GWh")
    legend()
    title("Power generated")

 end

 # TODO - add clustering, check on paper what to do next


text:
printstyled("  convert\n", color=:light_green)
C = randSeisChannel()
C1 = convert(EventChannel, C)
TD = convert(EventTraceData, C1)
@test sizeof(TD) > 136

C2 = convert(SeisChannel, C)
@test C == C2

EC = convert(EventChannel, randSeisChannel())

TD = unsafe_convert(EventTraceData, randSeisData(10))
for f in datafields
  @test length(getfield(TD, f)) == 10
end

S = unsafe_convert(SeisData, randSeisEvent(10).data)
for f in datafields
  @test length(getfield(S, f)) == 10
end

printstyled("  show\n", color=:light_green)
redirect_stdout(out) do
  show(SeisHdr())
  show(SeisSrc())
  show(SeisEvent())
  show(EventTraceData())
  show(EventChannel())
  show(randSeisHdr())
  show(randSeisEvent())
  summary(randSeisEvent())
  summary(randSeisHdr())
  show(EQMag())
  show(EQLoc())
  show(SourceTime())
end

# EQMag
Δ = 75.3
m1 = EQMag(3.2f0, "Ml", 23, Δ, "localmag")
m2 = EQMag(3.2f0, "Ml", 23, Δ, "localmag")
@test hash(m1) == hash(m2)

# SourceTime
@test isempty(SourceTime())
ST1 = SourceTime()
@test hash(ST1) == hash(SourceTime())


text:
"""
    AnimationTask

A thin wrapper around a `Task` together with an interrupt_switch that signals
the animation loop to exit.
"""
struct AnimationTask
    task::Task
    interrupt_switch::Ref{Bool}
end

macro async_showerr(ex)
    quote
        animationtask = @async try
            eval($(esc(ex)))
        catch err
            bt = catch_backtrace()
            println("Asynchronous animation errored:")
            showerror(stderr, err, bt)
        end
    end
end

"""
    animate_async(f::Function, anims::FiniteLengthAnimation...; duration::Real, fps::Int = 30)

Start an asynchronous animation where in each frame `f` is called with the current
animation time as well as the current value of each `Animation` in `anims`.

Returns an `AnimationTask` which can be stopped with `stop(animationtask)`.

Example:

    animate_async(anim1, anim2) do t, a1, a2
        # do something (e.g. with a plot or other visual object)
    end
"""
function animate_async(f::Function, anims::FiniteLengthAnimation...;
        duration = maximum(duration, anims),
        fps::Int = 30)

    frameduration = 1 / fps

    t_start = time()
    t_target = t_start

    interrupt_switch = Ref(false)

    task = @async_showerr while !interrupt_switch[]

        t_current = time()
        t_relative = t_current - t_start

        f(t_relative, (a(t_relative) for a in anims)...)

        if t_relative >= duration
            break
        end

        # always try to hit the next target exactly one frame duration away from
        # the last to avoid drift
        t_target += frameduration
        sleeptime = t_target - time()
        sleep(max(0.001, sleeptime))
    end

    AnimationTask(task, interrupt_switch)
end

"""
    stop(at::AnimationTask)

Stop a running `AnimationTask`. This only sets a flag for the animation loop to
exit, it won't kill a task that is stuck. You can manipulate the `Task` stored
in the `AnimationTask` directly if you need more control.
"""
stop(at::AnimationTask) = at.interrupt_switch[] = true

Base.wait(at::AnimationTask) = wait(at.task)


text:
function part1(input)
    credentials = collect_credentials(readlines(input))
    return count(valid_credential, credentials)
end

function part2(input)
    credentials = collect_credentials(readlines(input))
    return count(valid_credential_values, credentials)
end

function collect_credentials(lines)
    credentials = Dict{String, String}[]
    credential = Dict{String, String}()
    for line in lines
        if isempty(line)
            push!(credentials, credential)
            credential = Dict{String, String}()
        end
        for pair in split(line, " ", keepempty = false)
            key, value = split(pair, ":")
            credential[key] = value
        end
    end
    !isempty(credential) && push!(credentials, credential)
    return credentials
end

function valid_credential(credential)
    required_fields = ["byr", "iyr", "eyr", "hgt", "hcl", "ecl", "pid"]
    return isempty(setdiff(required_fields, keys(credential)))
end

function valid_credential_values(credential)
    valid_credential(credential) || return false

    byr = credential["byr"]
    is_n_digit_number(byr, 4) || return false
    1920 <= parse(Int, byr) <= 2002 || return false

    iyr = credential["iyr"]
    is_n_digit_number(iyr, 4) || return false
    2010 <= parse(Int, iyr) <= 2020 || return false

    eyr = credential["eyr"]
    is_n_digit_number(eyr, 4) || return false
    2020 <= parse(Int, eyr) <= 2030 || return false

    hgt = credential["hgt"]
    if endswith(hgt, "in")
        n = tryparse(Int, first(split(hgt, "in")))
        isnothing(n) && return false
        59 <= n <= 76 || return false
    elseif endswith(hgt, "cm")
        n = tryparse(Int, first(split(hgt, "cm")))
        isnothing(n) && return false
        150 <= n <= 193 || return false
    else
        return false
    end

    hcl = credential["hcl"]
    isnothing(match(r"#[0-9a-f][0-9a-f][0-9a-f][0-9a-f][0-9a-f][0-9a-f]", hcl)) && return false

    ecl = credential["ecl"]
    ecl in ["amb", "blu", "brn", "gry", "grn", "hzl", "oth"] || return false

    pid = credential["pid"]
    is_n_digit_number(pid, 9) || return false

    return true
end

function is_n_digit_number(s, n)
    return !isnothing(match(Regex("^" * "\\d" ^ n * "\$"), s))
end


text:
include("instructions.jl")

function manager(instruction::T, model::GeneModel) where {T <: IndividualInstruction}
    if hasfield(T, :tag)
        tag = instruction.tag
        for id in collect(keys(model.individuals))
            if tag = nothing
                instruction.name(model, instruction, id)
            end
        end
    end
end

instruction = RandomGenome(2,"binary",1.0, nothing, nothing, false, generate_genotype)

individuals = Dict(1,BasicIndividual(zeros(2,2)))

model = GeneModel(individuals)

manager(instruction, model)


text: invmod(a, b)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
#!/usr/bin/env julia
import Celeste: main
exit(main())


text:
# single board
include("simple.jl")
include("waveforms.jl")
include("slowDAC.jl")
include("slowDACMultiChannelWithOffset.jl")

# cluster
#include("cluster.jl")

max_stars_repo_name: tknopp/RedPitayaDAQServer

max_issues_repo_name: Fypsilonn/RedPitayaDAQServer

max_forks_repo_name: Fypsilonn/RedPitayaDAQServer


text:
# Sampler for von Mises-Fisher

immutable VonMisesFisherSampler
    p::Int          # the dimension
    κ::Float64
    b::Float64
    x0::Float64
    c::Float64
    Q::Matrix{Float64}
end

function VonMisesFisherSampler(μ::Vector{Float64}, κ::Float64)
    p = length(μ)
    b = _vmf_bval(p, κ)
    x0 = (1.0 - b) / (1.0 + b)
    c = κ * x0 + (p - 1) * log1p(-abs2(x0))
    Q = _vmf_rotmat(μ)
    VonMisesFisherSampler(p, κ, b, x0, c, Q)
end

function _rand!(spl::VonMisesFisherSampler, x::AbstractVector, t::AbstractVector)
    w = _vmf_genw(spl)
    p = spl.p
    t[1] = w
    s = 0.0
    for i = 2:p
        t[i] = ti = randn()
        s += abs2(ti)
    end

    # normalize t[2:p]
    r = sqrt((1.0 - abs2(w)) / s)
    for i = 2:p
        t[i] *= r
    end

    # rotate
    A_mul_B!(x, spl.Q, t)
    return x
end

_rand!(spl::VonMisesFisherSampler, x::AbstractVector) = _rand!(spl, x, Array(Float64, length(x)))

function _rand!(spl::VonMisesFisherSampler, x::AbstractMatrix)
    t = Array(Float64, size(x, 1))
    for j = 1:size(x, 2)
        _rand!(spl, view(x,:,j), t)
    end
    return x
end

### Core computation

_vmf_bval(p::Int, κ::Real) = (p - 1) / (2.0κ + sqrt(4 * abs2(κ) + abs2(p - 1)))

function _vmf_genw(p, b, x0, c, κ)
    # generate the W value -- the key step in simulating vMF
    #
    #   following movMF's document
    #

    r = (p - 1) / 2.0
    betad = Beta(r, r)
    z = rand(betad)
    w = (1.0 - (1.0 + b) * z) / (1.0 - (1.0 - b) * z)
    while κ * w + (p - 1) * log(1 - x0 * w) - c < log(rand())
        z = rand(betad)
        w = (1.0 - (1.0 + b) * z) / (1.0 - (1.0 - b) * z)
    end
    return w::Float64
end

_vmf_genw(s::VonMisesFisherSampler) = _vmf_genw(s.p, s.b, s.x0, s.c, s.κ)

function _vmf_rotmat(u::Vector{Float64})
    # construct a rotation matrix Q
    # s.t. Q * [1,0,...,0]^T --> u
    #
    # Strategy: construct a full-rank matrix
    # with first column being u, and then
    # perform QR factorization
    #

    p = length(u)
    A = zeros(p, p)
    copy!(view(A,:,1), u)

    # let k the be index of entry with max abs
    k = 1
    a = abs(u[1])
    for i = 2:p
        @inbounds ai = abs(u[i])
        if ai > a
            k = i
            a = ai
        end
    end

    # other columns of A will be filled with
    # indicator vectors, except the one
    # that activates the k-th entry
    i = 1
    for j = 2:p
        if i == k
            i += 1
        end
        A[i, j] = 1.0
    end

    # perform QR factorization
    Q = full(qrfact!(A)[:Q])
    if dot(view(Q,:,1), u) < 0.0  # the first column was negated
        for i = 1:p
            @inbounds Q[i,1] = -Q[i,1]
        end
    end
    return Q
end