text:
#  Copyright 2016, Iain Dunning, Joey Huchette, Miles Lubin, and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at http://mozilla.org/MPL/2.0/.
#############################################################################
# JuMP
# An algebraic modelling langauge for Julia
# See http://github.com/JuliaOpt/JuMP.jl
#############################################################################
# test/operator.jl
# Testing operator overloading is correct
#############################################################################
using JuMP, FactCheck

# To ensure the tests work on both Windows and Linux/OSX, we need
# to use the correct comparison operators in the strings
const  leq = JuMP.repl[:leq]
const  geq = JuMP.repl[:geq]
const   eq = JuMP.repl[:eq]
const Vert = JuMP.repl[:Vert]
const sub2 = JuMP.repl[:sub2]

facts("[operator] Testing basic operator overloads") do
    m = Model()
    @variable(m, w)
    @variable(m, x)
    @variable(m, y)
    @variable(m, z)
    aff = 7.1 * x + 2.5
    @fact string(aff) --> "7.1 x + 2.5"
    aff2 = 1.2 * y + 1.2
    @fact string(aff2) --> "1.2 y + 1.2"
    q = 2.5 * y * z + aff
    @fact string(q) --> "2.5 y*z + 7.1 x + 2.5"
    q2 = 8 * x * z + aff2
    @fact string(q2) --> "8 x*z + 1.2 y + 1.2"
    q3 = 2 * x * x + 1 * y * y + z + 3
    @fact string(q3) --> "2 x² + y² + z + 3"

    nrm = norm([w,1-w])
    @fact string(nrm) --> "$Vert[w,-w + 1]$Vert$sub2"
    socexpr = 1.5*nrm - 2 - w
    @fact string(socexpr) --> "1.5 $Vert[w,-w + 1]$Vert$sub2 - w - 2"

    @fact isequal(3w + 2y, 3w +2y) --> true
    @fact isequal(3w + 2y + 1, 3w + 2y) --> false

    # Different objects that must all interact:
    # 1. Number
    # 2. Variable
    # 3. AffExpr
    # 4. QuadExpr

    # 1. Number tests
    context("Number--???") do
    # 1-1 Number--Number - nope!
    # 1-2 Number--Variable
    @fact string(4.13 + w) --> "w + 4.13"
    @fact string(3.16 - w) --> "-w + 3.16"
    @fact string(5.23 * w) --> "5.23 w"
    @fact_throws  2.94 / w
    @fact string(@LinearConstraint(2.1 ≤ w)) --> "-w $leq -2.1"
    @fact string(@LinearConstraint(2.1 == w)) --> "-w $eq -2.1"
    @fact string(@LinearConstraint(2.1 ≥ w)) --> "-w $geq -2.1"
    # 1-3 Number--Norm
    @fact string(4.13 + nrm) --> "$Vert[w,-w + 1]$Vert$sub2 + 4.13"
    @fact string(3.16 - nrm) --> "-1.0 $Vert[w,-w + 1]$Vert$sub2 + 3.16"
    @fact string(5.23 * nrm) --> "5.23 $Vert[w,-w + 1]$Vert$sub2"
    @fact_throws 2.94 / nrm
    @fact_throws @SOCConstraint(2.1 ≤ nrm)
    @fact_throws @SOCConstraint(2.1 == nrm)
    @fact string(@SOCConstraint(2.1 ≥ nrm)) --> "$Vert[w,-w + 1]$Vert$sub2 $leq 2.1"
    # 1-4 Number--AffExpr
    @fact string(1.5 + aff) --> "7.1 x + 4"
    @fact string(1.5 - aff) --> "-7.1 x - 1"
    @fact string(2 * aff) --> "14.2 x + 5"
    @fact_throws  2 / aff
    @fact string(@LinearConstraint(1 ≤ aff)) --> "-7.1 x $leq 1.5"
    @fact string(@LinearConstraint(1 == aff)) --> "-7.1 x $eq 1.5"
    @fact string(@LinearConstraint(1 ≥ aff)) --> "-7.1 x $geq 1.5"
    # 1-5 Number--QuadExpr
    @fact string(1.5 + q) --> "2.5 y*z + 7.1 x + 4"
    @fact string(1.5 - q) --> "-2.5 y*z - 7.1 x - 1"
    @fact string(2 * q) --> "5 y*z + 14.2 x + 5"
    @fact_throws  2 / q
    @fact string(@QuadConstraint(1 ≤ q)) --> "-2.5 y*z - 7.1 x - 1.5 $leq 0"
    @fact string(@QuadConstraint(1 == q)) --> "-2.5 y*z - 7.1 x - 1.5 $eq 0"
    @fact string(@QuadConstraint(1 ≥ q)) --> "-2.5 y*z - 7.1 x - 1.5 $geq 0"
    # 1-6 Number--SOCExpr
    @fact string(1.5 + socexpr) --> "1.5 $Vert[w,-w + 1]$Vert$sub2 - w - 0.5"
    @fact string(1.5 - socexpr) --> "-1.5 $Vert[w,-w + 1]$Vert$sub2 + w + 3.5"
    @fact string(1.5 * socexpr) --> "2.25 $Vert[w,-w + 1]$Vert$sub2 - 1.5 w - 3"
    @fact_throws 1.5 / socexpr
    @fact_throws @SOCConstraint(2.1 ≤ socexpr)
    @fact_throws @SOCConstraint(2.1 == socexpr)
    @fact string(@SOCConstraint(2.1 ≥ socexpr)) --> "1.5 $Vert[w,-w + 1]$Vert$sub2 $leq w + 4.1"
    end

    # 2. Variable tests
    context("Variable--???") do
    # 2-0 Variable unary
    @fact (+x) --> exactly(x)
    @fact string(-x) --> "-x"
    # 2-1 Variable--Number
    @fact string(w + 4.13) --> "w + 4.13"
    @fact string(w - 4.13) --> "w - 4.13"
    @fact string(w * 4.13) --> "4.13 w"
    @fact string(w / 2.00) --> "0.5 w"
    @fact w == w --> true
    @fact string(@LinearConstraint(w ≤ 1)) --> "w $leq 1"
    @fact string(@LinearConstraint(w == 1)) --> "w $eq 1"
    @fact string(@LinearConstraint(w ≥ 1)) --> "w $geq 1"
    @fact string(@QuadConstraint(x*y ≤ 1)) --> "x*y - 1 $leq 0"
    @fact string(@QuadConstraint(x*y == 1)) --> "x*y - 1 $eq 0"
    @fact string(@QuadConstraint(x*y ≥ 1)) --> "x*y - 1 $geq 0"
    # 2-2 Variable--Variable
    @fact string(w + x) --> "w + x"
    @fact string(w - x) --> "w - x"
    @fact string(w * x) --> "w*x"
    @fact string(x - x) --> "0"
    @fact_throws  w / x
    @fact string(@LinearConstraint(w ≤ x)) --> "w - x $leq 0"
    @fact string(@LinearConstraint(w == x)) --> "w - x $eq 0"
    @fact string(@LinearConstraint(w ≥ x)) --> "w - x $geq 0"
    @fact string(@QuadConstraint(y*z ≤ x)) --> "y*z - x $leq 0"
    @fact string(@QuadConstraint(y*z == x)) --> "y*z - x $eq 0"
    @fact string(@QuadConstraint(y*z ≥ x)) --> "y*z - x $geq 0"
    @fact string(@LinearConstraint(x ≤ x)) --> "0 $leq 0"
    @fact string(@LinearConstraint(x == x)) --> "0 $eq 0"
    @fact string(@LinearConstraint(x ≥ x)) --> "0 $geq 0"
    # 2-3 Variable--Norm
    @fact string(w + nrm) --> "$Vert[w,-w + 1]$Vert$sub2 + w"
    @fact string(w - nrm) --> "-1.0 $Vert[w,-w + 1]$Vert$sub2 + w"
    @fact_throws w * nrm
    @fact_throws w / nrm
    @fact_throws @SOCConstraint(w ≤ nrm)
    @fact_throws @SOCConstraint(w == nrm)
    @fact string(@SOCConstraint(w ≥ nrm)) --> "$Vert[w,-w + 1]$Vert$sub2 $leq w"
    # 2-4 Variable--AffExpr
    @fact string(z + aff) --> "7.1 x + z + 2.5"
    @fact string(z - aff) --> "-7.1 x + z - 2.5"
    @fact string(z * aff) --> "7.1 x*z + 2.5 z"
    @fact_throws  z / aff
    @fact_throws  z ≤ aff
    @fact string(@LinearConstraint(z ≤ aff)) --> "z - 7.1 x $leq 2.5"
    @fact string(@LinearConstraint(z == aff)) --> "z - 7.1 x $eq 2.5"
    @fact string(@LinearConstraint(z ≥ aff)) --> "z - 7.1 x $geq 2.5"
    @fact string(@LinearConstraint(7.1 * x - aff ≤ 0)) --> "0 $leq 2.5"
    @fact string(@LinearConstraint(7.1 * x - aff == 0)) --> "0 $eq 2.5"
    @fact string(@LinearConstraint(7.1 * x - aff ≥ 0)) --> "0 $geq 2.5"
    # 2-5 Variable--QuadExpr
    @fact string(w + q) --> "2.5 y*z + 7.1 x + w + 2.5"
    @fact string(w - q) --> "-2.5 y*z - 7.1 x + w - 2.5"
    @fact_throws  w*q
    @fact_throws  w/q
    @fact string(@QuadConstraint(w ≤ q)) --> "-2.5 y*z + w - 7.1 x - 2.5 $leq 0"
    @fact string(@QuadConstraint(w == q)) --> "-2.5 y*z + w - 7.1 x - 2.5 $eq 0"
    @fact string(@QuadConstraint(w ≥ q)) --> "-2.5 y*z + w - 7.1 x - 2.5 $geq 0"
    # 2-6 Variable--SOCExpr
    @fact string(y + socexpr) --> "1.5 $Vert[w,-w + 1]$Vert$sub2 - w + y - 2"
    @fact string(y - socexpr) --> "-1.5 $Vert[w,-w + 1]$Vert$sub2 + w + y + 2"
    @fact_throws y * socexpr
    @fact_throws y / socexpr
    @fact_throws @SOCConstraint(y ≤ socexpr)
    @fact_throws @SOCConstraint(y == socexpr)
    @fact string(@SOCConstraint(y ≥ socexpr)) --> "1.5 $Vert[w,-w + 1]$Vert$sub2 $leq y + w + 2"
    end

    # 3. Norm tests
    context("Norm--???") do
    # 3-0 Norm unary
    @fact string(+nrm) --> "$Vert[w,-w + 1]$Vert$sub2"
    @fact string(-nrm) --> "-1.0 $Vert[w,-w + 1]$Vert$sub2"
    # 3-1 Norm--Number
    @fact string(nrm + 1.5) --> "$Vert[w,-w + 1]$Vert$sub2 + 1.5"
    @fact string(nrm - 1.5) --> "$Vert[w,-w + 1]$Vert$sub2 - 1.5"
    @fact string(nrm * 1.5) --> "1.5 $Vert[w,-w + 1]$Vert$sub2"
    @fact string(nrm / 1.5) --> "0.6666666666666666 $Vert[w,-w + 1]$Vert$sub2"
    @fact string(@SOCConstraint(nrm ≤ 1.5)) --> "$Vert[w,-w + 1]$Vert$sub2 $leq 1.5"
    @fact_throws @SOCConstraint(nrm == 1.5)
    @fact_throws @SOCConstraint(nrm ≥ 1.5)
    # 3-2 Norm--Variable
    @fact string(nrm + w) --> "$Vert[w,-w + 1]$Vert$sub2 + w"
    @fact string(nrm - w) --> "$Vert[w,-w + 1]$Vert$sub2 - w"
    @fact_throws nrm * w
    @fact_throws nrm / w
    @fact string(@SOCConstraint(nrm ≤ w)) --> "$Vert[w,-w + 1]$Vert$sub2 $leq w"
    @fact_throws @SOCConstraint(nrm == w)
    @fact_throws @SOCConstraint(nrm ≥ w)
    # 3-3 Norm--Norm
    @fact_throws nrm + nrm
    @fact_throws nrm - nrm
    @fact_throws nrm * nrm
    @fact_throws nrm / nrm
    @fact_throws @SOCConstraint(nrm ≤ nrm)
    @fact_throws @SOCConstraint(nrm == nrm)
    @fact_throws @SOCConstraint(nrm ≥ nrm)
    # 3-4 Norm--AffExpr
    @fact string(nrm + aff) --> "$Vert[w,-w + 1]$Vert$sub2 + 7.1 x + 2.5"
    @fact string(nrm - aff) --> "$Vert[w,-w + 1]$Vert$sub2 - 7.1 x - 2.5"
    @fact_throws nrm * aff
    @fact_throws nrm / aff
    @fact string(@SOCConstraint(nrm ≤ aff)) --> "$Vert[w,-w + 1]$Vert$sub2 $leq 7.1 x + 2.5"
    @fact_throws @SOCConstraint(nrm == aff)
    @fact_throws @SOCConstraint(nrm ≥ aff)
    # 3-5 Norm--QuadExpr
    @fact_throws nrm + q
    @fact_throws nrm - q
    @fact_throws nrm * q
    @fact_throws nrm / q
    @fact_throws @SOCConstraint(nrm ≤ q)
    @fact_throws @SOCConstraint(nrm == q)
    @fact_throws @SOCConstraint(nrm ≥ q)
    # 3-6 Norm--SOCExpr
    @fact_throws nrm + socexpr
    @fact_throws nrm - socexpr
    @fact_throws nrm * socexpr
    @fact_throws nrm / socexpr
    @fact_throws @SOCConstraint(nrm ≤ socexpr)
    @fact_throws @SOCConstraint(nrm == socexpr)
    @fact_throws @SOCConstraint(nrm ≥ socexpr)
    end

    # 4. AffExpr tests
    context("AffExpr--???") do
    # 4-0 AffExpr unary
    @fact string(+aff) --> "7.1 x + 2.5"
    @fact string(-aff) --> "-7.1 x - 2.5"
    # 4-1 AffExpr--Number
    @fact string(aff + 1.5) --> "7.1 x + 4"
    @fact string(aff - 1.5) --> "7.1 x + 1"
    @fact string(aff * 2) --> "14.2 x + 5"
    @fact string(aff / 2) --> "3.55 x + 1.25"
    @fact_throws aff ≤ 1
    @fact aff == aff --> true
    @fact_throws aff ≥ 1
    @fact string(@LinearConstraint(aff ≤ 1)) --> "7.1 x $leq -1.5"
    @fact string(@LinearConstraint(aff == 1)) --> "7.1 x $eq -1.5"
    @fact string(@LinearConstraint(aff ≥ 1)) --> "7.1 x $geq -1.5"
    # 4-2 AffExpr--Variable
    @fact string(aff + z) --> "7.1 x + z + 2.5"
    @fact string(aff - z) --> "7.1 x - z + 2.5"
    @fact string(aff * z) --> "7.1 x*z + 2.5 z"
    @fact_throws  aff/z
    @fact string(@LinearConstraint(aff ≤ z)) --> "7.1 x - z $leq -2.5"
    @fact string(@LinearConstraint(aff == z)) --> "7.1 x - z $eq -2.5"
    @fact string(@LinearConstraint(aff ≥ z)) --> "7.1 x - z $geq -2.5"
    @fact string(@LinearConstraint(aff - 7.1 * x ≤ 0)) --> "0 $leq -2.5"
    @fact string(@LinearConstraint(aff - 7.1 * x == 0)) --> "0 $eq -2.5"
    @fact string(@LinearConstraint(aff - 7.1 * x ≥ 0)) --> "0 $geq -2.5"
    # 4-3 AffExpr--Norm
    @fact string(aff + nrm) --> "$Vert[w,-w + 1]$Vert$sub2 + 7.1 x + 2.5"
    @fact string(aff - nrm) --> "-1.0 $Vert[w,-w + 1]$Vert$sub2 + 7.1 x + 2.5"
    @fact_throws aff * nrm
    @fact_throws aff / nrm
    @fact_throws @SOCConstraint(aff ≤ nrm)
    @fact_throws @SOCConstraint(aff == nrm)
    @fact string(@SOCConstraint(aff ≥ nrm)) --> "$Vert[w,-w + 1]$Vert$sub2 $leq 7.1 x + 2.5"
    # 4-4 AffExpr--AffExpr
    @fact string(aff + aff2) --> "7.1 x + 1.2 y + 3.7"
    @fact string(aff - aff2) --> "7.1 x - 1.2 y + 1.3"
    @fact string(aff * aff2) --> "8.52 x*y + 3 y + 8.52 x + 3"
    @fact string((x+x)*(x+3)) --> string((x+3)*(x+x))  # Issue #288
    @fact_throws  aff/aff2
    @fact string(@LinearConstraint(aff ≤ aff2)) --> "7.1 x - 1.2 y $leq -1.3"
    @fact string(@LinearConstraint(aff == aff2)) --> "7.1 x - 1.2 y $eq -1.3"
    @fact string(@LinearConstraint(aff ≥ aff2)) --> "7.1 x - 1.2 y $geq -1.3"
    @fact string(@LinearConstraint(aff-aff ≤ 0)) --> "0 $leq 0"
    @fact string(@LinearConstraint(aff-aff == 0)) --> "0 $eq 0"
    @fact string(@LinearConstraint(aff-aff ≥ 0)) --> "0 $geq 0"
    # 4-5 AffExpr--QuadExpr
    @fact string(aff2 + q) --> "2.5 y*z + 1.2 y + 7.1 x + 3.7"
    @fact string(aff2 - q) --> "-2.5 y*z + 1.2 y - 7.1 x - 1.3"
    @fact_throws  aff2 * q
    @fact_throws  aff2 / q
    @fact string(@QuadConstraint(aff2 ≤ q)) --> "-2.5 y*z + 1.2 y - 7.1 x - 1.3 $leq 0"
    @fact string(@QuadConstraint(aff2 == q)) --> "-2.5 y*z + 1.2 y - 7.1 x - 1.3 $eq 0"
    @fact string(@QuadConstraint(aff2 ≥ q)) --> "-2.5 y*z + 1.2 y - 7.1 x - 1.3 $geq 0"
    # 4-6 AffExpr--SOCExpr
    @fact string(aff + socexpr) --> "1.5 $Vert[w,-w + 1]$Vert$sub2 + 7.1 x - w + 0.5"
    @fact string(aff - socexpr) --> "-1.5 $Vert[w,-w + 1]$Vert$sub2 + 7.1 x + w + 4.5"
    @fact_throws aff * socexpr
    @fact_throws aff / socexpr
    @fact_throws @SOCConstraint(aff ≤ socexpr)
    @fact_throws @SOCConstraint(aff == socexpr)
    @fact string(@SOCConstraint(aff ≥ socexpr)) --> "1.5 $Vert[w,-w + 1]$Vert$sub2 $leq 7.1 x + w + 4.5"
    end

    # 5. QuadExpr
    context("QuadExpr--???") do
    # 5-0 QuadExpr unary
    @fact string(+q) --> "2.5 y*z + 7.1 x + 2.5"
    @fact string(-q) --> "-2.5 y*z - 7.1 x - 2.5"
    # 5-1 QuadExpr--Number
    @fact string(q + 1.5) --> "2.5 y*z + 7.1 x + 4"
    @fact string(q - 1.5) --> "2.5 y*z + 7.1 x + 1"
    @fact string(q * 2) --> "5 y*z + 14.2 x + 5"
    @fact string(q / 2) --> "1.25 y*z + 3.55 x + 1.25"
    @fact q == q --> true
    @fact string(@QuadConstraint(aff2 ≤ q)) --> "-2.5 y*z + 1.2 y - 7.1 x - 1.3 $leq 0"
    @fact string(@QuadConstraint(aff2 == q)) --> "-2.5 y*z + 1.2 y - 7.1 x - 1.3 $eq 0"
    @fact string(@QuadConstraint(aff2 ≥ q)) --> "-2.5 y*z + 1.2 y - 7.1 x - 1.3 $geq 0"
    # 5-2 QuadExpr--Variable
    @fact string(q + w) --> "2.5 y*z + 7.1 x + w + 2.5"
    @fact string(q - w) --> "2.5 y*z + 7.1 x - w + 2.5"
    @fact_throws q*w
    @fact_throws q/w
    @fact string(@QuadConstraint(q ≤ w)) --> "2.5 y*z + 7.1 x - w + 2.5 $leq 0"
    @fact string(@QuadConstraint(q == w)) --> "2.5 y*z + 7.1 x - w + 2.5 $eq 0"
    @fact string(@QuadConstraint(q ≥ w)) --> "2.5 y*z + 7.1 x - w + 2.5 $geq 0"
    # 5-3 QuadExpr--Norm
    @fact_throws q + nrm
    @fact_throws q - nrm
    @fact_throws q * nrm
    @fact_throws q / nrm
    @fact_throws @SOCConstraint(q ≤ nrm)
    @fact_throws @SOCConstraint(q == nrm)
    @fact_throws @SOCConstraint(q ≥ nrm)
    # 5-4 QuadExpr--AffExpr
    @fact string(q + aff2) --> "2.5 y*z + 7.1 x + 1.2 y + 3.7"
    @fact string(q - aff2) --> "2.5 y*z + 7.1 x - 1.2 y + 1.3"
    @fact_throws  q * aff2
    @fact_throws  q / aff2
    @fact string(@QuadConstraint(q ≤ aff2)) --> "2.5 y*z + 7.1 x - 1.2 y + 1.3 $leq 0"
    @fact string(@QuadConstraint(q == aff2)) --> "2.5 y*z + 7.1 x - 1.2 y + 1.3 $eq 0"
    @fact string(@QuadConstraint(q ≥ aff2)) --> "2.5 y*z + 7.1 x - 1.2 y + 1.3 $geq 0"
    # 5-5 QuadExpr--QuadExpr
    @fact string(q + q2) --> "8 x*z + 2.5 y*z + 7.1 x + 1.2 y + 3.7"
    @fact string(q - q2) --> "-8 x*z + 2.5 y*z + 7.1 x - 1.2 y + 1.3"
    @fact_throws  q * q2
    @fact_throws  q / q2
    @fact string(@QuadConstraint(q ≤ q2)) --> "-8 x*z + 2.5 y*z + 7.1 x - 1.2 y + 1.3 $leq 0"
    @fact string(@QuadConstraint(q == q2)) --> "-8 x*z + 2.5 y*z + 7.1 x - 1.2 y + 1.3 $eq 0"
    @fact string(@QuadConstraint(q ≥ q2)) --> "-8 x*z + 2.5 y*z + 7.1 x - 1.2 y + 1.3 $geq 0"
    # 4-6 QuadExpr--SOCExpr
    @fact_throws q + socexpr
    @fact_throws q - socexpr
    @fact_throws q * socexpr
    @fact_throws q / socexpr
    @fact_throws @SOCConstraint(q ≤ socexpr)
    @fact_throws @SOCConstraint(q == socexpr)
    @fact_throws @SOCConstraint(q ≥ socexpr)
    end

    # 6. SOCExpr tests
    context("SOCExpr--???") do
    # 6-0 SOCExpr unary
    @fact string(+socexpr) --> "1.5 $Vert[w,-w + 1]$Vert$sub2 - w - 2"
    @fact string(-socexpr) --> "-1.5 $Vert[w,-w + 1]$Vert$sub2 + w + 2"
    # 6-1 SOCExpr--Number
    @fact string(socexpr + 1.5) --> "1.5 $Vert[w,-w + 1]$Vert$sub2 - w - 0.5"
    @fact string(socexpr - 1.5) --> "1.5 $Vert[w,-w + 1]$Vert$sub2 - w - 3.5"
    @fact string(socexpr * 1.5) --> "2.25 $Vert[w,-w + 1]$Vert$sub2 - 1.5 w - 3"
    @fact string(socexpr / 1.5) --> "$Vert[w,-w + 1]$Vert$sub2 - 0.6666666666666666 w - 1.3333333333333333"
    @fact string(@SOCConstraint(socexpr ≤ 1.5)) --> "1.5 $Vert[w,-w + 1]$Vert$sub2 $leq w + 3.5"
    @fact_throws @SOCConstraint(socexpr == 1.5)
    @fact_throws @SOCConstraint(socexpr ≥ 1.5)
    # 6-2 SOCExpr--Variable
    @fact string(socexpr + y) --> "1.5 $Vert[w,-w + 1]$Vert$sub2 - w + y - 2"
    @fact string(socexpr - y) --> "1.5 $Vert[w,-w + 1]$Vert$sub2 - w - y - 2"
    @fact_throws socexpr * y
    @fact_throws socexpr / y
    @fact string(@SOCConstraint(socexpr ≤ y)) --> "1.5 $Vert[w,-w + 1]$Vert$sub2 $leq w + y + 2"
    @fact_throws @SOCConstraint(socexpr == y)
    @fact_throws @SOCConstraint(socexpr ≥ y)
    # 6-3 SOCExpr--Norm
    @fact_throws socexpr + nrm
    @fact_throws socexpr - nrm
    @fact_throws socexpr * nrm
    @fact_throws socexpr / nrm
    @fact_throws @SOCConstraint(socexpr ≤ nrm)
    @fact_throws @SOCConstraint(socexpr == nrm)
    @fact_throws @SOCConstraint(socexpr ≥ nrm)
    # 6-4 SOCExpr--AffExpr
    @fact string(socexpr + aff) --> "1.5 $Vert[w,-w + 1]$Vert$sub2 - w + 7.1 x + 0.5"
    @fact string(socexpr - aff) --> "1.5 $Vert[w,-w + 1]$Vert$sub2 - w - 7.1 x - 4.5"
    @fact_throws socexpr * aff
    @fact_throws socexpr / aff
    @fact string(@SOCConstraint(socexpr ≤ aff)) --> "1.5 $Vert[w,-w + 1]$Vert$sub2 $leq w + 7.1 x + 4.5"
    @fact_throws @SOCConstraint(socexpr == aff)
    @fact_throws @SOCConstraint(socexpr ≥ aff)
    # 6-5 SOCExpr--QuadExpr
    @fact_throws socexpr + q
    @fact_throws socexpr - q
    @fact_throws socexpr * q
    @fact_throws socexpr / q
    @fact_throws @SOCConstraint(socexpr ≤ q)
    @fact_throws @SOCConstraint(socexpr == q)
    @fact_throws @SOCConstraint(socexpr ≥ q)
    # 6-6 SOCExpr--SOCExpr
    @fact_throws socexpr + socexpr
    @fact_throws socexpr - socexpr
    @fact_throws socexpr * socexpr
    @fact_throws socexpr / socexpr
    @fact_throws @SOCConstraint(socexpr ≤ socexpr)
    @fact_throws @SOCConstraint(socexpr == socexpr)
    @fact_throws @SOCConstraint(socexpr ≥ socexpr)
    end
end

facts("[operator] Higher-level operators") do
context("sum") do
    sum_m = Model()
    @variable(sum_m, 0 ≤ matrix[1:3,1:3] ≤ 1, start = 1)
    # sum(j::JuMPArray{Variable})
    @fact string(sum(matrix)) --> "matrix[1,1] + matrix[2,1] + matrix[3,1] + matrix[1,2] + matrix[2,2] + matrix[3,2] + matrix[1,3] + matrix[2,3] + matrix[3,3]"
    # sum(j::JuMPArray{Variable}) in a macro
    @objective(sum_m, Max, sum(matrix))
    @fact string(sum_m.obj) --> "matrix[1,1] + matrix[2,1] + matrix[3,1] + matrix[1,2] + matrix[2,2] + matrix[3,2] + matrix[1,3] + matrix[2,3] + matrix[3,3]"

    # sum{T<:Real}(j::JuMPArray{T})
    @fact sum(getvalue(matrix)) --> roughly(9, 1e-6)
    # sum(j::Array{Variable})
    @fact string(sum(matrix[1:3,1:3])) --> string(sum(matrix))
    # sum(affs::Array{AffExpr})
    @fact string(sum([2*matrix[i,j] for i in 1:3, j in 1:3])) --> "2 matrix[1,1] + 2 matrix[2,1] + 2 matrix[3,1] + 2 matrix[1,2] + 2 matrix[2,2] + 2 matrix[3,2] + 2 matrix[1,3] + 2 matrix[2,3] + 2 matrix[3,3]"

    S = [1,3]
    @variable(sum_m, x[S], start=1)
    # sum(j::JuMPDict{Variable})
    @fact length(string(sum(x))) --> 11 # order depends on hashing
    @fact contains(string(sum(x)),"x[1]") --> true
    @fact contains(string(sum(x)),"x[3]") --> true
    # sum{T<:Real}(j::JuMPDict{T})
    @fact sum(getvalue(x)) --> 2
end

context("dot") do
    dot_m = Model()
    @variable(dot_m, 0 ≤ x[1:3] ≤ 1)
    c = vcat(1:3)
    @fact string(dot(c,x)) --> "x[1] + 2 x[2] + 3 x[3]"
    @fact string(dot(x,c)) --> "x[1] + 2 x[2] + 3 x[3]"

    A = [1 3 ; 2 4]
    @variable(dot_m, 1 ≤ y[1:2,1:2] ≤ 1)
    @fact string(vecdot(A,y)) --> "y[1,1] + 2 y[2,1] + 3 y[1,2] + 4 y[2,2]"
    @fact string(vecdot(y,A)) --> "y[1,1] + 2 y[2,1] + 3 y[1,2] + 4 y[2,2]"

    B = ones(2,2,2)
    @variable(dot_m, 0 ≤ z[1:2,1:2,1:2] ≤ 1)
    @fact string(vecdot(B,z)) --> "z[1,1,1] + z[2,1,1] + z[1,2,1] + z[2,2,1] + z[1,1,2] + z[2,1,2] + z[1,2,2] + z[2,2,2]"
    @fact string(vecdot(z,B)) --> "z[1,1,1] + z[2,1,1] + z[1,2,1] + z[2,2,1] + z[1,1,2] + z[2,1,2] + z[1,2,2] + z[2,2,2]"

    @objective(dot_m, Max, dot(x, ones(3)) - vecdot(y, ones(2,2)) )
    #solve(dot_m)
    for i in 1:3
        setvalue(x[i], 1)
    end
    for i in 1:2, j in 1:2
        setvalue(y[i,j], 1)
    end
    @fact dot(c, getvalue(x)) --> roughly( 6, 1e-6)
    @fact vecdot(A, getvalue(y)) --> roughly(10, 1e-6)

    # https://github.com/JuliaOpt/JuMP.jl/issues/656
    issue656 = Model()
    @variable(issue656, x)
    floats = Float64[i for i in 1:2]
    anys   = Array(Any, 2)
    anys[1] = 10
    anys[2] = 20 + x
    @fact dot(floats, anys) --> 10 + 40 + 2x
end
end

module TestHelper # weird scoping behavior with FactCheck...
    using JuMP
    _lt(x,y) = (x.col < y.col)
    function sort_expr!(x::AffExpr)
        idx = sortperm(x.vars, lt=_lt)
        x.vars = x.vars[idx]
        x.coeffs = x.coeffs[idx]
        return x
    end

    vec_eq(x,y) = vec_eq([x;], [y;])

    function vec_eq(x::AbstractArray, y::AbstractArray)
        size(x) == size(y) || return false
        for i in 1:length(x)
            v, w = convert(AffExpr,x[i]), convert(AffExpr,y[i])
            sort_expr!(v)
            sort_expr!(w)
            string(v) == string(w) || return false
        end
        return true
    end

    function vec_eq(x::Array{QuadExpr}, y::Array{QuadExpr})
        size(x) == size(y) || return false
        for i in 1:length(x)
            string(x[i]) == string(y[i]) || return false
        end
        return true
    end
end

facts("Vectorized operations") do

context("Transpose") do
    m = Model()
    @variable(m, x[1:3])
    @variable(m, y[1:2,1:3])
    @variable(m, z[2:5])
    @fact TestHelper.vec_eq(x', [x[1] x[2] x[3]]) --> true
    @fact TestHelper.vec_eq(transpose(x), [x[1] x[2] x[3]]) --> true
    @fact TestHelper.vec_eq(y', [y[1,1] y[2,1]
                      y[1,2] y[2,2]
                      y[1,3] y[2,3]]) --> true
    @fact TestHelper.vec_eq(transpose(y),
                [y[1,1] y[2,1]
                 y[1,2] y[2,2]
                 y[1,3] y[2,3]]) --> true
    @fact_throws z'
    @fact_throws transpose(z)
end

context("Vectorized arithmetic") do
    m = Model()
    @variable(m, x[1:3])
    A = [2 1 0
         1 2 1
         0 1 2]
    B = sparse(A)
    @fact TestHelper.vec_eq(A*x, [2x[1] +  x[2]
                        2x[2] + x[1] + x[3]
                        x[2] + 2x[3]]) --> true
    @fact TestHelper.vec_eq(A*x, B*x) --> true
    @fact TestHelper.vec_eq(A*x, @JuMP.Expression(B*x)) --> true
    @fact TestHelper.vec_eq(@JuMP.Expression(A*x), @JuMP.Expression(B*x)) --> true
    @fact TestHelper.vec_eq(x'*A, [2x[1]+x[2]; 2x[2]+x[1]+x[3]; x[2]+2x[3]]') --> true
    @fact TestHelper.vec_eq(x'*A, x'*B) --> true
    @fact TestHelper.vec_eq(x'*A, @JuMP.Expression(x'*B)) --> true
    @fact TestHelper.vec_eq(@JuMP.Expression(x'*A), @JuMP.Expression(x'*B)) --> true
    @fact TestHelper.vec_eq(x'*A*x, [2x[1]*x[1] + 2x[1]*x[2] + 2x[2]*x[2] + 2x[2]*x[3] + 2x[3]*x[3]]) --> true
    @fact TestHelper.vec_eq(x'A*x, x'*B*x) --> true
    @fact TestHelper.vec_eq(x'*A*x, @JuMP.Expression(x'*B*x)) --> true
    @fact TestHelper.vec_eq(@JuMP.Expression(x'*A*x), @JuMP.Expression(x'*B*x)) --> true

    y = A*x
    @fact TestHelper.vec_eq(-x, [-x[1], -x[2], -x[3]]) --> true
    @fact TestHelper.vec_eq(-y, [-2x[1] -  x[2]
                       -x[1] - 2x[2] -  x[3]
                               -x[2] - 2x[3]]) --> true
    @fact TestHelper.vec_eq(y + 1, [2x[1] +  x[2]         + 1
                          x[1] + 2x[2] +  x[3] + 1
                                  x[2] + 2x[3] + 1]) --> true
    @fact TestHelper.vec_eq(y - 1, [2x[1] +  x[2]         - 1
                          x[1] + 2x[2] +  x[3] - 1
                                  x[2] + 2x[3] - 1]) --> true
    @fact TestHelper.vec_eq(y + 2ones(3), [2x[1] +  x[2]         + 2
                                 x[1] + 2x[2] +  x[3] + 2
                                         x[2] + 2x[3] + 2]) --> true
    @fact TestHelper.vec_eq(y - 2ones(3), [2x[1] +  x[2]         - 2
                                 x[1] + 2x[2] +  x[3] - 2
                                         x[2] + 2x[3] - 2]) --> true
    @fact TestHelper.vec_eq(2ones(3) + y, [2x[1] +  x[2]         + 2
                                 x[1] + 2x[2] +  x[3] + 2
                                         x[2] + 2x[3] + 2]) --> true
    @fact TestHelper.vec_eq(2ones(3) - y, [-2x[1] -  x[2]         + 2
                                 -x[1] - 2x[2] -  x[3] + 2
                                         -x[2] - 2x[3] + 2]) --> true
    @fact TestHelper.vec_eq(y + x, [3x[1] +  x[2]
                          x[1] + 3x[2] +  x[3]
                                  x[2] + 3x[3]]) --> true
    @fact TestHelper.vec_eq(x + y, [3x[1] +  x[2]
                          x[1] + 3x[2] +  x[3]
                                  x[2] + 3x[3]]) --> true
    @fact TestHelper.vec_eq(2y + 2x, [6x[1] + 2x[2]
                           2x[1] + 6x[2] + 2x[3]
                                   2x[2] + 6x[3]]) --> true
    @fact TestHelper.vec_eq(y - x, [ x[1] + x[2]
                          x[1] + x[2] + x[3]
                                 x[2] + x[3]]) --> true
    @fact TestHelper.vec_eq(x - y, [-x[1] - x[2]
                         -x[1] - x[2] - x[3]
                                -x[2] - x[3]]) --> true
    @fact TestHelper.vec_eq(y + x[:], [3x[1] +  x[2]
                             x[1] + 3x[2] +  x[3]
                                     x[2] + 3x[3]]) --> true
    @fact TestHelper.vec_eq(x[:] + y, [3x[1] +  x[2]
                             x[1] + 3x[2] +  x[3]
                                     x[2] + 3x[3]]) --> true

    @fact TestHelper.vec_eq(@JuMP.Expression(A*x/2), A*x/2) --> true
end

context("Dot-ops") do
    m = Model()
    @variable(m, x[1:2,1:2])
    A = [1 2;
         3 4]
    B = sparse(A)
    y = SparseMatrixCSC(2, 2, copy(B.colptr), copy(B.rowval), vec(x))
    @fact TestHelper.vec_eq(A.+x, [1+x[1,1]  2+x[1,2];
                                   3+x[2,1]  4+x[2,2]]) --> true
    @fact TestHelper.vec_eq(A.+x, B.+x) --> true
    # @fact TestHelper.vec_eq(A.+x, A.+y) --> true
    # @fact TestHelper.vec_eq(A.+y, B.+y) --> true
    @fact TestHelper.vec_eq(x.+A, [1+x[1,1]  2+x[1,2];
                                   3+x[2,1]  4+x[2,2]]) --> true
    @fact TestHelper.vec_eq(x.+A, x.+B) --> true
    @fact TestHelper.vec_eq(x.+A, y.+A) --> true
    @fact TestHelper.vec_eq(x .+ x, [2x[1,1] 2x[1,2]; 2x[2,1] 2x[2,2]]) --> true
    # @fact TestHelper.vec_eq(y.+A, y.+B) --> true
    @fact TestHelper.vec_eq(A.-x, [1-x[1,1]  2-x[1,2];
                                   3-x[2,1]  4-x[2,2]]) --> true
    @fact TestHelper.vec_eq(A.-x, B.-x) --> true
    @fact TestHelper.vec_eq(A.-x, A.-y) --> true
    @fact TestHelper.vec_eq(x .- x, [zero(AffExpr) for _1 in 1:2, _2 in 1:2]) --> true
    # @fact TestHelper.vec_eq(A.-y, B.-y) --> true
    @fact TestHelper.vec_eq(x.-A, [-1+x[1,1]  -2+x[1,2];
                                   -3+x[2,1]  -4+x[2,2]]) --> true
    @fact TestHelper.vec_eq(x.-A, x.-B) --> true
    @fact TestHelper.vec_eq(x.-A, y.-A) --> true
    # @fact TestHelper.vec_eq(y.-A, y.-B) --> true
    @fact TestHelper.vec_eq(A.*x, [1*x[1,1]  2*x[1,2];
                                   3*x[2,1]  4*x[2,2]]) --> true
    @fact TestHelper.vec_eq(A.*x, B.*x) --> true
    @fact TestHelper.vec_eq(A.*x, A.*y) --> true
    # @fact TestHelper.vec_eq(A.*y, B.*y) --> true

    @fact TestHelper.vec_eq(x.*A, [1*x[1,1]  2*x[1,2];
                                   3*x[2,1]  4*x[2,2]]) --> true
    @fact TestHelper.vec_eq(x.*A, x.*B) --> true
    @fact TestHelper.vec_eq(x.*A, y.*A) --> true
    # @fact TestHelper.vec_eq(y.*A, y.*B) --> true

    @fact TestHelper.vec_eq(x .* x, [x[1,1]^2 x[1,2]^2; x[2,1]^2 x[2,2]^2]) --> true
    @fact_throws TestHelper.vec_eq(A./x, [1*x[1,1]  2*x[1,2];
                                          3*x[2,1]  4*x[2,2]])
    @fact TestHelper.vec_eq(x./A, [1/1*x[1,1]  1/2*x[1,2];
                                   1/3*x[2,1]  1/4*x[2,2]]) --> true
    @fact TestHelper.vec_eq(x./A, x./B) --> true
    @fact TestHelper.vec_eq(x./A, y./A) --> true
    # @fact TestHelper.vec_eq(A./y, B./y) --> true

    @fact TestHelper.vec_eq((2*x) / 3, full((2*y) / 3)) --> true
    @fact TestHelper.vec_eq(2 * (x/3), full(2 * (y/3))) --> true
    @fact TestHelper.vec_eq(x[1,1] * A, full(x[1,1] * B)) --> true
end

context("Vectorized comparisons") do
    m = Model()
    @variable(m, x[1:3])
    A = [1 2 3
         0 4 5
         6 0 7]
    B = sparse(A)
    @constraint(m, x'*A*x .>= 1)
    @fact TestHelper.vec_eq(m.quadconstr[1].terms, [x[1]*x[1] + 2x[1]*x[2] + 4x[2]*x[2] + 9x[1]*x[3] + 5x[2]*x[3] + 7x[3]*x[3] - 1]) --> true
    @fact m.quadconstr[1].sense --> :(>=)
    @constraint(m, x'*A*x .>= 1)
    @fact TestHelper.vec_eq(m.quadconstr[1].terms, m.quadconstr[2].terms) --> true

    mat = [ 3x[1] + 12x[3] +  4x[2]
            2x[1] + 12x[2] + 10x[3]
           15x[1] +  5x[2] + 21x[3]]

    @constraint(m, (x'A)' + 2A*x .<= 1)
    terms = map(v->v.terms, m.linconstr[1:3])
    lbs   = map(v->v.lb,    m.linconstr[1:3])
    ubs   = map(v->v.ub,    m.linconstr[1:3])
    @fact TestHelper.vec_eq(terms, mat) --> true
    @fact lbs --> fill(-Inf, 3)
    @fact ubs --> fill(   1, 3)
    @fact TestHelper.vec_eq((x'A)' + 2A*x, (x'A)' + 2B*x) --> true
    @fact TestHelper.vec_eq((x'A)' + 2A*x, (x'B)' + 2A*x) --> true
    @fact TestHelper.vec_eq((x'A)' + 2A*x, (x'B)' + 2B*x) --> true
    @fact TestHelper.vec_eq((x'A)' + 2A*x, @JuMP.Expression((x'A)' + 2A*x)) --> true
    @fact TestHelper.vec_eq((x'A)' + 2A*x, @JuMP.Expression((x'B)' + 2A*x)) --> true
    @fact TestHelper.vec_eq((x'A)' + 2A*x, @JuMP.Expression((x'A)' + 2B*x)) --> true
    @fact TestHelper.vec_eq((x'A)' + 2A*x, @JuMP.Expression((x'B)' + 2B*x)) --> true

    @constraint(m, -1 .<= (x'A)' + 2A*x .<= 1)
    terms = map(v->v.terms, m.linconstr[4:6])
    lbs   = map(v->v.lb,    m.linconstr[4:6])
    ubs   = map(v->v.ub,    m.linconstr[4:6])
    @fact TestHelper.vec_eq(terms, mat) --> true
    @fact lbs --> fill(-1, 3)
    @fact ubs --> fill( 1, 3)

    @constraint(m, -[1:3;] .<= (x'A)' + 2A*x .<= 1)
    terms = map(v->v.terms, m.linconstr[7:9])
    lbs   = map(v->v.lb,    m.linconstr[7:9])
    ubs   = map(v->v.ub,    m.linconstr[7:9])
    @fact TestHelper.vec_eq(terms, mat) --> true
    @fact lbs --> -[1:3;]
    @fact ubs --> fill( 1, 3)

    @constraint(m, -[1:3;] .<= (x'A)' + 2A*x .<= [3:-1:1;])
    terms = map(v->v.terms, m.linconstr[10:12])
    lbs   = map(v->v.lb,    m.linconstr[10:12])
    ubs   = map(v->v.ub,    m.linconstr[10:12])
    @fact TestHelper.vec_eq(terms, mat) --> true
    @fact lbs --> -[1:3;]
    @fact ubs --> [3:-1:1;]

    @constraint(m, -[1:3;] .<= (x'A)' + 2A*x .<= 3)
    terms = map(v->v.terms, m.linconstr[13:15])
    lbs   = map(v->v.lb,    m.linconstr[13:15])
    ubs   = map(v->v.ub,    m.linconstr[13:15])
    @fact TestHelper.vec_eq(terms, mat) --> true
    @fact lbs --> -[1:3;]
    @fact ubs --> fill(3,3)
end

end

facts("[operator] JuMPArray concatenation") do
    m = Model()
    @variable(m, x[1:3])
    @variable(m, y[1:3,1:3])
    @variable(m, z[1:1])
    @variable(m, w[1:1,1:3])

    @fact TestHelper.vec_eq([x y], [x[1] y[1,1] y[1,2] y[1,3]
                                    x[2] y[2,1] y[2,2] y[2,3]
                                    x[3] y[3,1] y[3,2] y[3,3]]) --> true

    @fact TestHelper.vec_eq([x 2y+1], [x[1] 2y[1,1]+1 2y[1,2]+1 2y[1,3]+1
                                       x[2] 2y[2,1]+1 2y[2,2]+1 2y[2,3]+1
                                       x[3] 2y[3,1]+1 2y[3,2]+1 2y[3,3]+1]) --> true

    @fact TestHelper.vec_eq([1 x'], [1 x[1] x[2] x[3]]) --> true
    @fact TestHelper.vec_eq([2x;x], [2x[1],2x[2],2x[3],x[1],x[2],x[3]]) --> true
    # vcat on JuMPArray
    @fact TestHelper.vec_eq([x;x], [x[1],x[2],x[3],x[1],x[2],x[3]]) --> true
    # hcat on JuMPArray
    @fact TestHelper.vec_eq([x x], [x[1] x[1]
                                    x[2] x[2]
                                    x[3] x[3]]) --> true
    # hvcat on JuMPArray
    tmp1 = [z w; x y]
    tmp2 = [z[1] w[1,1] w[1,2] w[1,3]
            x[1] y[1,1] y[1,2] y[1,3]
            x[2] y[2,1] y[2,2] y[2,3]
            x[3] y[3,1] y[3,2] y[3,3]]
    @fact TestHelper.vec_eq(tmp1, tmp2) --> true
    tmp3 = [1 2x'
            x 2y-x*x']
    tmp4 = [1    2x[1]               2x[2]               2x[3]
            x[1] -x[1]*x[1]+2y[1,1]  -x[1]*x[2]+2y[1,2]  -x[1]*x[3] + 2y[1,3]
            x[2] -x[1]*x[2]+2y[2,1]  -x[2]*x[2]+2y[2,2]  -x[2]*x[3] + 2y[2,3]
            x[3] -x[1]*x[3]+2y[3,1]  -x[2]*x[3]+2y[3,2]  -x[3]*x[3] + 2y[3,3]]
    @fact TestHelper.vec_eq(tmp3, tmp4) --> true

    A = sprand(3, 3, 0.2)
    B = full(A)
    @fact TestHelper.vec_eq([A y], [B y]) --> true
end

# The behavior in this test is no longer well-defined
#let
#    dot_m = Model()
#    @variable(dot_m, 0 <= v[1:3,1:4] <= 1)
#    @variable(dot_m, 0 <= w[3:2:7,7:-2:1] <= 1)
#    @fact string(dot(v,w)) --> "v[1,1]*w[3,7] + v[1,2]*w[3,5] + v[1,3]*w[3,3] + v[1,4]*w[3,1] + v[2,1]*w[5,7] + v[2,2]*w[5,5] + v[2,3]*w[5,3] + v[2,4]*w[5,1] + v[3,1]*w[7,7] + v[3,2]*w[7,5] + v[3,3]*w[7,3] + v[3,4]*w[7,1]"
#
#    @constraint(dot_m, sum(v) + sum(w) --> 2)
#    @objective(dot_m, :Max, v[2,3] + w[5,3])
#    solve(dot_m)
#    @fact dot(getvalue(v),getvalue(w)) --> 1.0
#end

# Ditto
#let
#    slice_m = Model()
#    C = [:cat,:dog]
#    @variable(slice_m, x[-1:1,C])
#    catcoef = [1,2,3]
#    dogcoef = [3,4,5]
#
#    @fact string(dot(catcoef, x[:,:cat])) --> "x[-1,cat] + 2 x[0,cat] + 3 x[1,cat]"
#    @fact string(dot(dogcoef, x[:,:dog])) --> "3 x[-1,dog] + 4 x[0,dog] + 5 x[1,dog]"
#end


text:
abstract type HamiltonianModel end

struct Ising <: HamiltonianModel
    Ni::Int
    Nj::Int
end
Ising() = Ising(1,1)

struct Ising22 <: HamiltonianModel
    Ni::Int
    Nj::Int
    r::Real
end
Ising22(r) = Ising22(2,2,r)

struct Ising33 <: HamiltonianModel
    Ni::Int
    Nj::Int
    r::Real
end
Ising33(r) = Ising33(3,3,r)

ext: jl


text:
type T1
   f1::Int64
   f2::Array{Float64,1}
end

function t1Cons1(a::Int64,n::Int64)
   T1(a,[i for i=1:n])
end

function t1Cons2(a::Int64,n::Int64)
   T1(a,rand(n))
end

function t1Mult(t::T1)
   return t.f1*t.f2;
end

max_stars_repo_name: Alfheim/node-julia

max_issues_repo_name: Alfheim/node-julia

max_forks_repo_name: Alfheim/node-julia


text:
"""
	MaxOutputs()

A model attribute a solver should use the specify the maximum number of scalar-valued
outputs it can handle in multiobjective optimization.
`MOI.get( model, MaxOutputs )` defaults to 1 and should be overwritten by multiobjective
solvers.
"""
struct MaxOutputs <: MOI.AbstractModelAttribute end
MOI.attribute_value_type( ::MaxOutputs ) = Int64

"""
    NumberOfObjectives()
A model attribute for the number of objective functions in the model.
As each objective function might have multiple outputs use `NumberOfOutputs()` to get the total number of scalar-valued outputs.
"""
struct NumberOfObjectives <: MOI.AbstractModelAttribute end
MOI.attribute_value_type( ::NumberOfObjectives ) = Int64

"""
	NumberOfOutputs()

A model attribute for the number of scalar-valued objectives in the model.

If, e.g., the model describes a multiobjective optimization problem and there
are two objective functions, the first one being scalar-valued and the second
returning a vector with 2 entries, then the NumberOfOutputs is 3.
"""
struct NumberOfOutputs <: MOI.AbstractModelAttribute end
MOI.attribute_value_type( ::NumberOfOutputs ) = Int64

"""
	ListOfObjectiveIndices()

A model attribute for the `Vector{ObjectiveIndex}` of all objective indices present in the model
(i.e., of length equal to the value of `NumberOfObjectives()`) in the order in
which they were added.
"""
struct ListOfObjectiveIndices <: MOI.AbstractModelAttribute end

"""
    ObjectiveValue(result_index::Int = 1, objective_index = nothing)
A model attribute for the objective value of the primal solution `result_index`.
See [`ResultCount`](@ref) for information on how the results are ordered.

If `objective_index` is an `ObjectiveIndex` then the value(s) of that objective
are concatenated and returned. If `objective_index` is `nothing` then the values of all objectives
are returned.
"""
struct ObjectiveValue{T<:Union{ObjectiveIndex,Nothing}} <: MOI.AbstractModelAttribute
    result_index::Int
	objective_index::T
    function ObjectiveValue(result_index::Int = 1, objective_index::Nothing=nothing)
		return new{Nothing}(result_index, objective_index)
	end
end

"""
    AbstractObjectiveAttribute
Abstract supertype for attribute objects that can be used to set or get attributes (properties) of objective functions in the model.
"""
abstract type AbstractObjectiveAttribute end

# NOTE: all the attributes below are `AbstractModelAttribute`s in MOI
# I want them to be `AbstractObjectiveAttribute`s instead.

"$(Docs.doc(MOI.ObjectiveSense()))"
struct ObjectiveSense <: AbstractObjectiveAttribute end
MOI.attribute_value_type(::ObjectiveSense) = MOI.OptimizationSense

"""
    ObjectiveFunction{F<:Union{AbstractScalarFunction, AbstractVectorFunction}}()
A model attribute for the objective function which has a type `F`.
`F` should be guaranteed to be equivalent but not necessarily identical to the function type provided by the user.
Throws an `InexactError` if the objective function cannot be converted to `F`,
e.g. the objective function is quadratic and `F` is `ScalarAffineFunction{Float64}` or
it has non-integer coefficient and `F` is `ScalarAffineFunction{Int}`.
"""
struct ObjectiveFunction{F<:MOI.AbstractFunction} <: AbstractObjectiveAttribute end
MOI.attribute_value_type(::ObjectiveFunction{F}) where {F} = F

"""
    ObjectiveFunctionType()
A model attribute for the type `F` of the objective function set using the
`ObjectiveFunction{F}` attribute.

## Examples
In the following code, `attr` should be equal to `MOI.VariableIndex`:
```julia
x = MOI.add_variable(model)
objf_1 = MOI.add_objective(model, x)
attr = MOI.get(model, MOI.ObjectiveFunctionType(), objf_1)
```
"""
struct ObjectiveFunctionType <: AbstractObjectiveAttribute end
MOI.attribute_value_type(::ObjectiveFunctionType) = Type{<:MOI.AbstractFunction}

struct NumberOfObjectiveOutputs <: AbstractObjectiveAttribute end
MOI.attribute_value_type(::NumberOfObjectiveOutputs) = Int64


text:
#!/usr/bin/env julia

const DATA_DIR   = "../17_data"
const INPUT_FILE = "input.txt"

@enum Activity active inactive

# ----------------------------------------
macro assert(bool_expr)
    message = string("Assertion: ", bool_expr, " failed")
    :($(esc(bool_expr)) || error($message))
end
# ----------------------------------------

struct BBox3
    xmin::Int
    xmax::Int
    ymin::Int
    ymax::Int
    zmin::Int
    zmax::Int
    # --------------------
    function BBox3(grid)
        xmin = xmax = nothing
        ymin = ymax = nothing
        zmin = zmax = nothing
        for p in keys(grid)
            if grid[p] == inactive
                continue
            end
            x = p[1]
            y = p[2]
            z = p[3]
            if isnothing(xmin) || x < xmin
               xmin = x
            end
            if isnothing(xmax) || x > xmax
               xmax = x
            end
            if isnothing(ymin) || y < ymin
               ymin = y
            end
            if isnothing(ymax) || y > ymax
               ymax = y
            end
            if isnothing(zmin) || z < zmin
               zmin = z
            end
            if isnothing(zmax) || z > zmax
               zmax = z
            end
        end
        return new(xmin, xmax, ymin, ymax, zmin, zmax)
    end
end

struct BBox4
    xmin::Int
    xmax::Int
    ymin::Int
    ymax::Int
    zmin::Int
    zmax::Int
    wmin::Int
    wmax::Int
    # --------------------
    function BBox4(grid)
        xmin = xmax = nothing
        ymin = ymax = nothing
        zmin = zmax = nothing
        wmin = wmax = nothing

        for p in keys(grid)
            if grid[p] == inactive
                continue
            end
            x = p[1]
            y = p[2]
            z = p[3]
            w = p[4]

            (isnothing(xmin) || x < xmin) && (xmin = x)
            (isnothing(xmax) || x > xmax) && (xmax = x)

            (isnothing(ymin) || y < ymin) && (ymin = y)
            (isnothing(ymax) || y > ymax) && (ymax = y)

            (isnothing(zmin) || z < zmin) && (zmin = z)
            (isnothing(zmax) || z > zmax) && (zmax = z)

            (isnothing(wmin) || w < wmin) && (wmin = w)
            (isnothing(wmax) || w > wmax) && (wmax = w)
        end
        return new(xmin, xmax, ymin, ymax, zmin, zmax, wmin, wmax)
    end
end

# ----------------------------------------

function Base.show(io::IO, bb::BBox3)
    print(io, "BBox3: $(bb.xmin)->x->$(bb.xmax), $(bb.ymin)->y->$(bb.ymax), $(bb.zmin)->z->$(bb.zmax)")
end

function Base.show(io::IO, bb::BBox4)
    print(io, "BBox4: $(bb.xmin)->x->$(bb.xmax), $(bb.ymin)->y->$(bb.ymax), $(bb.zmin)->z->$(bb.zmax), $(bb.wmin)->w->$(bb.wmax)")
end

function banner(txt::String)
    println('='^10, ' ', txt, ' ', '='^10)
end

const dxyzs = [ [dx,dy,dz]
                for dx in -1:1
                for dy in -1:1
                for dz in -1:1
                if !(dx==dy==dz==0)
              ]

const dps = [ [dx,dy,dz,dw]
              for dx in -1:1
              for dy in -1:1
              for dz in -1:1
              for dw in -1:1
              if !(dx==dy==dz==dw==0)
            ]

function get_states_count(grid)
    on  = length([1 for (p, state) in grid if state == active])
    off = length([1 for (p, state) in grid if state == inactive])
    return (on, off)
end

function get_nbr_count(grid, p, deltas)
    nbrs = [p+dp for dp in deltas if p+dp in keys(grid) && grid[p+dp] == active]
    return length(nbrs)
end

function get_nextgrid(prevgrid, deltas)
    states_prev = get_states_count(prevgrid)

    nextgrid = deepcopy(prevgrid)
    states_next = get_states_count(nextgrid)

    # ==================== Expand ====================
    next_keys = [p for (p, state) in nextgrid if state == active]
    for p in next_keys
        for dp in deltas
            if !(p+dp in keys(nextgrid))
                # println("@p=$(p): dp=$(dp). Adding $(p+dp)")
                nextgrid[p+dp] = inactive
            end
        end
    end

    # ==================== Cycle ====================
    # Note: We already copied, so here we only track changes
    for p in keys(nextgrid)
        nbr_count = get_nbr_count(prevgrid, p, deltas)  # active(prevgrid) == active(nextgrid)
        if get_state(prevgrid, p) == active && (nbr_count < 2 || nbr_count > 3)
            # println("INFO: $(p) --> inactive ($(nbr_count) nbrs)")
            nextgrid[p] = inactive
        elseif get_state(prevgrid, p) == inactive && nbr_count == 3
            # println("INFO: $(p) --> active   ($(nbr_count) nbrs)")
            nextgrid[p] = active
        end
    end

    # ==================== Prune ====================
    result = Dict([(p, state) for (p, state) in nextgrid if state == active])
    return result
end

function get_state(grid, p)
    return (p in keys(grid)) ? grid[p] : inactive
end

function popcount(grid)
    return length([1 for (p, state) in grid if state == active])
end

function print_grid3(label, grid)
    bb = BBox3(grid)
    println("vvvvv $(label): BBox3=$(bb) vvvvv")
    for z in bb.zmin:bb.zmax
        println("----- $(label): z=$(z) -----")
        for y in bb.ymin:bb.ymax
            for x in bb.xmin:bb.xmax
                if [x,y,z] in keys(grid)
                    print(grid[[x,y,z]] == active ? "#" : ".")
                else
                    print(".")
                end
            end
            println()
        end
        println()
    end
    println("^^^^^^^^^^^^^^^^^^^^")
end

function print_grid4(label, grid)
    bb = BBox4(grid)
    println("vvvvv $(label): BBox=$(bb) vvvvv")
    for w in bb.wmin:bb.wmax
        for z in bb.zmin:bb.zmax
            println("----- $(label): z=$(z), w=$(w)  -----")
            for y in bb.ymin:bb.ymax
                for x in bb.xmin:bb.xmax
                    state = get_state(grid, [x,y,z,w])
                    print(state == active ? "#" : ".")
                end
                println()
            end
            println()
        end
    end
    println("^^^^^^^^^^^^^^^^^^^^")
end

# ----------------------------------------

function part1(lines::Array{String}, verbose=false)
    banner("PART 1")

    xyzs = [ [x, y, 0]
             for (y, row) in enumerate(lines)
             for (x, c) in enumerate(row)
             if c == '#'
           ]
    grid0 = Dict([(p, active) for p in xyzs])
    verbose && print_grid3("Original grid", grid0)

    grid1 = get_nextgrid(grid0, dxyzs)
    verbose && print_grid3("After 1 cycle", grid1)

    grid2 = get_nextgrid(grid1, dxyzs)
    verbose && print_grid3("After 2 cycles", grid2)

    grid3 = get_nextgrid(grid2, dxyzs)
    verbose && print_grid3("After 3 cycles", grid3)

    grid4 = get_nextgrid(grid3, dxyzs)
    verbose && print_grid3("After 4 cycles", grid4)

    grid5 = get_nextgrid(grid4, dxyzs)
    verbose && print_grid3("After 5 cycles", grid5)

    grid6 = get_nextgrid(grid5, dxyzs)
    verbose && print_grid3("After 6 cycles", grid6)

    println("Popcount after 6 cycles: $(popcount(grid6))")
end

function part2(lines::Array{String})
    banner("PART 2")

    ps = [[x,y,0,0] for (y, row) in enumerate(lines) for (x, c) in enumerate(row) if c == '#']
    grid0 = Dict([(p, active) for p in ps])
    # print_grid4("Original grid", grid0)

    grid1 = get_nextgrid(grid0, dps)
    grid2 = get_nextgrid(grid1, dps)
    grid3 = get_nextgrid(grid2, dps)
    grid4 = get_nextgrid(grid3, dps)
    grid5 = get_nextgrid(grid4, dps)
    grid6 = get_nextgrid(grid5, dps)
    println("Popcount after 6 cycles: $(popcount(grid6))")
end

function main(filename::String)
    lines::Array{String} = readlines(filename)
    part1(lines)
    println()
    part2(lines)
end

main("$DATA_DIR/$INPUT_FILE")

max_stars_repo_name: jaycoskey/AdventOfCode

max_issues_repo_name: jaycoskey/AdventOfCode

max_forks_repo_name: jaycoskey/AdventOfCode


text:
"Exception type based on PROJ error handling"
struct PROJError <: Exception
    msg::String
    function PROJError(msg)
        # reset PROJ's error stack on construction
        proj_errno_reset(C_NULL)
        new(msg)
    end
end

function Base.showerror(io::IO, err::PROJError)
    err = string("PROJError: ", err.msg)
    println(io, err)
end

"Custom error handler, automatically set with `proj_log_func`"
function log_func(user_data::Ptr{Cvoid}, level::Cint, msg::Cstring)
    if level == PJ_LOG_ERROR
        throw(PROJError(unsafe_string(msg)))
    end
    return C_NULL
end

"Prevent an error converting a null pointer to a string, returns `nothing` instead"
aftercare(x::Cstring) = x == C_NULL ? nothing : unsafe_string(x)
aftercare(x::Ptr{Cstring}) = unsafe_loadstringlist(x)


text:
# Methods for site counting in sequences using a naive method
# ============================================================
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

# Methods for a single sequence.

"""
    count_sites_naive{T<:Site,A<:FourBitAlphs}(::Type{T}, a::BioSequence{A})
"""
function count_sites_naive{T<:Site,A<:FourBitAlphs}(::Type{T}, a::BioSequence{A})
    k = 0
    @inbounds for idx in eachindex(a)
        k += issite(T, a[idx])
    end
    return k
end

count_sites_naive{T<:Union{Gap,Ambiguous},A<:TwoBitAlphs}(::Type{T}, a::BioSequence{A}) = 0
count_sites_naive{A<:TwoBitAlphs}(::Type{Certain}, a::BioSequence{A}) = length(a)

"""
    count_sites_naive{T<:Site,A<:DNA_OR_RNA}(::Type{T}, a::BioSequence{A}, b::BioSequence{A})
This method counts the number of sites between a pair of aligned sequences of
type `T`.
"""
function count_sites_naive{T<:Site,A<:Alphabet,B<:Alphabet}(::Type{T},
                                                       a::BioSequence{A},
                                                       b::BioSequence{B})
    k = 0
    @inbounds for idx in 1:min(endof(a), endof(b))
        k += issite(T, a[idx], b[idx])
    end
    return k
end

"""
    count_sites_naive{T<:Mutation,A<:DNA_OR_RNA}(::Type{T}, a::BioSequence{A}, b::BioSequence{A})
This method counts the number of sites between a pair of aligned sequences of
SiteCase `T`.
Since the types of SiteCase counted by this function may not always be
unambiguously determined, this method returns a tuple of the count as well as
the number of sites which could not be unambiguously determined.
This second count is important for some downstream purposes, for example
evolutionary/genetic distance computations in which pairwise deletion of
ambiguous sites is nessecery.
"""
function count_sites_naive{T<:Mutation,A<:Alphabet,B<:Alphabet}(::Type{T},
                                                             a::BioSequence{A},
                                                             b::BioSequence{B})
    k = 0
    j = 0
    @inbounds for idx in 1:min(endof(a), endof(b))
        an = a[idx]
        bn = b[idx]
        k += issite(T, an, bn)
        j += !issite(Certain, an, bn)
    end
    return k, j
end

count_sites_naive{T<:Union{Gap,Ambiguous},A<:TwoBitAlphs}(::Type{T}, a::BioSequence{A}, b::BioSequence{A}) = 0
count_sites_naive{A<:TwoBitAlphs}(::Type{Certain}, a::BioSequence{A}, b::BioSequence{A}) = min(length(a), length(b))


text:
mutable struct Submarine
  depth::Int
  position::Int
  aim::Int
end

function move(self::Submarine, instruction::String)
  action, n = split(strip(instruction))
  n = parse(Int, n)
  if action == "forward"
    self.position += n
  elseif action == "down"
    self.depth += n
  elseif action == "up"
    self.depth -= n
  else
    throw(RuntimeError("Wrong instruction"))
  end
end

function updated_movement(self::Submarine, instruction::String)
  action, n = split(strip(instruction))
  n = parse(Int, n)
  if action == "forward"
    self.position += n
    self.depth += self.aim * n
  elseif action == "down"
    self.aim += n
  elseif action == "up"
    self.aim-= n
  else
    throw(RuntimeError("Wrong instruction"))
  end
end

function codify_position(self::Submarine)
  return self.depth * self.position
end

function run(instructions, mode, init_depth=0, init_pos=0, init_aim=0)
  sub = Submarine(init_depth, init_pos, init_aim)
  for instruction in instructions
    if mode == 1
      move(sub, instruction)
    else
      updated_movement(sub, instruction)
    end
  end
  result = codify_position(sub)
  println("The result (in mode $mode) is $result")
end

test_input = ["forward 5", "down 5", "forward 8", "up 3", "down 8", "forward 2"]
test_output = run(test_input, 1)
test_output = run(test_input, 2)

open("input.txt", "r") do f
  data = [s for s in readlines(f)]
  run(data, 1)
  run(data, 2)
end


text:
type MultiProductMultiStepCommon
  rateTimes::Vector{Float64}
  evolution::EvolutionDescription
end

function MultiProductMultiStepCommon(rateTimes::Vector{Float64})
  n = length(rateTimes) - 1
  evolutionTimes = zeros(n)
  relevanceRates = Vector{Pair{Int, Int}}(n)
  for i in eachindex(evolutionTimes)
    evolutionTimes[i] = rateTimes[i]
    relevanceRates[i] = Pair(i, i+1)
  end

  evolution = EvolutionDescription(rateTimes, evolutionTimes, relevanceRates)

  return MultiProductMultiStepCommon(rateTimes, evolution)
end

clone(mpmsc::MultiProductMultiStepCommon) = MultiProductMultiStepCommon(copy(mpmsc.rateTimes), clone(mpmsc.evolution))

rate_times(mpms::MultiProductMultiStep) = mpms.common.rateTimes
get_evolution(mpms::MultiProductMultiStep) = mpms.common.evolution

type MultiStepInverseFloater <: MultiProductMultiStep
  common::MultiProductMultiStepCommon
  fixedAccruals::Vector{Float64}
  floatingAccruals::Vector{Float64}
  fixedStrikes::Vector{Float64}
  fixedMultipliers::Vector{Float64}
  floatingSpreads::Vector{Float64}
  paymentTimes::Vector{Float64}
  payer::Bool
  multiplier::Float64
  lastIndex::Int
  currentIndex::Int
end

function MultiStepInverseFloater(rateTimes::Vector{Float64},
                                fixedAccruals::Vector{Float64},
                                floatingAccruals::Vector{Float64},
                                fixedStrikes::Vector{Float64},
                                fixedMultipliers::Vector{Float64},
                                floatingSpreads::Vector{Float64},
                                paymentTimes::Vector{Float64},
                                payer::Bool = true)

  common = MultiProductMultiStepCommon(rateTimes)
  multiplier = payer ? -1.0 : 1.0

  return MultiStepInverseFloater(common, fixedAccruals, floatingAccruals, fixedStrikes, fixedMultipliers, floatingSpreads, paymentTimes,
                                payer, multiplier, length(rateTimes), 1)
end

type ExerciseAdapter{M <: MarketModelExerciseValue} <: MultiProductMultiStep
  common::MultiProductMultiStepCommon
  exercise::M
  numberOfProducts::Int
  isExerciseTime::BitArray{1}
  currentIndex::Int
end

function ExerciseAdapter(exercise::MarketModelExerciseValue, numberOfProducts::Int = 1)
  ex = clone(exercise)

  common = MultiProductMultiStepCommon(ex.evolution.rateTimes)
  isExerciseTime = ex.isExerciseTime

  return ExerciseAdapter(common, ex, numberOfProducts, isExerciseTime, 1)
end

number_of_products(::MultiStepInverseFloater) = 1

possible_cash_flow_times(msif::MultiStepInverseFloater) = msif.paymentTimes
possible_cash_flow_times(ea::ExerciseAdapter) = possible_cash_flow_times(ea.exercise)

max_number_of_cashflows_per_product_per_step(::ExerciseAdapter) = 1
max_number_of_cashflows_per_product_per_step(::MultiStepInverseFloater) = 1

reset!(msif::MultiStepInverseFloater) = msif.currentIndex = 1

function reset!(ea::ExerciseAdapter)
  reset!(ea.exercise)
  ea.currentIndex = 1

  return ea
end

function next_time_step!(msif::MultiStepInverseFloater, currentState::CurveState, numberCashFlowsThisStep::Vector{Int}, genCashFlows::Vector{Vector{MarketModelCashFlow}})
  liborRate = forward_rate(currentState, msif.currentIndex)
  inverseFloatingCoupon = max((msif.fixedStrikes[msif.currentIndex] - msif.fixedMultipliers[msif.currentIndex] * liborRate), 0.0) * msif.fixedAccruals[msif.currentIndex]
  floatingCoupon = (liborRate + msif.floatingSpreads[msif.currentIndex]) * msif.floatingAccruals[msif.currentIndex]

  genCashFlows[1][1].timeIndex = msif.currentIndex
  genCashFlows[1][1].amount = msif.multiplier * (inverseFloatingCoupon - floatingCoupon)

  numberCashFlowsThisStep[1] = 1
  msif.currentIndex += 1

  return (msif.currentIndex == msif.lastIndex)
end

function next_time_step!(ea::ExerciseAdapter, currentState::CurveState, numberCashFlowsThisStep::Vector{Int}, genCashFlows::Vector{Vector{MarketModelCashFlow}})
  fill!(numberCashFlowsThisStep, 0)
  isDone = false

  next_step!(ea.exercise, currentState)
  if ea.isExerciseTime[ea.currentIndex]
    cf = get_value(ea.exercise, currentState)
    numberCashFlowsThisStep[1] = 1
    genCashFlows[1][1] = cf
    isDone = true
  end

  ea.currentIndex += 1
  return isDone || ea.currentIndex == length(ea.isExerciseTime)+1
end

## Clone ##
clone(msif::MultiStepInverseFloater) = MultiStepInverseFloater(clone(msif.common), copy(msif.fixedAccruals), copy(msif.floatingAccruals), copy(msif.fixedStrikes),
                                        copy(msif.fixedMultipliers), copy(msif.floatingSpreads), copy(msif.paymentTimes), msif.payer, msif.multiplier,
                                        msif.lastIndex, msif.currentIndex)

clone(ea::ExerciseAdapter) = ExerciseAdapter{typeof(ea.exercise)}(clone(ea.common), clone(ea.exercise), ea.numberOfProducts, ea.isExerciseTime, ea.currentIndex)


text:
using CSV, DataFrames, Plots, LaTeXStrings, LUComptonScatter

df = DataFrame(CSV.File("examples/exampleScatter.csv"))
x = df[!,1]
y = df[!,2]
#inverse of fit function, f to energy
f_inverse(x) = (x + 0.1221) / 0.739
#calculate uncertainty from propogation
energy_uncertainty(x) = sqrt((0.02627 * x)^2 + 16.4079^2)
#get them uncerainty values in energy
σy = energy_uncertainty.(f_inverse.(y))

E, σE, m, σm, chiSq = scan_box(610.0, 700.0, 460.0, 540.0, x, y, σy, 0.08)

println("E = $E ± $σE; mₑ = $m ± $σm")

#data for fit plotting
xFit = LinRange(0, x[length(x)], 1000)
yFit = compton(xFit, [E, m])
chiOut = round(chiSq, digits=2)

#plot results
scatter(x, y, yerror = σy,
        label = "Data",
        xlabel = "Scattering Angle [degrees]",
        ylabel = "Energy [keV]",
        dpi = 500, size = (350, 350))
plot!(xFit, yFit, label = L"Fit: $\chi^2/\textrm{dof} = %$chiOut$")


text:
function main(n)
  n = round(Int, n / 2 * 2)
  a = rand(0:100_000_000, n, n)
  b = a
  c = a * b
  v = round(Int, n/2) + 1
  println(c[v, v])
end

max_stars_repo_name: terasakisatoshi/nimProgram

max_issues_repo_name: terasakisatoshi/nimProgram

max_forks_repo_name: terasakisatoshi/nimProgram


text:
include("DM1_1.jl")
include("DM2_1.jl")
include("DM1_2.jl")

#======================================================================#
function ReactiveGRASP(
    matrix, #matrice de taille m*n representant les condion de notre SPP
    cost, #liste des couts de taille n de notre SPP
    n, #nombre de variable de notre SPP
    m, #nombre de condtion de notre SPP
    ite, #nombre de tirage au sort de alpha avant recalcul des probabilités
    coupe, #nombre de coupe de notre segmen [0,1], et le nombre de alpha
    alphaset, #permet de passer une liste de alpha personnelle
    temps #ressource en seconde alloué a notre reactive-GRASP
    )

    #initialisation
    (p,nb_iteration,z_cumul,zBest,zWorst) = intialiser(matrix,cost, n, m, ite, coupe, alphaset)
    evol_p=Float64[]
    t=time()
    z_global = zeros(Int64, length(p))
    ite_global = zeros(Int64, length(p))

    while (time()-t <= temps)
        append!(evol_p,p) #MaJ de notre historique de probabilité
        cpt=1
        #On s'assure que chaque alpha s'exprime au moins une fois.
        for i in 1:length(p)
            (SOL,z, crts) = GRASP(cost, matrix, n, m, p[i])
            z_cumul[i] += z
        end #fin for
        while (cpt <= ite)

            #Selection du alpha
            prob = rand(Float64)
            alpha_choisit = choix_alpha(p,prob)
            #
            (SOL, z, crts) = GRASP(cost, matrix, n, m, p[alpha_choisit])
            #Amelioration
            (SOL, z) = exchange1_1(SOL,n,m,cost,crts,matrix)
            #Réinitialisation
            nb_iteration[alpha_choisit] += 1
            z_cumul[alpha_choisit] += z
            #Mise à jours
            if z > zBest
                zBest = z
            end #fin if
            if z < zWorst
                zWorst = z
            end #fin if
            #Incrément
            cpt = cpt+1
        end #Fin while

        recalcul_p!(p,z_cumul,zBest,zWorst,nb_iteration,evol_p)

        for i in 1:length(p)
            z_global[i] = z_global[i] + z_cumul[i]
            ite_global[i] = ite_global[i] + nb_iteration[i]
        end

        nb_iteration = ones(Int64, length(p))
        z_cumul = zeros(Int64, length(p))

    end #fin while

    moyenne_global=zeros(Float64,length(p))

    for i in 1:length(p)
        moyenne_global[i] = z_global[i]/ite_global[i]
    end
    zAvg = sum(moyenne_global)/length(p)
    println("zBest: ", zBest,"   zAvg:  ", zAvg, "    zWorst: ", zWorst)
    #println(evol_p) #A decommenter si l'on souhaite afficher l'évolution des probabilités
    return(evol_p)
end #fin reactive-GRASP

#=======================================#
function intialiser(matrix, cost, n, m, ite, coupe, alphaset) #initialisation des variables
    p = Float64[]
    if alphaset == 0 #On verifie si on posséde ou non un set
        if coupe >= 1 #On verifie que les coupes sont admissibles
            for i in 1:coupe
                push!(p,i/coupe)
            end #fin for
        end #fin if
    else #on initialise le alphaset
        p = alphaset
    end #fin if

    nb_iteration = ones(Int64, length(p)) #On s'assurera plus tard que tous les alpha s'esprime au moins une fois
    z_cumul = zeros(Int64, length(p)) #Sera une liste de stockage pour le recalcul des probabilités

    (SOL,z) =Glouton(cost, matrix, n, m) #Initialisation d'une solution de base
    zBest = z
    zWorst = z

    return(p,nb_iteration,z_cumul,zBest,zWorst)
end #fin initialiser

#=======================================#
function recalcul_p!(p,z_cumul,zBest,zWorst,nb_iteration,evol_p)
    #initialisation
    q=zeros(Float64,length(p))
    somme_q=0

    for i in 1:length(p)
        moyenne = (z_cumul[i]/nb_iteration[i])
        q[i] = (moyenne-zWorst)/(zBest-zWorst)
        somme_q += q[i]
    end

    for i in 1:length(q)
        p[i] = q[i]/somme_q
    end

    for i in 2:length(p)-1 #On additionne les proabilités entre elles
        p[i] += p[i-1]
    end
    p[length(p)] = 1 #On s'assure qu'il n'y ai pas d'erreur d'arrondie pour 1
end #fin recalcul_p!

#===========================================#
function choix_alpha(p,prob)
    for i in 1:length(p)
        if prob < p[i]
            return(i)
        elseif i==length(p)
            return(i)
        end #fin if
    end #fin for
end #fin choix_alpha


text:
import LibXLS
using Test

const DATA_FOLDER = joinpath(@__DIR__, "..", "data")
@assert isdir(DATA_FOLDER)

fp_book1 = joinpath(DATA_FOLDER, "book1.xls")
@assert isfile(fp_book1)

fp_book1_1904 = joinpath(DATA_FOLDER, "book1_1904.xls")
@assert isfile(fp_book1_1904)

fp_xlsx = joinpath(DATA_FOLDER, "blank.xlsx")
@assert isfile(fp_xlsx)

# Checks wether `matrix` equals `test_data`, where test_data is a vector of columns
function check_test_data(matrix, test_data::Vector)

    function size_of_data(d::Vector)
        isempty(d) && return (0, 0)
        return length(d[1]), length(d)
    end

    rows, cols = size_of_data(test_data)

    for row in 1:rows, col in 1:cols
        test_value = test_data[col][row]
        value = matrix[row, col]
        if ismissing(test_value) || ( isa(test_value, AbstractString) && isempty(test_value) )
            @test ismissing(value) || ( isa(value, AbstractString) && isempty(value) )
        else
            if isa(test_value, Float64)
                @test isapprox(value, test_value)
            else
                @test value == test_value
            end
        end
    end

    nothing
end

function debug_cell_data(cell::LibXLS.st_cell_data)
    println("Cell Data Debug")
    println("record = ", cell.id)
    println("xf = ", cell.xf)
    if cell.str != C_NULL
        println("str = ", unsafe_string(cell.str))
    end

    println("d = ", cell.d)
    println("l = ", Int(cell.l))
end

@testset "Workbook" begin

    @testset "Valid XLS file" begin
        wb = LibXLS.openxls(fp_book1)
        LibXLS.close(wb)
    end

    @testset "Invalid XLS file" begin
        @test_throws ErrorException LibXLS.openxls(fp_xlsx)
    end

    @testset "Workbook info" begin
        LibXLS.openxls(fp_book1) do wb
            @test !LibXLS.is1904(wb)
            @test LibXLS.sheetcount(wb) == 2
            @test LibXLS.sheetnames(wb) == [ "Plan1", "Plan2" ]
            @test LibXLS.sheetname(wb, 1) == "Plan1"
            @test LibXLS.sheetname(wb, 2) == "Plan2"
            @test LibXLS.sheetindex(wb, "Plan1") == 1
            @test LibXLS.sheetindex(wb, "Plan2") == 2

            # returns false, is that really it?
            # @test LibXLS.isvisible(wb, 1))
            # @test LibXLS.isvisible(wb, "Plan1")
        end

        LibXLS.openxls(fp_book1_1904) do wb
            @test LibXLS.is1904(wb)
            @test LibXLS.sheetcount(wb) == 2
            @test LibXLS.sheetnames(wb) == [ "Plan1", "Plan2" ]
        end
    end
end

@testset "Worksheet" begin
    LibXLS.openxls(fp_book1) do wb
        @testset "open/close" begin
            let
                ws = wb["Plan2"]

                # closing ws will cause segfault or errors
                # on the remaining tests
                # LibXLS.close(ws)
            end

            let
                ws = wb["Plan1"]
                ws = wb[1]
                ws = wb[2]
            end
        end

        @testset "bounds" begin
            @test_throws AssertionError wb["invalid_sheetname"]
            @test_throws AssertionError wb[0]
            @test_throws AssertionError wb[3]
        end

        @testset "size" begin
            @test size(wb[1]) == (6, 7)
            @test size(wb[2]) == (4, 6)
        end

        @testset "sheetname" begin
            ws = wb[1]
            @test LibXLS.sheetindex(ws) == 1
            @test LibXLS.sheetname(ws) == "Plan1"
        end

        @testset "row data" begin
            let
                ws = wb["Plan1"]
                @test ws[2, 2] == 1

                test_data = [ [missing for i in 1:6],
                              [missing, 1, 2, 3, missing, 5],
                              [missing, 1000.1, 1000.2, 1000.3, missing, 1000.5],
                              [missing, "abc", "def", "ghi", missing, "xyz"],
                              # [missing, Date(2018, 12, 1), Date(2018, 12, 31), Date(2019, 1, 1), missing, Date(2019, 2, 26)]
                            ]
                check_test_data(ws, test_data)
            end

            let
                ws = wb["Plan2"]
                @test ws[2,2] == "A"

                test_data = [ [missing for i in 1:4],
                              [missing, "A", "B", "C"],
                              ["A", 1, 0.2, 0.3],
                              ["B", 0.2, 1, 0.4],
                              ["C", 0.3, 0.4, 1]
                            ]
                check_test_data(ws, test_data)
            end
        end
    end

    LibXLS.openxls(fp_book1) do wb
        ws = wb["Plan2"]
        @test ws[2,2] == "A"
    end

    LibXLS.openxls(fp_book1) do wb
        ws1 = wb["Plan1"]
        ws2 = wb["Plan2"]
        @test ws1[2,2] == 1
        @test ws2[2,2] == "A"
    end
end


text:
#Importing Neccessary Packages
using DifferentialEquations,ModelingToolkit #Solving ODEs
using Turing,Distributions #Running statistics on MCMC
using CSV, DataFrames #Handle the data
using StatsPlots,ProgressMeter #Plotting and Monitoring

include("HelperFunctions.jl")
include("MCMCRun.jl")

cd("./UserFriendly") #Change this to your model's folder

#General workflow for MCMC parameter fitting

#=
1. Set up the ODE
2. Import and process the data
3. Transform the ODE output to match the data
4. Run the MCMC algorithm
5. Postprocessing
=#

###############################################################
                    # 1. Set up the ODE
###############################################################

#Some Model Information needed for MCMC
parNames = [:Kab, :Kbc, :Kad]
varNames = [:a, :b, :c,:d]
#Set the number of parameters and states
parNum = length(parNames)
varNum = length(varNames)

modelInfo = ModelInformation(parNames,varNames)

## Define any constants for the model
const V0 = 4/7 #Can be defined here or just as a number in the equations

#Write the rate equations
function Model(du,u,p,t)
  Kab, Kbc, Kad, = p
  a,b,c,d = u

  du[1] = V0*(10.0-a) - Kab*a - Kad*a^2 #a
  du[2] = -V0*b - Kbc*b + Kab*a #b
  du[3] = -V0*c + Kbc*b #c
  du[4] = -V0*d + (1/2)*Kad*a^2 #d
end

#Define Values for ODE model
p = rand(parNum) #Initial arameter values (in order they are defined above)
u0 = zeros(varNum) #Initial Conditions (number of Differential Equations)
tspan = (0.0,10.0) #Time (start, end)

#Contruct the ODE Problem
prob = ODEProblem(Model,u0,tspan,p)
alg = Vern7() #ODE solver

###############################################################
                    # 2. Import the data
###############################################################

#Read in the CSV
data = CSV.read("./VVdata.csv")

#Convert the data into a workable form
dataTransform = ConvertData(data)

###############################################################
        # 3. Transform the ODE output to match the data
###############################################################

#This function takes in the solved ODE, and outputs the observable species
function ObserveTransform(sol,measuredTime)
  #Create a container to hold the solution at the desired time points
  obversedSpecies = sol(measuredTime)
  #What species are being measured?
  return obversedSpecies
end

###############################################################
                    # 4. Run the MCMC algorithm
###############################################################

#Provide any prior knowledge for the parameters
priors = fill(FlatPos(0.0),length(prob.p))

#How many MCMC sample do you want
mcmcSamples = 100_000

#Gather all the information to one structure
sampleProblem = MCMCSetup(modelInfo,prob,alg,dataTransform,mcmcSamples,priors)

#Run the MCMC
result = MCMCRun(sampleProblem)

#Save the parameter chains and information about the MCMC run
chainParameters = DataFrame(result,:parameters)
CSV.write("./Parameters.csv",chainParameters)

chainInternals = DataFrame(result,:internals)
CSV.write("./Internals.csv",chainInternals)
###############################################################
                    # 5. Postprocessing
###############################################################

#True parameter values
parTure = [0.833, 1.667, 0.167] #Parameter values

#Make a plot of the chains and density plots
chainsPlot = plot(result)
savefig(chainsPlot,"./ChainsODE.pdf")

#Corner plot for correlations
corPlot = autocorplot(result)
savefig(corPlot,"./AutoCorr.pdf")

#Running average Plot
runAvePlot = meanplot(result)
savefig(runAvePlot,"./RunAve.pdf")

#Running average Plot
cornerPlot = corner(result)
savefig(cornerPlot,"./Corner.pdf")

#Create a plot fitting "best" parameter values
@df data plot(:Time,[:A,:B,:C,:D],markershape=:auto,layout=4,legend=false)

#Retrieve the chains and information about the chains
#Get best parameter set
bestPar = BestParSet(result,parNames) #ignore the std

#Rerun the problem with these parameters
newProb = remake(prob, p=bestPar)
newSol = solve(newProb,alg)

plot!(newSol,layout=4)
savefig("./DataFit.pdf")

max_stars_repo_name: RobertGregg/TuringODE

max_issues_repo_name: RobertGregg/TuringODE

max_forks_repo_name: RobertGregg/TuringODE


text:
using WienerHopf, ApproxFun, RiemannHilbert, SingularIntegralEquations, Test
import ApproxFun: mappoint, setcanonicaldomain
import RiemannHilbert: orientedleftendpoint, fpcauchymatrix, fpstieltjesmatrix, collocationpoints, evaluationmatrix!, undirected, evaluationmatrix
import WienerHopf: fourtoone, ifourtoone, sqrtline_tocanonical

@testset "SqrtLine" begin
    @test ifourtoone(0.001-eps()) ≈ ifourtoone(0.001+eps()) # taylor v explicit
    @test ifourtoone(-0.001-eps()) ≈ ifourtoone(-0.001+eps()) # taylor v explicit
    @test ifourtoone(-0.001im-eps()im) ≈ ifourtoone(-0.001im+eps()im) # taylor v explicit
    @test sqrtline_tocanonical(-eps()) ≈ sqrtline_tocanonical(eps()) atol=1E-15
    @test sqrtline_tocanonical(-eps()) ≈ sqrtline_tocanonical(0) atol=1E-15
    @test 0 ≈ sqrtline_tocanonical(0) atol=1E-15

    @test ifourtoone(fourtoone(0.1)) ≈ ifourtoone(1,fourtoone(0.1)) ≈ 0.1
    @test ifourtoone(fourtoone(-0.1)) ≈ ifourtoone(1,fourtoone(-0.1)) ≈ -0.1
    @test fourtoone(ifourtoone(1,0.1+0.1im)) ≈ fourtoone(ifourtoone(2,0.1+0.1im)) ≈
        fourtoone(ifourtoone(3,0.1+0.1im)) ≈ fourtoone(ifourtoone(4,0.1+0.1im)) ≈ 0.1+0.1im
        @test fourtoone(ifourtoone(1,-0.1-0.1im)) ≈ fourtoone(ifourtoone(2,-0.1-0.1im)) ≈
        fourtoone(ifourtoone(3,-0.1-0.1im)) ≈ fourtoone(ifourtoone(4,-0.1-0.1im)) ≈ -0.1-0.1im

    @test angle(fromcanonical(SqrtLine{1/4}(), 0.1)) ≈ 0.7853981633974483
    @test angle(fromcanonical(SqrtLine{1/4}(), -0.1)) ≈ 0.7853981633974483-π
    @test isreal(tocanonical(SqrtLine{1/4}(), 0.1exp(im*π/4)))
    @test isreal(tocanonical(SqrtLine{1/4}(), -0.1exp(im*π/4)))

    f = Fun(x -> 1/sqrt(x-im) + 1/(x-im) + 1/sqrt(x+im-1) + 1/(x+im-1), SqrtLine())
    @test ncoefficients(f) ≤ 250 # spectral convergence
    @test f(0.1) ≈ 1/sqrt(0.1-im) + 1/(0.1-im) + 1/sqrt(0.1+im-1) + 1/(0.1+im-1)
end

@testset "SqrtLine Cauchy" begin
    f = Fun(x -> 1/sqrt(x-im) + 1/(x-im) + 1/sqrt(x+im-1) + 1/(x+im-1), SqrtLine())
    @test cauchy(f, 1+im) ≈ 1/sqrt((1+im)+im-1) + 1/((1+im)+im-1)
    @test cauchy(f, -1+im) ≈ 1/sqrt((-1+im)+im-1) + 1/((-1+im)+im-1)
    @test cauchy(f, 1-im) ≈ -1/sqrt((1-im)-im) - 1/((1-im)-im)
    @test cauchy(f, -1-im) ≈ -1/sqrt((-1-im)-im) - 1/((-1-im)-im)

    f = Fun(x -> 1/sqrt(x-im) + 1/(x-im) + 1/sqrt(x+im-1) + 1/(x+im-1), Legendre(SqrtLine()))
    @test cauchy(f, 1+im) ≈ 1/sqrt((1+im)+im-1) + 1/((1+im)+im-1)
    @test cauchy(f, 1-im) ≈ -1/sqrt((1-im)-im) - 1/((1-im)-im)

    C = fpcauchymatrix(space(f), ncoefficients(f), ncoefficients(f))
    pts = collocationpoints(space(f), ncoefficients(f))
    @test (C*f.coefficients)[5] ≈ cauchy(f,pts[5]-eps()im)

    k = 1
    γ = α -> im*sqrt(k-α)*sqrt(α+k)
    θ₀ = 0.1
    f = Fun(α -> -k*sin(θ₀)/(γ(α)*(α-k*cos(θ₀))), SqrtLine{-1/4}())
    @test ncoefficients(f) ≤ 300
    @test f(exp(-im*π/4)) ≈ (α->(-k*sin(θ₀)/(γ(α)*(α-k*cos(θ₀)))))(exp(-im*π/4))
    @test cauchy(f,exp(-im*π/4)+10eps())-cauchy(f,exp(-im*π/4)-10eps()) ≈ f(exp(-im*π/4))

    S = Legendre(SqrtLine{-1/4}())
    f = Fun(α -> -k*sin(θ₀)/(γ(α)*(α-k*cos(θ₀))), S)
    C = fpcauchymatrix(S, ncoefficients(f), ncoefficients(f))
    @test norm(C) ≤ 1500
    pts = collocationpoints(space(f), ncoefficients(f))
    @test (C*f.coefficients)[50] ≈ cauchy(f,pts[50]-10eps()im) ≈ cauchy(f,pts[50]⁻)
    @test (C*f.coefficients)[100] ≈ cauchy(f,pts[100]-10eps()im) ≈ cauchy(f,pts[100]⁻)
    @test (C*f.coefficients)[229] ≈ cauchy(f,pts[229]-10eps()im) ≈ cauchy(f,pts[229]⁻)
end

@testset "RHP" begin
    S = Legendre(SqrtLine{-1/4}())
    n = 256
    pts = collocationpoints(S, n)
    Cm = fpcauchymatrix(S,n,n)[2:end-1,:]
    E = evaluationmatrix(S,pts[2:end-1],n)
    Cp = E + Cm
    k = 1
    θ₀ = 0.1
    γ = α -> isinf(α) ? complex(Inf) : im*sqrt(k-α)*sqrt(α+k)
    f = α -> isinf(α) ? zero(α) : -k*sin(θ₀)/(γ(α)*(α-k*cos(θ₀)))

    Hp = α -> k*sin(θ₀)/(α-k*cos(θ₀)) * (1/sqrt(α+k) - 1/sqrt(k+k*cos(θ₀)))
    Hm = α -> k*sin(θ₀) / (sqrt(k+k*cos(θ₀))*(α-k*cos(θ₀)))

    Φdp = α -> -sqrt(α+k)Hp(α)
    Dm = α -> -Hm(α)/(im*sqrt(k-α))

    @test Φdp(0.0) + γ(0.0)Dm(0.0) ≈ -k*sin(θ₀)/(0.0-k*cos(θ₀))

    v = Fun(α -> Φdp(α*exp(-im*π/4))-Dm(α*exp(-im*π/4)), SqrtLine())
    @which evaluate(v.coefficients,space(v),0.0)
    @which tocanonical(space(v),0.000001)
    @test v(0.0) ≈ Φdp(0.0)-Dm(0.0)
    @test cauchy(v,0.1+10eps()im)-cauchy(v,0.1-10eps()im) ≈ v(0.1)
    @test cauchy(v,-0.1+10eps()im)-cauchy(v,-0.1-10eps()im) ≈ v(-0.1)
    @test cauchy(v,100eps()im)-cauchy(v,-100eps()im) ≈ v(0.0)

    v = Fun(α -> Φdp(α)-Dm(α), S)
    @test cauchy(v,2.0) ≈ Φdp(2.0)
    @test cauchy(v,-2.0) ≈ Dm(-2.0)

    @test Φdp(im) ≈ cauchy(v,im)
    h = 0.00001; @test Φdp(h) ≈ cauchy(v,h)
    h = 0.00001im; @test Φdp(h) ≈ cauchy(v,h)

    @test cauchy(v,pts[50]+100eps()im)-cauchy(v,pts[50]-100eps()im) ≈ v(pts[50])

    @test cauchy(v,0.0+10eps()im) + γ(0.0)*Dm(0.0) ≈ Φdp(0.0) + γ(0.0)*Dm(0.0) ≈ -k*sin(θ₀)/(0.0-k*cos(θ₀))
    @test (Cm*v.coefficients)[50] ≈ cauchy(v,pts[50]-100eps()im)
    @test (Cp*v.coefficients)[50] ≈ cauchy(v,pts[50]+100eps()im)

    L = Diagonal(inv.(γ.(pts[2:end-1])))*Cp .+ Cm
    norm(L*v.coefficients .- f.(pts[2:end-1]))
    u = [fill(1.0,1,n); L; ((-1.0).^(0:n-1))'] \ f.(pts)
    @test u ≈ v.coefficients

    L = Cp .+ Diagonal(γ.(pts[2:end-1]))*Cm
    g = α -> isinf(α) ? zero(α) : -k*sin(θ₀)/(α-k*cos(θ₀))
    u = [fill(1.0,1,n); L; ((-1.0).^(0:n-1))'] \ g.(pts)
    @test u ≈ v.coefficients

    L = E .+ Diagonal(γ.(pts[2:end-1]) .+ 1)*Cm
    u = [fill(1.0,1,n); L; ((-1.0).^(0:n-1))'] \ g.(pts)
    @test u ≈ v.coefficients
end


text:
using SimulateMotif
using Test

@testset "SimulateMotif.jl" begin
    # Write your tests here.
end


text:
import UncertainData: RandomSequences

x = rand(100)
y = rand(100)

test_ccm = ConvergentCrossMappingTest(timeseries_lengths = [45, 50], n_reps = 10)
test_cm = CrossMappingTest(n_reps = 10)
test_vf = VisitationFrequencyTest(binning = RectangularBinning(5), ηs = 1:5)
test_tog = TransferOperatorGridTest(binning = RectangularBinning(5), ηs = 1:5)
test_jdd = JointDistanceDistributionTest()
test_jddt = JointDistanceDistributionTTest()
predtest = VisitationFrequencyTest(binning = RectangularBinning(5), ηs = -5:5)
test_pa = PredictiveAsymmetryTest(predictive_test = predtest)
test_npa = NormalisedPredictiveAsymmetryTest(predtest, f = 1.0)
test_NNmi = NearestNeighbourMITest(ηs = -5:5)

Nseq = 80
nchunks = 2
rtest_cm = RandomSequencesTest(test_cm, RandomSequences(nchunks, Nseq))
rtest_ccm = RandomSequencesTest(test_ccm, RandomSequences(nchunks, Nseq))
rtest_vf = RandomSequencesTest(test_vf, RandomSequences(nchunks, Nseq))
rtest_tog = RandomSequencesTest(test_tog, RandomSequences(nchunks, Nseq))
rtest_jdd = RandomSequencesTest(test_jdd, RandomSequences(nchunks, Nseq))
rtest_jddt = RandomSequencesTest(test_jddt, RandomSequences(nchunks, Nseq))
rtest_pa = RandomSequencesTest(test_pa, RandomSequences(nchunks, Nseq))
rtest_npa = RandomSequencesTest(test_npa, RandomSequences(nchunks, Nseq))
rtest_NNmi = RandomSequencesTest(test_NNmi, RandomSequences(nchunks, Nseq))

@test causality(test_cm, x, y) isa CausalAnalysis
@test causality(test_ccm, x, y) isa CausalAnalysis
@test causality(test_vf, x, y) isa CausalAnalysis
@test causality(test_tog, x, y) isa CausalAnalysis
@test causality(test_jdd, x, y) isa CausalAnalysis
@test causality(test_jddt, x, y) isa CausalAnalysis
@test causality(test_pa, x, y) isa CausalAnalysis
@test causality(test_npa, x, y) isa CausalAnalysis
@test causality(test_NNmi, x, y) isa CausalAnalysis


text:
module UiPathOrchestratorJobSchedulingPlanCreate1

# Write your package code here.
"""
    createtimeset()

# 処理概要
時刻のデータセットを作成

# 引数
* なし

# 結果（戻り値）
* 時間の組み合わせ（文字列）
"""
end


text:
@testset "Core" begin
e2 = Edge(1,3)
e3 = Edge(1,4)
@test @inferred(is_ordered(e2))
@test @inferred(!is_ordered(reverse(e3)))

gx = Graph(10)
for g in testgraphs(gx)
    add_vertices!(g, 5)
    @test @inferred(nv(g)) == 15
end

g5w = WheelGraph(5)
for g in testgraphs(g5w)
    @test @inferred(indegree(g,1)) == @inferred(outdegree(g,1)) == @inferred(degree(g,1)) == 4
    @test @inferred(indegree(g)) == @inferred(outdegree(g)) == @inferred(degree(g)) == [4,3,3,3,3]

    @test @inferred(Δout(g)) == @inferred(Δin(g)) == @inferred(Δ(g)) == 4
    @test @inferred(δout(g)) == @inferred(δin(g)) == @inferred(δ(g)) == 3

    z1 = @inferred(degree_histogram(g))
    z2 = @inferred(degree_histogram(g, indegree))
    z3 = @inferred(degree_histogram(g, outdegree))

    @test z1 == z2 == z3 == Dict(4=>1, 3=>4)

    @test @inferred(neighbors(g, 2)) == @inferred(all_neighbors(g, 2)) == [1,3,5]
    @test @inferred(common_neighbors(g, 1, 5)) == [2, 4]

    gsl = copy(g)
    add_edge!(gsl, 3, 3)
    add_edge!(gsl, 2, 2)

    @test @inferred(!has_self_loops(g))
    @test @inferred(has_self_loops(gsl))
    @test @inferred(num_self_loops(g)) == 0
    @test @inferred(num_self_loops(gsl)) == 2

    @test @inferred(density(g)) == 0.8

    @test eltype(squash(g)) == UInt8
    end

    g5wd = WheelDiGraph(5)
    for g in testdigraphs(g5wd)
        @test @inferred(indegree(g,2)) == 2
        @test @inferred(outdegree(g,2)) == 1
        @test @inferred(degree(g,2)) == 3
        @test @inferred(indegree(g)) == [0,2,2,2,2]
        @test @inferred(outdegree(g)) == [4,1,1,1,1]
        @test @inferred(degree(g)) == [4,3,3,3,3]

        @test @inferred(Δout(g)) == @inferred(Δ(g)) == 4
        @test @inferred(Δin(g)) == 2
        @test @inferred(δout(g)) == 1
        @test @inferred(δin(g)) == 0
        @test @inferred(δ(g)) == 3

        z1 = @inferred(degree_histogram(g))
        z2 = @inferred(degree_histogram(g, indegree))
        z3 = @inferred(degree_histogram(g, outdegree))
        @test z1 == Dict(4=>1, 3=>4)
        @test z2 == Dict(0=>1, 2=>4)
        @test z3 == Dict(4=>1, 1=>4)

        @test @inferred(neighbors(g, 2)) == [3]
        @test Set(@inferred(all_neighbors(g, 2))) == Set([1,3,5])
        @test @inferred(common_neighbors(g, 1, 5)) == [2]

        gsl = copy(g)
        add_edge!(gsl, 3, 3)
        add_edge!(gsl, 2, 2)

        @test @inferred(!has_self_loops(g))
        @test @inferred(has_self_loops(gsl))
        @test @inferred(num_self_loops(g)) == 0
        @test @inferred(num_self_loops(gsl)) == 2

        @test @inferred(density(g)) == 0.4
        @test eltype(squash(g)) == UInt8
        end
    end

max_stars_repo_name: AsileBcd/JohnsonCycles

max_issues_repo_name: AsileBcd/JohnsonCycles

max_forks_repo_name: AsileBcd/JohnsonCycles


text:
import AstroTime: timescale

export AbstractCoord,
    AbstractState,
    AbstractTrajectory,
    centralbody,
    epoch,
    keplerian,
    position,
    state,
    refframe,
    timescale,
    velocity

abstract type AbstractCoord{S, F, B, T, N} <: AbstractArray{T, N} end
abstract type AbstractState{S, F, B, T} <: AbstractCoord{S, F, B, T, 1} end
abstract type AbstractTrajectory{S, F, B, T} <: AbstractCoord{S, F, B, T, 2} end

function epoch end

centralbody(::AbstractCoord{S, F, B}) where {S, F, B} = B
refframe(::AbstractCoord{S, F}) where {S, F} = F
timescale(::AbstractCoord{S}) where {S} = S
Base.eltype(::AbstractCoord{S, F, B, T}) where {S, F, B, T} = T

Base.size(::AbstractState) = (6,)
Base.length(::AbstractState) = 6
Base.IndexStyle(::Type{<:AbstractState}) = IndexLinear()


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    Parallel(transforms)

A transform where `transforms` are applied in parallel.

# Examples

```julia
Scale(low=0.3, high=0.6) ⊔ EigenAnalysis(:VDV)
ZScore() ⊔ EigenAnalysis(:V)
```
"""
struct Parallel <: Transform
  transforms::Vector{Transform}
end

isrevertible(p::Parallel) = any(isrevertible, p.transforms)

function apply(p::Parallel, table)
  # apply transforms in parallel
  f(transform) = apply(transform, table)
  vals = tcollect(f(t) for t in p.transforms)

  # retrieve tables and caches
  tables = first.(vals)
  caches = last.(vals)

  # table with concatenated columns
  newtable = tablehcat(tables)

  # find first revertible transform
  ind = findfirst(isrevertible, p.transforms)

  # save info to revert transform
  rinfo = if isnothing(ind)
    nothing
  else
    tcols  = Tables.columns.(tables)
    tnames = Tables.columnnames.(tcols)
    ncols  = length.(tnames)
    nrcols = ncols[ind]
    start  = sum(ncols[1:ind-1]) + 1
    finish = start + nrcols - 1
    range  = start:finish
    (ind, range)
  end

  newtable, (caches, rinfo)
end

function revert(p::Parallel, newtable, cache)
  # retrieve cache
  caches = cache[1]
  rinfo  = cache[2]

  @assert !isnothing(rinfo) "transform is not revertible"

  # retrieve info to revert transform
  ind    = rinfo[1]
  range  = rinfo[2]
  rtrans = p.transforms[ind]
  rcache = caches[ind]

  # columns of transformed table
  cols  = Tables.columns(newtable)
  names = Tables.columnnames(cols)

  # retrieve subtable to revert
  rcols  = [Tables.getcolumn(cols, j) for j in range]
  rnames = names[range]
  𝒯 = (; zip(rnames, rcols)...)
  rtable = 𝒯 |> Tables.materializer(newtable)

  # revert transform on subtable
  revert(rtrans, rtable, rcache)
end

function reapply(p::Parallel, table, cache)
  # retrieve caches
  caches = cache[1]

  # reapply transforms in parallel
  f(t, c) = reapply(t, table, c)
  itr     = zip(p.transforms, caches)
  tables  = tcollect(f(t, c) for (t, c) in itr)

  # table with concatenated columns
  tablehcat(tables)
end

function tablehcat(tables)
  # concatenate columns
  allvars, allvals = [], []
  varsdict = Set{Symbol}()
  for 𝒯 in tables
    cols = Tables.columns(𝒯)
    vars = Tables.columnnames(cols)
    vals = [Tables.getcolumn(cols, var) for var in vars]
    for (var, val) in zip(vars, vals)
      while var ∈ varsdict
        var = Symbol(var,:_)
      end
      push!(varsdict, var)
      push!(allvars, var)
      push!(allvals, val)
    end
  end

  # table with concatenated columns
  𝒯 = (; zip(allvars, allvals)...)
  𝒯 |> Tables.materializer(first(tables))
end

"""
    transform₁ ⊔ transform₂ ⊔ ⋯ ⊔ transformₙ

Create a [`Parallel`](@ref) transform with
`[transform₁, transform₂, …, transformₙ]`.
"""
⊔(t1::Transform, t2::Transform) = Parallel([t1, t2])
⊔(t1::Transform, t2::Parallel)  = Parallel([t1; t2.transforms])
⊔(t1::Parallel, t2::Transform)  = Parallel([t1.transforms; t2])
⊔(t1::Parallel, t2::Parallel)   = Parallel([t1.transforms; t2.transforms])


text:
#=
    Examples file
=#

using Nemo

include("../src/ExactODEReduction.jl")
using .ExactODEReduction: construct_jacobians, invariant_subspace, polynormalize

# we are using __randomized_wiedemann_minpoly to speed up computations, so,
# an example may fail occasionally

#------------------------------------------------------------------------------

# Example 1

R, (x₁, x₂, x₃, x₄) = QQ["x₁", "x₂", "x₃", "x₄"]

system = [
    x₁^2 + 2x₁*x₂,  # derivative of x1
    x₂^2 + x₃ + x₄, # derivative of x2
    x₂ + x₄,        # derivative of x3
    x₁ + x₃         # derivative of x4
]

# construct jacobians of the system to obtain matrices
matrices = construct_jacobians(system)

# find an invariant subspaces of the matrices
invariants = invariant_subspace(matrices)

println(invariants)

# transform it back ti polynomial form
transformation = polynormalize(invariants, R)

@assert transformation == [x₁ + x₂, x₃ + x₄]

#------------------------------------------------------------------------------

# Example 2
R, (x₁, x₂, x₃) = QQ["x₁", "x₂", "x₃"]

system = [
    x₂^2 + 4x₂*x₃ + 4x₃^2,  # derivative of x1
    4x₃ - 2x₁,              # derivative of x2
    x₁ + x₂                 # derivative of x3
]

# construct jacobians of the system to obtain matrices
matrices = construct_jacobians(system)

# find an invariant subspaces of the matrices
invariants = invariant_subspace(matrices)

# transform it back ti polynomial form
transformation = polynormalize(invariants, R)

# maybe it will be good to normalize, say,
# by the coef of the first variable
@assert transformation == [x₂ + 2x₃]

#------------------------------------------------------------------------------


text:
#####
##### `AbstractDifferential`
#####

"""
The subtypes of `AbstractDifferential` define a custom \"algebra\" for chain
rule evaluation that attempts to factor various features like complex derivative
support, broadcast fusion, zero-elision, etc. into nicely separated parts.

All subtypes of `AbstractDifferential` implement the following operations:

`+(a, b)`: linearly combine differential `a` and differential `b`

`*(a, b)`: multiply the differential `a` by the differential `b`

`Base.conj(x)`: complex conjugate of the differential `x`

`extern(x)`: convert `x` into an appropriate non-`AbstractDifferential` type for
use outside of `ChainContext`.

Valid arguments to these operations are `T` where `T<:AbstractDifferential`, or
where `T` has proper `+` and `*` implementations.

Additionally, all subtypes of `AbstractDifferential` support `Base.iterate` and
`Base.Broadcast.broadcastable(x)`.
"""
abstract type AbstractDifferential end

Base.:+(x::AbstractDifferential) = x

"""
    extern(x)

Return `x` converted to an appropriate non-`AbstractDifferential` type, for use
with external packages that might not handle `AbstractDifferential` types.

Note that this function may return an alias (not necessarily a copy) to data
wrapped by `x`, such that mutating `extern(x)` might mutate `x` itself.
"""
@inline extern(x) = x

@inline Base.conj(x::AbstractDifferential) = x

#####
##### `Wirtinger`
#####

"""
    Wirtinger(primal::Union{Number,AbstractDifferential},
              conjugate::Union{Number,AbstractDifferential})

Returns a `Wirtinger` instance representing the complex differential:

```
df = ∂f/∂z * dz + ∂f/∂z̄ * dz̄
```

where `primal` corresponds to `∂f/∂z * dz` and `conjugate` corresponds to `∂f/∂z̄ * dz̄`.

The two fields of the returned instance can be accessed generically via the
[`wirtinger_primal`](@ref) and [`wirtinger_conjugate`](@ref) methods.
"""
struct Wirtinger{P,C} <: AbstractDifferential
    primal::P
    conjugate::C
    function Wirtinger(primal::Union{Number,AbstractDifferential},
                       conjugate::Union{Number,AbstractDifferential})
        return new{typeof(primal),typeof(conjugate)}(primal, conjugate)
    end
end

wirtinger_primal(x::Wirtinger) = x.primal
wirtinger_primal(x) = x

wirtinger_conjugate(x::Wirtinger) = x.conjugate
wirtinger_conjugate(::Any) = Zero()

extern(x::Wirtinger) = throw(ArgumentError("`Wirtinger` cannot be converted to an external type."))

Base.Broadcast.broadcastable(w::Wirtinger) = Wirtinger(broadcastable(w.primal),
                                                       broadcastable(w.conjugate))

Base.iterate(x::Wirtinger) = (x, nothing)
Base.iterate(::Wirtinger, ::Any) = nothing

# TODO: define `conj` for` `Wirtinger`
Base.conj(x::Wirtinger) = throw(MethodError(conj, x))

#####
##### `Casted`
#####

"""
    Casted(v)

This differential wraps another differential (including a number-like type)
to indicate that it should be lazily broadcast.
"""
struct Casted{V} <: AbstractDifferential
    value::V
end

cast(x) = Casted(x)
cast(f, args...) = Casted(broadcasted(f, args...))

extern(x::Casted) = materialize(broadcasted(extern, x.value))

Base.Broadcast.broadcastable(x::Casted) = x.value

Base.iterate(x::Casted) = iterate(x.value)
Base.iterate(x::Casted, state) = iterate(x.value, state)

Base.conj(x::Casted) = cast(conj, x.value)

#####
##### `Zero`
#####

"""
    Zero()
The additive identity for differentials.
This is basically the same as `0`.
"""
struct Zero <: AbstractDifferential end

extern(x::Zero) = false  # false is a strong 0. E.g. `false * NaN = 0.0`

Base.Broadcast.broadcastable(::Zero) = Ref(Zero())

Base.iterate(x::Zero) = (x, nothing)
Base.iterate(::Zero, ::Any) = nothing

#####
##### `DNE`
#####

"""
    DNE()

This differential indicates that the derivative Does Not Exist (D.N.E).
This is not the cast that it is not implemented, but rather that it mathematically
is not defined.
"""
struct DNE <: AbstractDifferential end

function extern(x::DNE)
    throw(ArgumentError("Derivative does not exit. Cannot be converted to an external type."))
end

Base.Broadcast.broadcastable(::DNE) = Ref(DNE())

Base.iterate(x::DNE) = (x, nothing)
Base.iterate(::DNE, ::Any) = nothing

#####
##### `One`
#####

"""
     One()
The Differential which is the multiplicative identity.
Basically, this represents `1`.
"""
struct One <: AbstractDifferential end

extern(x::One) = true  # true is a strong 1.

Base.Broadcast.broadcastable(::One) = Ref(One())

Base.iterate(x::One) = (x, nothing)
Base.iterate(::One, ::Any) = nothing

#####
##### `Thunk`
#####

"""
    Thunk(()->v)
A thunk is a deferred computation.
It wraps a zero argument closure that when invoked returns a differential.
"""
struct Thunk{F} <: AbstractDifferential
    f::F
end

macro thunk(body)
    return :(Thunk(() -> $(esc(body))))
end

@inline extern(x::Thunk) = x.f()

Base.Broadcast.broadcastable(x::Thunk) = broadcastable(extern(x))

@inline function Base.iterate(x::Thunk)
    externed = extern(x)
    element, state = iterate(externed)
    return element, (externed, state)
end

@inline function Base.iterate(::Thunk, (externed, state))
    element, new_state = iterate(externed, state)
    return element, (externed, new_state)
end

Base.conj(x::Thunk) = @thunk(conj(extern(x)))


text:
#
# SIRD model input data
#

@with_kw struct SIRDInput

  #
  # SIRD simulation parameters
  #

  Si :: Float64 = 0.01 # Initial fraction of sick individuals
  Di :: Float64 = 0. # Initial fraction of dead indidivuals
  Ii :: Float64 = 0. # Initial fraction of immune individuals

  kc :: Float64  = 0.02 # Rate constant for contamination
  kd :: Float64 = 0.001 # Rate constant for deceases
  ki :: Float64 = 0.005 # Rate constant for immunization
  kiu :: Float64 = 0.00001 # Rate constant for immunity loss

  dt :: Float64 = 0.01 # time step
  tmax :: Float64 = 1000. # maximum time

  tol :: Float64 = 1.e-4 # stop if the fraction of sick people is smaller than this

  err :: Float64 = 1.e-3 # error tolerance for integrator

  nsave :: Int64 = 200 # Number of data points be saved

  print :: Bool  = true # Print or not data on the screen

end


text:
#!/usr/bin/env julia

import GR

function main()
    x = LinRange(0, 20, 50)
    y = x .^ 2

    GR.setviewport(0.1, 0.95, 0.1, 0.95)
    GR.setwindow(0, 20, 0.2, 400)
    GR.setscale(GR.OPTION_Y_LOG2)
    GR.settextfontprec(233, 3)

    GR.setmarkersize(1)
    GR.setmarkertype(GR.MARKERTYPE_SOLID_CIRCLE)
    GR.polymarker(x, y)

    GR.axes(2, 2, 0, 0.2, 1, 2, 0.005)

    GR.updatews()
end

main()


text:
let
    start_pt = VecE2(0.0,0.0)
    start_time = 0.0
    cell_width = 1.0
    transition_time = 1.0

    action_sequence = [WEST,WEST,NORTH,NORTH,WAIT,EAST,NORTH,EAST]
    grid_path = construct_grid_world_path(start_pt,start_time,
        action_sequence,cell_width,transition_time)
    traj = construct_trajectory(grid_path)
    verify(traj)

    controller = SwitchingController()
    state = [0.0,0.0,0.0]
    t = 0.1
    u = get_action(controller,traj,state,t)

    model = UnicycleModel()
    dt = 0.1
    simulate(model,controller,traj,state,0.0,2.0,dt)
end
# Test for singulatiries
let
    state = [0.0,0.0,0.0]
    target = state
    ff = [0.0,0.0]
    t = 0.0
    let
        controller = TrackingController()
        @test !any(isnan, get_action(controller,target,ff,state,t))
    end
    let
        controller = StabilizeController()
        @test !any(isnan, get_action(controller,target,ff,state,t))
    end
    let
        controller = PivotController()
        @test !any(isnan, get_action(controller,target,ff,state,t))
    end
end
let
    t1 = 0.0
    t2 = 1.0
    interval = TimeInterval(t1,t2)
    t = 0.5
    dt = get_Δt(interval,t)
    p1 = VecE2(0.0,0.0)
    p2 = VecE2(0.0,1.0)
    heading = (p2-p1)/norm(p2-p1)
    arclength = norm(p2-p1)
    seg1 = StraightTrajectory(p1,heading,arclength,TimeInterval(t1,t2))

    pos_vec = p1+dt*(p2-p1)
    velocity = (p2-p1)/(t2-t1)
    yaw_rate = 0.0

    t3 = 2.0
    t2_mid = (t3+t2)/2
    center = VecE2(0.0,1.0)
    radius = 0.0
    θ1 = 0.0
    Δθ = π/2
    θ2 = θ1 + Δθ
    seg2 = ArcTrajectory(center,p2,get_heading(seg1,get_end_time(seg1)),Δθ,TimeInterval(t2,t3))

    traj_ = Trajectory()
    push!(traj_,seg1)
    push!(traj_,seg2)
    verify(traj_)

    # Construct a DenseTrajectory that encodes the exact same information as the
    # piecewise constant Trajectory
    traj = Trajectory()
    for seg in traj_.segments
        t_vec = [get_start_time(seg), get_end_time(seg)]
        accel = [0.0]
        vel = [norm(get_vel(seg,get_start_time(seg))), norm(get_vel(seg,get_end_time(seg)))]
        pos = [0.0, get_length(seg)]
        push!(traj,DenseTrajectory(seg,t_vec,accel,vel,pos))
    end

    controller = SwitchingController()
    state = [0.001,0.001,0.001]
    t = 0.1
    u = get_action(controller,traj,state,t)

    model = UnicycleModel()
    dt = 0.1
    simulate(model,controller,traj,state,0.01,2.0,dt)
end
# Test Driving in reverse mode
let
    t1 = 0.0
    t2 = 1.0
    t = (t1+t2)/2.0
    p1 = VecE2(1.0,0.0)
    p2 = VecE2(0.0,1.0)
    heading = (p2-p1)/norm(p2-p1)
    arclength = norm(p2-p1)

    traj = StraightTrajectory(p1,heading,arclength,TimeInterval(t1,t2))
    reversed_traj = ReverseTrajectory(traj)

    controller = SwitchingController()

    state1 = [0.1, 0.5, 0.1]
    state2 = [0.1, 0.5, pi+0.1]
    a1 = get_action(controller,traj,state1,t)
    a2 = get_action(controller,reversed_traj,state2,t)
    @test isapprox(a1[1], a2[1])
    @test isapprox(a1[2],-a2[2])
end