text:
using Base.Test
using EGR

println("TestSAGA")

createOracleOutputLevel = 1
numEquivalentPasses = 10
algOutputLevel = 0
maxOutputNum=20
numChunks=5
stepOutputLevel = 0

sLikeGd(k,numTP) = k==0 ? 0 : numTP
uLikeGd(k,numTP) = k==0 ? numTP : 0

myREfunction(problem, opts, sd, expIndices,n) = false
myWriteFunction(problem, sd, opts, k, gnum, origWant, fromOutputsFunction) = false

for (gradientOracle, numVars, numTP, csDataType, LossFunctionString, myOutputter, L2reg, thisDataName, thisProblem) in Oracles

	println(" $thisDataName $LossFunctionString L2reg = $L2reg")

	maxG  = Int(round(numEquivalentPasses*numTP))

	expIndices=unique(round(Int, logspace(0,log10(maxG+1),maxOutputNum)))-1

	myOutputOpts =  OutputOpts(myOutputter,expIndices)
	myOptss = Opts(zeros(numVars,1); stepSizePower=-10, maxG=maxG, outputLevel=algOutputLevel, stepOutputLevel=0)

	for numChunks in [1,round(Int,numTP^(1/4)),round(Int,numTP^(2/4)),round(Int,numTP^(3/4)), numTP]

		println("numChunks = $numChunks")
		problem = thisProblem(Task(() -> getSequentialFinite(numTP, gradientOracle)))
		y=csDataType[]
		for i in 1:numChunks
			# println("Working on chunk $i")
			for j in (i-1)*round(Int,floor(numTP/numChunks))+1:i*round(Int,floor(numTP/numChunks))
				(f,sampleG) = (problem.getSampleFunctionAt(j))(myOptss.init)
				push!(y,sampleG)
			end
		end

		alg(problem,  myOptss, SAGA(numVars, y, numTP,numChunks,1), myOutputOpts, myWriteFunction)

	end
end

println("TestSAGA successful!")
println()


text:
function acrobot_benchmarks!(suite::BenchmarkGroup)

    T = Float64

    # options
    max_con_viol = 1.0e-8
    verbose=false

    opts_ilqr = iLQRSolverOptions{T}(verbose=verbose,
        live_plotting=:off)

    opts_al = AugmentedLagrangianSolverOptions{T}(verbose=verbose,
        opts_uncon=opts_ilqr,
        cost_tolerance=1.0e-5,
        cost_tolerance_intermediate=1.0e-2,
        constraint_tolerance=max_con_viol,
        penalty_scaling=100.,
        penalty_initial=1.)

    opts_pn = ProjectedNewtonSolverOptions{T}(verbose=verbose,
        feasibility_tolerance=max_con_viol)

    opts_altro = ALTROSolverOptions{T}(verbose=verbose,
        opts_al=opts_al,
        opts_pn=opts_pn,
        projected_newton=true,
        projected_newton_tolerance=1.0e-4);

    opts_ipopt = DIRCOLSolverOptions{T}(verbose=verbose,
        nlp=Ipopt.Optimizer(),
        feasibility_tolerance=max_con_viol)

    opts_snopt = DIRCOLSolverOptions{T}(verbose=verbose,
        nlp=SNOPT7.Optimizer(),
        feasibility_tolerance=max_con_viol)

    run_benchmarks!(suite, Problems.acrobot_problem, [opts_ilqr, opts_al, opts_altro, opts_ipopt, opts_snopt])
end

ext: jl

max_issues_repo_name: MehrdadZareian/GMARA

max_forks_repo_name: MehrdadZareian/GMARA


text:
import Dates
import TimeZones

"""
    datetime_intervals(input_dt_list::Vector{Dates.DateTime};
                       round_to::Dates.Period)

## Arguments
- `input_dt_list::Vector{Dates.DateTime}`: A vector of `DateTime`s.

## Keyword Arguments
- `round_to::Dates.Period`: Resolution to which all intervals should be rounded.

## Example

```jldoctest
julia> using Dates

julia> using DateShifting

julia> dates = [
           DateTime("2000-01-01T00:00:00"),
           DateTime("2000-02-01T00:00:00"),
           DateTime("2000-01-05T00:00:00"),
           DateTime("2000-01-02T04:05:06"),
           DateTime("2000-01-02T01:02:03"),
       ]
5-element Array{DateTime,1}:
 2000-01-01T00:00:00
 2000-02-01T00:00:00
 2000-01-05T00:00:00
 2000-01-02T04:05:06
 2000-01-02T01:02:03

julia> sequence, intervals = datetime_intervals(dates; round_to = Day(1))
([1, 5, 4, 3, 2], Dates.Day[0 days, 31 days, 4 days, 1 day, 1 day])

julia> sequence
5-element Array{Int64,1}:
 1
 5
 4
 3
 2

julia> intervals
5-element Array{Day,1}:
 0 days
 31 days
 4 days
 1 day
 1 day
```
"""
function datetime_intervals(input_dt_list::Vector{Dates.DateTime};
                            round_to::Dates.Period)
    input_zdt_list = [TimeZones.ZonedDateTime(x, Dates.TimeZone("UTC")) for x in input_dt_list]
    return datetime_intervals(input_zdt_list; round_to = round_to)
end

"""
    datetime_intervals(input_zdt_list::Vector{TimeZones.ZonedDateTime};
                       round_to::Dates.Period)

## Arguments
- `input_zdt_list::Vector{TimeZones.ZonedDateTime}`: A vector of `ZonedDateTime`s.

## Keyword Arguments
- `round_to::Dates.Period`: Resolution to which all intervals should be rounded.

## Example

```jldoctest
julia> using Dates

julia> using DateShifting

julia> using TimeZones

julia> dates = [
           ZonedDateTime(DateTime("2000-01-01T00:00:00"), tz"America/New_York"),
           ZonedDateTime(DateTime("2000-02-01T00:00:00"), tz"America/New_York"),
           ZonedDateTime(DateTime("2000-01-05T00:00:00"), tz"America/New_York"),
           ZonedDateTime(DateTime("2000-01-02T03:05:06"), tz"America/Chicago"),
           ZonedDateTime(DateTime("2000-01-02T01:02:03"), tz"America/New_York"),
       ]
5-element Array{ZonedDateTime,1}:
 2000-01-01T00:00:00-05:00
 2000-02-01T00:00:00-05:00
 2000-01-05T00:00:00-05:00
 2000-01-02T03:05:06-06:00
 2000-01-02T01:02:03-05:00

julia> sequence, intervals = datetime_intervals(dates; round_to = Day(1))
([1, 5, 4, 3, 2], Dates.Day[0 days, 31 days, 4 days, 1 day, 1 day])

julia> sequence
5-element Array{Int64,1}:
 1
 5
 4
 3
 2

julia> intervals
5-element Array{Day,1}:
 0 days
 31 days
 4 days
 1 day
 1 day
```
"""
function datetime_intervals(input_zdt_list::Vector{TimeZones.ZonedDateTime};
                            round_to::Dates.Period)
    sequence = sortperm(input_zdt_list)

    value, index_of_input_start_zoned_datetime = findmin(sequence)
    always_assert(value == 1)
    input_start_zoned_datetime = input_zdt_list[index_of_input_start_zoned_datetime]

    intervals_nonrounded = [_compute_interval_nonrounded(input_start_zoned_datetime, x) for x in input_zdt_list]
    intervals = [round(interval, round_to) for interval in intervals_nonrounded]

    return sequence, intervals
end


text:
using HypergeoMat
using Test
using LinearAlgebra

@testset "Hypergeo0F0 is exponential of trace" begin
  X1 = [0.3 0.2 0.1; 0.2 0.3 0.2; 0.1 0.2 0.3]
  X2 = [0.3im 0.2 0.1; 0.2 0.3im 0.2; 0.1 0.2 0.3im]
  @test isapprox(
    hypergeomPQ(10, Float64[], Float64[], X1),
    exp(tr(X1))
  )
  @test isapprox(
    hypergeomPQ(10, Float64[], Float64[], X2),
    exp(tr(X2))
  )
end

@testset "Hypergeo1F0 is det(I-X)^(-a)" begin
  X1 = [0.03 0.02 0.01; 0.02 0.03 0.02; 0.01 0.02 0.03]
  X2 = [0.02 0.02 0.01; 0.02 0.02 0.02; 0.01 0.02 0.02]
  X3 = [0.03im 0.02 0.01; 0.02 0.03im 0.02; 0.01 0.02 0.03im]
  @test isapprox(
    hypergeomPQ(15, [3.0], Float64[], X1),
    det(I - X1)^(-3)
  )
  @test isapprox(
    hypergeomPQ(15, [4.0im], Float64[], X2),
    det(I - X2)^(-4.0im)
  )
  @test isapprox(
    hypergeomPQ(15, [3.0], Float64[], X3),
    det(I - X3)^(-3)
  )
end

@testset "Hypergeo2F1" begin
  @testset "some values" begin
    @test isapprox(
      hypergeomPQ(10, [1.0; 2.0], [3.0], [0.2; 0.5]),
      1.79412894456143
    )
    @test isapprox(
      hypergeomPQ(10, [1.0im; 2.0], [3.0im], [0.2; 0.5]),
      1.677558924 - 0.183004016im
    )
    @test isapprox(
      hypergeomPQ(10, [1.0; 2.0], [3.0], [0.2im; 0.5]),
      1.513810425 + 0.20576184im
    )
    @test isapprox(
      hypergeomPQ(10, [1.0; 2.0im], [3.0], [0.2im; 0.5]),
      0.7733140719 + 0.3092059749im
    )
  end
  @testset "Gauss formula" begin
    a = 1.1
    b = 2.0im
    c = 9.0
    @test isapprox(
      hypergeomPQ(100, [a; b], [c], [1.0, 1.0, 1.0]),
      mvgamma(c, 3) * mvgamma(c-a-b, 3) / mvgamma(c-a, 3) / mvgamma(c-b, 3)
    )
  end
end

@testset "Hypergeo3F2: Kummer relation" begin
  a = [1.0; 2.0; 3.0]
  b = [9.0; 10.0]
  c = sum(b) - sum(a)
  p = 4
  K = mvgamma(b[2],p) * mvgamma(c,p) / mvgamma(b[2]-a[3],p) / mvgamma(c+a[3],p)
  o1 = K *
    hypergeomPQ(100, [b[1]-a[1]; b[1]-a[2]; a[3]], [b[1]; c+a[3]], [1.0; 1.0; 1.0; 1.0])
  o2 = hypergeomPQ(15, a, b, [1.0; 1.0; 1.0; 1.0])
  @test isapprox(o1, o2, rtol = 1e-3)
end


text:
module BLIS

using Base
using Libdl
using blis_jll
using LinearAlgebra

global libblis = C_NULL

__init__() = begin
    if length(get(ENV, "BLISDIR", "")) > 0
        # BLIS installation overriden by environmental variables.
        @info "Using custom defined BLIS installation instead of blis_jll."
        global libblis = dlopen(string(get(ENV, "BLISDIR", ""), "/lib/libblis"))
    else
        blis_path = blis_jll.blis_path
        # Use BinaryBuilder provided BLIS library.
        @info "blis_jll yields BLIS installation: $blis_path."
        global libblis = dlopen(blis_path)
    end
end

# Data types.
module Types
include("types.jl")
end

# Typed-API backend.
module TypedBackend
import ..libblis
using ..Libdl
using ..Types
using ..Types: cblas_typechar
using ..Types: BliCompatibleType
using ..Types: BliDim, BliInc, BliDoff, BliSiz
using ..Types: BliTrans, BliConj, BliUpLo, BliSide, BliDiag, BliInvDiag, BliStruc
include("backend_typed/common.jl")
include("backend_typed/level1v.jl")
include("backend_typed/level1d.jl")
include("backend_typed/level1m.jl")
include("backend_typed/level1f.jl")
include("backend_typed/level2.jl")
include("backend_typed/level3.jl")
include("backend_typed/utility.jl")
end

# Object-API backend.
module ObjectBackend
import ..libblis
using ..Libdl
using ..Base: RefValue, unsafe_convert
using ..Types: ctype_to_bli_num
using ..Types: BliCompatibleType
using ..Types: BliDim, BliInc, BliDoff, BliSiz
using ..Types: BliObjBits, BliAtomic, BliNum
using ..Types: BliTrans, BliConj, BliUpLo, BliSide, BliDiag, BliInvDiag, BliStruc
using ..Types: BLIS_CONJTRANS_BITS, BLIS_DATATYPE_BITS, BLIS_DOMAIN_BIT, BLIS_PRECISION_BIT
using ..Types: BLIS_CONJTRANS_BITS, BLIS_TRANS_BIT, BLIS_CONJ_BIT, BLIS_UPLO_BITS
using ..Types: BLIS_UPPER_BIT, BLIS_DIAG_BIT, BLIS_LOWER_BIT, BLIS_UNIT_DIAG_BIT
using ..Types: BLIS_INVERT_DIAG_BIT, BLIS_STRUC_BITS
include("backend_object/object.jl")
include("backend_object/common.jl")
include("backend_object/level1v.jl")
include("backend_object/level1d.jl")
include("backend_object/level1m.jl")
include("backend_object/level1f.jl")
include("backend_object/level2.jl")
include("backend_object/level3.jl")
include("backend_object/utility.jl")
end

# LinearAlgebra BLAS interface.
module BLASInterface
import ..LinearAlgebra.BLAS: gemm!, hemm!, symm!, trmm!, trsm!, herk!, syrk!, her2k!, syr2k!
import ..LinearAlgebra.BLAS: gemv!, hemv!, symv!, trmv!, trsv!, ger!, her!, syr!
import ..LinearAlgebra.BLAS: axpby!, axpy!, scal!, blascopy!
import ..LinearAlgebra: lapack_size, gemm_wrapper!
import ..LinearAlgebra: Adjoint, StridedVecOrMat
export gemm!, hemm!, symm!, trmm!, trsm!, herk!, syrk!, her2k!, syr2k!
export gemv!, hemv!, symv!, trmv!, trsv!, ger!, her!, syr!
export axpby!, axpy!, scal!, blascopy!
export gemm_wrapper!
using ..TypedBackend
using ..ObjectBackend
using ..Types
using ..Types: cblas_typechar
using ..Types: char_to_trans, char_to_conj, char_to_side, char_to_uplo, char_to_diag
using ..Types: BliTrans, BliConj, BliUpLo, BliSide, BliDiag, BliInvDiag, BliStruc
using ..Types: BliCompatibleType
using ..Types: BLIS_TRANS_BIT
using DelimitedFiles
# Load blacklist.
blacklist = []
if isfile(string(ENV["HOME"], "/.blis_jlbla_blacklist"))
    blacklist = readdlm(string(ENV["HOME"], "/.blis_jlbla_blacklist"))
end
include("interface_linalg/level3.jl")
include("interface_linalg/level2.jl")
include("interface_linalg/level1.jl")
end

end


text:
module Keldysh

# utility functions
export fermi

# Branch functions
export BranchEnum, forward_branch, backward_branch, imaginary_branch, BranchPoint, Branch

# Contour functions
export ContourEnum, full_contour, keldysh_contour, imaginary_contour, Contour, twist, heaviside, θ, get_branch, nbranches

# TimeGrid functions
export TimeGrid, TimeGridPoint, branch_bounds, integrate, step, realtimes, imagtimes

# TimeGF functions
export TimeGF, TimeGFTranspose, jump

# DOS functions
export dos_integrator, flat_dos, gaussian_dos, bethe_dos, chain_dos, square_dos

# generate_gf functions
export gf_1level, dos2gf

import Base.in, Base.size, Base.getindex, Base.setindex!, Base.IndexStyle, Base.step, Base.length, Base.similar, Base.adjoint, Base.read, Base.write

include("util.jl")
include("branch.jl")
include("contour.jl")
include("time_grid.jl")
include("time_gf.jl")
include("dos.jl")
include("generate_gf.jl")
include("hdf5.jl")
include("nessi.jl")

end # module


text:
"""
    abstract type AbstractPolynomialBasis end

Polynomial basis of a subspace of the polynomials [Section~3.1.5, BPT12].

[BPT12] Blekherman, G.; Parrilo, P. A. & Thomas, R. R.
*Semidefinite Optimization and Convex Algebraic Geometry*.
Society for Industrial and Applied Mathematics, **2012**.
"""
abstract type AbstractPolynomialBasis end

function MP.polynomial(coefs::Vector, basis::AbstractPolynomialBasis)
    return MP.polynomial(i -> coefs[i], basis)
end

"""
    maxdegree_basis(B::Type{<:AbstractPolynomialBasis}, variables, maxdegree::Int)

Return the basis of type `B` generating all polynomials of degree up to
`maxdegree` with variables `variables`.
"""
function maxdegree_basis end

"""
    basis_covering_monomials(B::Type{<:AbstractPolynomialBasis}, monos::AbstractVector{<:AbstractMonomial})

Return the minimal basis of type `B` that can generate all polynomials of the
monomial basis generated by `monos`.

## Examples

For example, to generate all the polynomials with nonzero coefficients for the
monomials `x^4` and `x^2`, we need three polynomials as otherwise, we generate
polynomials with nonzero constant term.
```jldoctest
julia> using DynamicPolynomials

julia> @polyvar x
(x,)

julia> basis_covering_monomials(ChebyshevBasis, [x^4, x^2])
ChebyshevBasisFirstKind{Polynomial{true,Float64}}(DynamicPolynomials.Polynomial{true,Float64}[8.0x⁴ - 8.0x² + 1.0, 2.0x² - 1.0, 1.0])
```
"""
function basis_covering_monomials end

max_stars_repo_name: JuliaAlgebra/PolynomialBases

max_issues_repo_name: JuliaAlgebra/PolynomialBases

max_forks_repo_name: JuliaAlgebra/PolynomialBases


text:
export iauLtecm
"""
ICRS equatorial to ecliptic rotation matrix, long-term.

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

Status:  support function.

Given:
   epj     double         Julian epoch (TT)

Returned:
   rm      double[3][3]   ICRS to ecliptic rotation matrix

Notes:

   1. The matrix is in the sense

      E_ep = rm x P_ICRS,

   where P_ICRS is a vector with respect to ICRS right ascension
   and declination axes and E_ep is the same vector with respect to
   the (inertial) ecliptic and equinox of epoch epj.

   2. P_ICRS is a free vector, merely a direction, typically of unit
   magnitude, and not bound to any particular spatial origin, such
   as the Earth, Sun or SSB.  No assumptions are made about whether
   it represents starlight and embodies astrometric effects such as
   parallax or aberration.  The transformation is approximately that
   between mean J2000.0 right ascension and declination and ecliptic
   longitude and latitude, with only frame bias (always less than
   25 mas) to disturb this classical picture.

   3. The Vondrak et al. (2011, 2012) 400 millennia precession model
   agrees with the IAU 2006 precession at J2000.0 and stays within
   100 microarcseconds during the 20th and 21st centuries.  It is
   accurate to a few arcseconds throughout the historical period,
   worsening to a few tenths of a degree at the end of the
   +/- 200,000 year time span.

Called:
   iauLtpequ    equator pole, long term
   iauLtpecl    ecliptic pole, long term
   iauPxp       vector product
   iauPn        normalize vector

References:

   Vondrak, J., Capitaine, N. and Wallace, P., 2011, New precession
   expressions, valid for long time intervals, Astron.Astrophys. 534,
   A22

   Vondrak, J., Capitaine, N. and Wallace, P., 2012, New precession
   expressions, valid for long time intervals (Corrigendum),
   Astron.Astrophys. 541, C1

This revision:  2015 December 6

SOFA release 2018-01-30

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

# void iauLtecm(double epj, double rm[3][3])
function iauLtecm(epj::Real)
   rm = zeros(Float64, 3, 3)
   ccall((:iauLtecm, libsofa_c), Cvoid,
         (Cdouble, Ref{Cdouble}),
          convert(Float64, epj),
          rm)

   return SMatrix{3,3}(rm')
end


text:
using Test
using Spandex

@testset "add 1" begin
    local g = DirectedGraph{Int64}(3)
    g[1, 2] = 5
    g[1, 3] = 6
    g[2, 1] = 1
    g[2, 2] = 4
    g[2, 3] = 7
    g[3, 1] = 2
    g[3, 2] = 3
    g[3, 3] = 8
    local a = from_graph(g, 3, 3)

    clear!(g)
    g[1, 1] = -1
    g[1, 3] = 1
    g[2, 1] = -1
    g[2, 3] = 1
    g[3, 1] = -1
    g[3, 3] = 1
    local b = from_graph(g, 3, 3)

    clear!(g)
    g[1, 1] = -1
    g[1, 2] = 5
    g[1, 3] = 7
    g[2, 1] = 0
    g[2, 2] = 4
    g[2, 3] = 8
    g[3, 1] = 1
    g[3, 2] = 3
    g[3, 3] = 9
    local e = from_graph(g, 3, 3)

    local c = add(a, b)

    @test !equals(a, c)
    @test !equals(b, c)
    @test equals(e, c)
end

@testset "mul 1" begin
    local g = DirectedGraph{Int64}(3)
    g[1, 1] = 6
    g[1, 2] = -4
    g[1, 3] = 7
    g[2, 1] = 4
    g[2, 2] = -3
    g[2, 3] = 0
    g[3, 1] = 1
    g[3, 2] = 2
    g[3, 3] = 5
    local a = from_graph(g, 3, 3)

    clear!(g)
    g[1, 1] = 6
    g[1, 2] = 3
    g[2, 1] = -1
    g[2, 2] = 4
    g[3, 1] = 8
    g[3, 2] = -3
    local b = from_graph(g, 3, 2)

    clear!(g)
    g[1, 1] = 96
    g[1, 2] = -19
    g[2, 1] = 27
    g[2, 2] = 0
    g[3, 1] = 44
    g[3, 2] = -4
    local e = from_graph(g, 3, 2)

    local c = mul(a, b)

    @test !equals(a, c)
    @test !equals(b, c)
    @test equals(e, c)
end

@testset "mul 2" begin
    local g = DirectedGraph{Int64}(3)
    g[1, 1] = 1
    g[1, 3] = 3
    g[2, 2] = 5
    g[3, 2] = 8
    local a = from_graph(g, 3, 3)

    clear!(g)
    g[1, 2] = 2
    g[1, 3] = 3
    g[2, 1] = 4
    g[2, 3] = 6
    g[3, 1] = 7
    g[3, 2] = 8
    local b = from_graph(g, 3, 3)

    clear!(g)
    g[1, 1] = 21
    g[1, 2] = 26
    g[1, 3] = 3
    g[2, 1] = 20
    g[2, 3] = 30
    g[3, 1] = 32
    g[3, 3] = 48
    local e = from_graph(g, 3, 3)

    local c = mul(a, b)

    @test !equals(a, c)
    @test !equals(b, c)
    @test equals(e, c)
end

@testset "mul 3" begin
    local g = DirectedGraph{Int64}(3)
    g[1, 1] = 1
    g[1, 3] = 3
    g[2, 2] = 5
    g[3, 2] = 8
    local a = from_graph(g, 3, 3)

    clear!(g)
    g[1, 1] = 6
    g[2, 1] = -1
    g[3, 1] = 8
    local b = from_graph(g, 3, 1)

    clear!(g)
    g[1, 1] = 30
    g[2, 1] = -5
    g[3, 1] = -8
    local e = from_graph(g, 3, 1)

    local c = mul(a, b)

    @test !equals(a, c)
    @test !equals(b, c)
    @test equals(e, c)
end

@testset "mul 4" begin
    local g = DirectedGraph{Int64}(3)
    g[1, 1] = 1
    g[1, 3] = 3
    g[2, 2] = 5
    g[3, 2] = 8
    local a = from_graph(g, 3, 3)

    local x = [0, 3, 2]
    local z = mul(a, x)
    @test 3 == length(z)
    @test 6 == z[1]
    @test 15 == z[2]
    @test 24 == z[3]

    local y = SparseArray{Int64}(3)
    y[2] = 3
    y[3] = 2
    z = mul(a, y)
    @test 3 == z.size
    @test 3 == z.nnz
    @test 6 == z[1]
    @test 15 == z[2]
    @test 24 == z[3]

    y = from_csr(3, 1, [1, 1, 2, 3], [1, 1], [3, 2])
    z = mul(a, y)
    @test 3 == z.nnz
    @test 6 == get_rowwise(z, 1, 1)
    @test 15 == get_rowwise(z, 2, 1)
    @test 24 == get_rowwise(z, 3, 1)
end

@testset "mul 5" begin
    local g = DirectedGraph{Int64}(4)
    g[1, 1] = 1
    g[1, 3] = 3
    g[2, 2] = 5
    g[3, 2] = 8
    g[4, 1] = 10
    g[4, 2] = 11
    g[4, 3] = 12
    local a = from_graph(g, 4, 3)

    local b1 = [0, 2, 0]
    local c1 = mul(a, b1)
    @test 4 == length(c1)
    @test 0 == c1[1]
    @test 10 == c1[2]
    @test 16 == c1[3]
    @test 22 == c1[4]

    local b2 = SparseArray{Int64}(3)
    b2[2] = 2
    local c2 = mul(a, b2)
    @test 4 == c2.size
    @test 3 == c2.nnz
    @test 10 == c2[2]
    @test 16 == c2[3]
    @test 22 == c2[4]

    local b3 = from_csr(3, 1, [1, 1, 2, 2], [1], [2])
    local c3 = mul(a, b3)
    @test 3 == c3.nnz
    @test 10 == get_rowwise(c3, 2, 1)
    @test 16 == get_rowwise(c3, 3, 1)
    @test 22 == get_rowwise(c3, 4, 1)
end

@testset "sqr 1" begin
    local g = DirectedGraph{Int64}(3)
    g[1, 2] = 1
    g[2, 1] = 2
    g[3, 2] = 3
    local a = from_graph(g, 3, 2)

    local ata = sqr_sym(a)

    @test 2 == ata.nnz
    @test contains(ata, 1, 1)
    @test contains(ata, 2, 2)
end

@testset "sqr 2" begin
    local g = DirectedGraph{Float64}(4)
    g[2, 5] = 0.02675
    g[3, 1] = 0.78664
    g[4, 1] = 0.26856
    g[4, 4] = 0.51423
    g[4, 6] = 0.46234
    local a = from_graph(g, 4, 6)

    resize!(g, 6)
    g[1, 1] = 0.69094
    g[4, 1] = 0.1381
    g[4, 4] = 0.26444
    g[5, 5] = 0.00072
    g[6, 1] = 0.12417
    g[6, 4] = 0.23775
    g[6, 6] = 0.21376
    local e = from_graph(g, 6, 6)

    local ata = sqr_sym(a)

    @test 12 == ata.nnz
    @test contains(ata, 1, 1)
    @test contains(ata, 2, 2)
    @test contains(ata, 3, 3)
    @test contains(ata, 4, 1)
    @test contains(ata, 4, 4)
    @test contains(ata, 5, 5)
    @test contains(ata, 6, 1)
    @test contains(ata, 6, 4)
    @test contains(ata, 6, 6)
    @test contains(ata, 1, 4)
    @test contains(ata, 1, 6)
    @test contains(ata, 4, 6)
end

@testset "sqr 3" begin
    local g = DirectedGraph{Int64}(3)
    g[1, 2] = 1
    g[2, 1] = 2
    g[3, 2] = 3
    local a = from_graph(g, 3, 2)

    clear!(g)
    g[1, 1] = 4
    g[2, 2] = 10
    local e = from_graph(g, 2, 2)

    local ata = sqr(a)

    @test equals(ata, e)

    local e2 = mul(transpose(a), a)

    @test equals(ata, e2)
end

@testset "sqr 4" begin
    local g = DirectedGraph{Float64}(4)
    g[2, 5] = 0.02674817948
    g[3, 1] = 0.7866442604
    g[4, 1] = 0.2685635172
    g[4, 4] = 0.5142332896
    g[4, 6] = 0.462339404
    local a = from_graph(g, 4, 6)

    resize!(g, 6)
    g[1, 1] = 0.6909355552
    g[2, 2] = 0.0
    g[3, 3] = 0.0
    g[4, 1] = 0.13810430091
    g[4, 4] = 0.2644358761
    g[5, 5] = 0.0007154651053
    g[6, 1] = 0.1241674965
    g[6, 4] = 0.2377503126
    g[6, 6] = 0.2137577245
    g[1, 4] = 0.13810430091
    g[1, 6] = 0.1241674965
    g[4, 6] = 0.2377503126
    local e = from_graph(g, 6, 6)

    local ata = sqr(a)

    ata.values = round.(ata.values, digits = 8)
    e.values = round.(ata.values, digits = 8)

    @test equals(ata, e)
end


text:
module Signals

using ..Network
export take_due_signals, connect, Signal, add_acceptor, add_donor, can_add_acceptor, can_add_donor,
    name_back_spike, BackSpike, name_delta_v, DeltaV, name_dc_instruction, DCInstruction,
    name_new_dc, NewDC, ForwardSignal, BackwardSignal

abstract type ForwardSignal <: Signal end
abstract type BackwardSignal <: Signal end
function amplify end

# struct MarkramDelta{T <: AbstractFloat} <: ForwardSignal
#     delta::T
# end

# struct ExctDeltaCond{T <: AbstractFloat} <: ForwardSignal
#     delta_cond::T
# end

# struct InhbtDeltaCond{T <: AbstractFloat} <: ForwardSignal
#     delta_cond::T
# end

# struct DeltaCond{T <: AbstractFloat} <: ForwardSignal
#     delta_cond::T
# end

const name_delta_v = "DeltaV"
struct DeltaV{T <: AbstractFloat} <: ForwardSignal
    delta_v::T
end

amplify(s::DeltaV{T}, w::T) where {T <: AbstractFloat} = DeltaV(s.delta_v * w)

## state_dict is for recording states, not for network-serialization. temporarily commented.
#state_dict(s::DeltaV{T}) where {T <: AbstractFloat} = Dict(["delta_v" => s.delta_v])

const name_dc_instruction = "DCInstruction"
struct DCInstruction{T <: AbstractFloat} <: ForwardSignal
    previous::T
    new::T
end

amplify(s::DCInstruction{T}, w::T) where {T <: AbstractFloat} = DCInstruction(s.previous * w, s.new * w)

const name_new_dc = "NewDC"
struct NewDC{T <: AbstractFloat} <: ForwardSignal
    current::T
end

const name_back_spike = "BackSpike"
struct BackSpike <: BackwardSignal
end

## state_dict is for recording states, not for network-serialization. temporarily commented.
#state_dict(s::DCInstruction{T}) where {T <: AbstractFloat} = Dict(["precious" => s.previous, "new" => s.new])

function can_add_acceptor end # (agent, signal_name) -> bool
function can_add_donor end
function add_acceptor end # (agent, signal_name, acptr_address) -> ()
function add_donor end # (agent, signal_name, dnr_address) -> ()

function connect(network::Dict{String, Population{T}},
                 signal_name::String,
                 donor_address::Address,
                 acceptor_address::Address) where {T <: Unsigned}
    dnr = get_agent(network, donor_address)
    acptr = get_agent(network, acceptor_address)
    if can_add_acceptor(dnr, signal_name) && can_add_donor(acptr, signal_name)
        add_acceptor(dnr, signal_name, acceptor_address)
        add_donor(acptr, signal_name, donor_address)
    else
        error(
            "connect failed. Donor: $donor_address; acceptor: $acceptor_address; signal_name: $signal_name"
        )
    end
end

end # Module end


text:
using Sundials

## f routine. Compute function f(t,y).

function f(t, y)
    ydot[1] = -0.04*y[1] + 1.0e4*y[2]*y[3]
    ydot[3] = 3.0e7*y[2]*y[2]
    ydot[2] = -ydot[1] - ydot[3]
    return ydot
end

## g routine. Event function whose roots cause integration to stop.

function g(t, y)
    gout[1] = y[1] - 0.0001
    gout[2] = y[3] - 0.01
    return(gout, [false,false], [0,0])
end

## Jacobian routine. Compute J(t,y) = df/dy.

# Using StrPack
# -- it works
# -- it's clunky and dependent on Sundials version
# Note: this is for Sundials v 2.4; the structure changed for v 2.5
# Because of this, I'm commenting this out.
## using StrPack
## @struct type J_DlsMat
##     typ::Int32
##     M::Int32
##     N::Int32
##     ldim::Int32
##     mu::Int32
##     ml::Int32
##     s_mu::Int32
##     data::Ptr{Float64}
##     ldata::Int32
##     cols::Ptr{Ptr{Float64}}
## end
# Note: Here is the (untested) structure for v 2.5
## using StrPack
## @struct type J_DlsMat
##     typ::Int32
##     M::Int
##     N::Int
##     ldim::Int
##     mu::Int
##     ml::Int
##     s_mu::Int
##     data::Ptr{Float64}
##     ldata::Int
##     cols::Ptr{Ptr{Float64}}
## end

# The following works if the code above is uncommented.
function Jac(N, t, y, fy, Jptr, user_data,
             tmp1, tmp2, tmp3)
    y = Sundials.asarray(y)
    dlsmat = unpack(IOString(pointer_to_array(convert(Ptr{Uint8}, Jptr),
                                              (sum(map(sizeof, J_DlsMat.types))+10,))),
                    J_DlsMat)
    J = pointer_to_array(unsafe_ref(dlsmat.cols), (int(neq), int(neq)), false)
    J[1,1] = -0.04
    J[1,2] = 1.0e4*y[3]
    J[1,3] = 1.0e4*y[2]
    J[2,1] = 0.04
    J[2,2] = -1.0e4*y[3] - 6.0e7*y[2]
    J[2,3] = -1.0e4*y[2]
    J[3,2] = 6.0e7*y[2]
    return int32(0)
end

y = [1.0,0.0,0.0]
reltol = 1e-4
abstol = [1e-8, 1e-14, 1e-6]

#cvode_mem = Sundials.CVodeCreate(Sundials.CV_BDF, Sundials.CV_NEWTON)
#flag = Sundials.CVodeInit(cvode_mem, f, t0, y)
#flag = Sundials.CVodeSVtolerances(cvode_mem, reltol, abstol)
#flag = Sundials.CVodeRootInit(cvode_mem, int32(2), g)
#flag = Sundials.CVDense(cvode_mem, neq)
## flag = Sundials.CVDlsSetDenseJacFn(cvode_mem, Jac)  # works, but clunky, see above

t = [0, 4*10 .^(-1.:10.)]
sol = Sundials.cvodej(f, t, y, events=g, reltol=reltol,abstol=abstol)


text:
using Printf
using JLD2
using Oceananigans.Utils

"""
    JLD2OutputWriter{I, T, O, IF, IN, KW} <: AbstractOutputWriter

An output writer for writing to JLD2 files.
"""
mutable struct JLD2OutputWriter{I, T, O, FS, D, IF, IN, KW} <: AbstractOutputWriter
              filepath :: String
               outputs :: O
    iteration_interval :: I
         time_interval :: T
          field_slicer :: FS
            array_type :: D
                  init :: IF
             including :: IN
              previous :: Float64
                  part :: Int
          max_filesize :: Float64
                 force :: Bool
               verbose :: Bool
               jld2_kw :: KW
end

noinit(args...) = nothing

"""
    JLD2OutputWriter(model, outputs; prefix,
                                                       dir = ".",
                                        iteration_interval = nothing,
                                             time_interval = nothing,
                                     time_averaging_window = nothing,
                                     time_averaging_stride = 1,
                                              max_filesize = Inf,
                                                     force = false,
                                                      init = noinit,
                                                 including = [:grid, :coriolis, :buoyancy, :closure],
                                                   verbose = false,
                                                      part = 1,
                                                array_type = Array{Float32},
                                                   jld2_kw = Dict{Symbol, Any}())

Construct a `JLD2OutputWriter` for an Oceananigans `model` that writes `label, output` pairs
in `outputs` to a JLD2 file.

The argument `outputs` may be a `Dict` or `NamedTuple`. The keys of `outputs` are symbols or
strings that "name" output data. The values of `outputs` are either `AbstractField`s, objects that
are called with the signature `output(model)`, or `WindowedTimeAverage`s of `AbstractFields`s,
functions, or callable objects.

Keyword arguments
=================

    ## Filenaming

    - `prefix`: Descriptive filename prefixed to all output files.

    - `dir`: Directory to save output to.
             Default: "." (current working directory).

    ## Output frequency and time-averaging

    - `iteration_interval`: Save output every `iteration_interval` model iterations.

    - `time_interval`: Save output every `time_interval` units of `model.clock.time`.

    - `time_averaging_window`: Specifies a time window over which each member of `output` is averaged before
                               being saved. For this each member of output is converted to
                               `Oceananigans.Diagnostics.WindowedTimeAverage`.
                               Default `nothing` indicates no averaging.

    - `time_averaging_stride`: Specifies a iteration 'stride' between the calculation of each `output` during
                               time-averaging. Longer strides means that output is calculated less frequently,
                               and that the resulting time-average is faster to compute, but less accurate.
                               Default: 1.

    ## Slicing and type conversion prior to output

    - `field_slicer`: An object for slicing field output in ``(x, y, z)``, including omitting halos.
                      Has no effect on output that is not a field. `field_slicer = nothing` means
                      no slicing occurs, so that all field data, including halo regions, is saved.
                      Default: FieldSlicer(), which slices halo regions.

    - `array_type`: The array type to which output arrays are converted to prior to saving.
                    Default: Array{Float32}.

    ## File management

    - `max_filesize`: The writer will stop writing to the output file once the file size exceeds `max_filesize`,
                      and write to a new one with a consistent naming scheme ending in `part1`, `part2`, etc.
                      Defaults to `Inf`.

    - `force`: Remove existing files if their filenames conflict.
               Default: `false`.

    ## Output file metadata management

    - `init`: A function of the form `init(file, model)` that runs when a JLD2 output file is initialized.
              Default: `noinit(args...) = nothing`.

    - `including`: List of model properties to save with every file.
                   Default: `[:grid, :coriolis, :buoyancy, :closure]`

    ## Miscellaneous keywords

    - `verbose`: Log what the output writer is doing with statistics on compute/write times and file sizes.
                 Default: `false`.

    - `part`: The starting part number used if `max_filesize` is finite.
              Default: 1.

    - `jld2_kw`: Dict of kwargs to be passed to `jldopen` when data is written.
"""
function JLD2OutputWriter(model, outputs; prefix,
                                                            dir = ".",
                                             iteration_interval = nothing,
                                                  time_interval = nothing,
                                          time_averaging_window = nothing,
                                          time_averaging_stride = 1,
                                                   field_slicer = FieldSlicer(),
                                                     array_type = Array{Float32},
                                                   max_filesize = Inf,
                                                          force = false,
                                                           init = noinit,
                                                      including = [:grid, :coriolis, :buoyancy, :closure],
                                                        verbose = false,
                                                           part = 1,
                                                        jld2_kw = Dict{Symbol, Any}())

    validate_intervals(iteration_interval, time_interval)

    # Convert each output to WindowedTimeAverage if time_averaging_window is specified
    if !isnothing(time_averaging_window)

        !isnothing(iteration_interval) && error("Cannot specify iteration_interval with time_averaging_window.")

        output_names = Tuple(keys(outputs))

        averaged_output =
            Tuple(
                  WindowedTimeAverage(outputs[name], model; time_interval = time_interval,
                                                              time_window = time_averaging_window,
                                                                   stride = time_averaging_stride,
                                                             field_slicer = field_slicer)
                  for name in output_names
                 )

        outputs = NamedTuple{output_names}(averaged_output)
    end

    mkpath(dir)
    filepath = joinpath(dir, prefix * ".jld2")
    force && isfile(filepath) && rm(filepath, force=true)

    jldopen(filepath, "a+"; jld2_kw...) do file
        init(file, model)
        saveproperties!(file, model, including)
    end

    return JLD2OutputWriter(filepath,
                            outputs,
                            iteration_interval,
                            time_interval,
                            field_slicer,
                            array_type,
                            init,
                            including,
                            0.0,
                            part,
                            max_filesize,
                            force,
                            verbose,
                            jld2_kw)
end

function write_output!(writer::JLD2OutputWriter, model)

    verbose = writer.verbose

    # Fetch JLD2 output and store in dictionary `data`
    verbose && @info @sprintf("Fetching JLD2 output %s...", keys(writer.outputs))

    tc = Base.@elapsed data = Dict((name, fetch_and_convert_output(output, model, writer)) for (name, output)
                                   in zip(keys(writer.outputs), values(writer.outputs)))

    verbose && @info "Fetching time: $(prettytime(tc))"

    # Start a new file if the filesize exceeds max_filesize
    filesize(writer.filepath) >= writer.max_filesize && start_next_file(model, writer)

    # Write output from `data`
    path = writer.filepath
    verbose && @info "Writing JLD2 output $(keys(writer.outputs)) to $path..."

    start_time, old_filesize = time_ns(), filesize(path)

    jld2output!(path, model.clock.iteration, model.clock.time, data, writer.jld2_kw)

    end_time, new_filesize = time_ns(), filesize(path)

    verbose && @info @sprintf("Writing done: time=%s, size=%s, Δsize=%s",
                              prettytime((start_time - end_time) / 1e9),
                              pretty_filesize(new_filesize),
                              pretty_filesize(new_filesize - old_filesize))

    return nothing
end

"""
    jld2output!(path, iter, time, data, kwargs)

Write the (name, value) pairs in `data`, including the simulation
`time`, to the JLD2 file at `path` in the `timeseries` group,
stamping them with `iter` and using `kwargs` when opening
the JLD2 file.
"""
function jld2output!(path, iter, time, data, kwargs)
    jldopen(path, "r+"; kwargs...) do file
        file["timeseries/t/$iter"] = time
        for (name, datum) in data
            file["timeseries/$name/$iter"] = datum
        end
    end
    return nothing
end

function start_next_file(model, writer::JLD2OutputWriter)
    verbose = writer.verbose
    sz = filesize(writer.filepath)
    verbose && @info begin
        "Filesize $(pretty_filesize(sz)) has exceeded maximum file size $(pretty_filesize(writer.max_filesize))."
    end

    if writer.part == 1
        part1_path = replace(writer.filepath, r".jld2$" => "_part1.jld2")
        verbose && @info "Renaming first part: $(writer.filepath) -> $part1_path"
        mv(writer.filepath, part1_path, force=writer.force)
        writer.filepath = part1_path
    end

    writer.part += 1
    writer.filepath = replace(writer.filepath, r"part\d+.jld2$" => "part" * string(writer.part) * ".jld2")
    writer.force && isfile(writer.filepath) && rm(writer.filepath, force=true)
    verbose && @info "Now writing to: $(writer.filepath)"

    jldopen(writer.filepath, "a+"; writer.jld2_kw...) do file
        writer.init(file, model)
        saveproperties!(file, model, writer.including)
    end

    return nothing
end


text:
function ising_local_MPO(::Type{T}, shift::Number=0.0; hz::Number=1.0, hx::Number=0.0)::MPO_open{T} where {T}
    E = Matrix{Float64}(I,2,2)
    X = [0.0 1.0; 1.0 0.0]
    Z = [1.0 0.0; 0.0 -1.0]

    h1_11 = -hz*Z - hx*X + shift*I

    h1 = zeros(eltype(h1_11), 2,2,1,2)
    h2 = zeros(eltype(h1_11), 2,2,2,1)

    h1[:,:,1,1] = h1_11
    h1[:,:,1,2] = -X

    h2[:,:,1,1] = E
    h2[:,:,2,1] = X

    h1 = permutedims(h1, (3,2,4,1)) #[m1,ket,m2,bra]
    h2 = permutedims(h2, (3,2,4,1)) #[m1,ket,m2,bra]

    h1 = convert(MPOTensor{T}, h1)
    h2 = convert(MPOTensor{T}, h2)

    MPOTensor{T}[h1, h2]
end

function ising_PBC_MPO(::Type{T}; hz::Real=1.0, hx::Real=0.0)::MPO_PBC_uniform{T} where {T}
    E = Matrix{Float64}(I,2,2)
    X = [0.0 1.0; 1.0 0.0]
    Z = [1.0 0.0; 0.0 -1.0]

    hM = zeros(Float64, 2,2,3,3)
    hB = zeros(Float64, 2,2,3,3)

    hM[:,:,1,1] = E
    hM[:,:,2,1] = X
    hM[:,:,3,1] = -hz*Z - hx*X
    hM[:,:,3,2] = -X
    hM[:,:,3,3] = E

    hB[:,:,1,:] = hM[:,:,3,:]
    hB[:,:,2,3] = X

    hM = permutedims(hM, (3,2,4,1)) #[m1,ket,m2,bra]
    hB = permutedims(hB, (3,2,4,1)) #[m1,ket,m2,bra]

    hM = convert(MPOTensor{T}, hM)
    hB = convert(MPOTensor{T}, hB)

    (hB, hM)
end

function ising_OBC_MPO(::Type{T}; hz::Real=1.0, hx::Real=0.0)::MPO_open_uniform{T} where {T}
    E = Matrix{Float64}(I,2,2)
    X = [0.0 1.0; 1.0 0.0]
    Z = [1.0 0.0; 0.0 -1.0]

    hM = zeros(Float64, 2,2,3,3)

    hM[:,:,1,1] = E
    hM[:,:,2,1] = X
    hM[:,:,3,1] = -hz*Z - hx*X
    hM[:,:,3,2] = -X
    hM[:,:,3,3] = E

    hL = zeros(Float64, 2,2,1,3)
    hR = zeros(Float64, 2,2,3,1)

    hL[:,:,1,:] = hM[:,:,3,:]
    hR[:,:,:,1] = hM[:,:,:,1]

    hM = permutedims(hM, (3,2,4,1)) #[m1,ket,m2,bra]
    hL = permutedims(hL, (3,2,4,1)) #[m1,ket,m2,bra]
    hR = permutedims(hR, (3,2,4,1)) #[m1,ket,m2,bra]

    hM = convert(MPOTensor{T}, hM)
    hL = convert(MPOTensor{T}, hL)
    hR = convert(MPOTensor{T}, hR)

    (hL, hM, hR)
end

function ising_PBC_MPO_split(::Type{T}; hz::Real=1.0, hx::Real=0.0)::MPO_PBC_uniform_split{T} where {T}
    X = [0.0 1.0; 1.0 0.0]
    hL = reshape(-X, (2,2,1,1))
    hR = reshape(X, (2,2,1,1))
    hL = convert(MPOTensor{T}, permutedims(hL, (3,2,4,1)))
    hR = convert(MPOTensor{T}, permutedims(hR, (3,2,4,1)))

    h_B = MPOTensor{T}[hL, hR]

    (ising_OBC_MPO(T, hz=hz, hx=hx), h_B)
end

function ising_Hn_MPO_split(::Type{T}, n::Integer, N::Integer; hz::Real=1.0, hx::Real=0.0) where {T}
    (hL, hM, hR), (hb1, hb2) = ising_PBC_MPO_split(T, hz=hz, hx=hx)

    hL[1,:,1,:] *= cis(n*2π/N)
    hL[1,:,2,:] *= cis(n*1.5*2π/N)

    hR[3,:,1,:] *= cis(n*N*2π/N) #I realise this is not strictly necessary, but it shows intent! :)

    get_hM = (j::Integer)->begin
        hM_j = copy(hM)
        hM_j[3,:,1,:] *= cis(n*j*2π/N)
        hM_j[3,:,2,:] *= cis(n*(j+0.5)*2π/N)
        hM_j
    end

    Hn_OBC = MPOTensor{T}[hL, (get_hM(n) for n in 2:N-1)..., hR]

    hb1 *= cis(n*(N+0.5)*2π/N)

    Hn_b = MPOTensor{T}[hb1, hb2]

    n, (Hn_OBC, Hn_b)
end

#------------------------------
function heis_local_MPO(::Type{T}; J::Number=1, spin::Real=1//2)::MPO_open{T} where {T}
    if spin == 1//2
        E = Matrix{ComplexF64}(I,2,2)
        X = complex([0.0 1.0; 1.0 0.0])
        Y = 1.0im*[0.0 -1.0; 1.0 0.0]
        Z = complex([1.0 0.0; 0.0 -1.0])
        zer = zero(E)
    else
        ValueError()
    end

    h1 = J*[X Y Z]
    h2 = [X;
          Y;
          Z]

    h1 = reshape(h1, (2,1,2,3))
    h2 = reshape(h2, (2,3,2,1))

    h1 = permutedims(h1, (2,3,4,1))
    h2 = permutedims(h2, (2,3,4,1))

    MPOTensor{T}[convert(MPOTensor{T}, h1), convert(MPOTensor{T}, h2)]
end

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

"""
-(lamda*X1*X2 + delta1*X1*X3 + delta2*Z1*Z2 + hz*Z1)
"""
function ANNNI_local_MPO(::Type{T}; hz::Real=1.0, delta1::Real=0.0, delta2::Real=0.0, lambda::Real=1.0)::MPO_open{T} where {T}
    E = Matrix{Float64}(I,2,2)
    X = [0.0 1.0; 1.0 0.0]
    Z = [1.0 0.0; 0.0 -1.0]

    h1 = zeros(Float64, 2,2,1,2)
    h2 = zeros(Float64, 2,2,2,2)
    h3 = zeros(Float64, 2,2,2,1)

    h1[:,:,1,1] = -X
    h1[:,:,1,2] = -Z

    h2[:,:,1,1] = lambda*X
    h2[:,:,1,2] = E
    h2[:,:,2,1] = delta2 * Z + I*hz

    h3[:,:,1,1] = E
    h3[:,:,2,1] = delta1 * X

    h1 = permutedims(h1, (3,2,4,1)) #[m1,ket,m2,bra]
    h2 = permutedims(h2, (3,2,4,1)) #[m1,ket,m2,bra]
    h3 = permutedims(h3, (3,2,4,1)) #[m1,ket,m2,bra]

    h1 = convert(MPOTensor{T}, h1)
    h2 = convert(MPOTensor{T}, h2)
    h3 = convert(MPOTensor{T}, h3)

    MPOTensor{T}[h1, h2, h3]
end

function ANNNI_OBC_MPO(::Type{T}; hz::Real=1.0, delta1::Real=0.0, delta2::Real=0.0, lambda::Real=1.0)::MPO_open_uniform{T} where {T}
    E = Matrix{Float64}(I,2,2)
    X = [0.0 1.0; 1.0 0.0]
    Z = [1.0 0.0; 0.0 -1.0]

    hM = zeros(Float64, 2,2,5,5)

    hM[:,:,1,1] = E
    hM[:,:,2,1] = X
    hM[:,:,3,1] = Z
    hM[:,:,4,2] = E
    hM[:,:,5,2] = -lambda*X
    hM[:,:,5,3] = -delta2*Z-hz*I
    hM[:,:,5,4] = -delta1*X
    hM[:,:,5,5] = E

    hL = zeros(Float64, 2,2,1,5)
    hR = zeros(Float64, 2,2,5,1)

    hL[:,:,1,:] = hM[:,:,5,:]
    hL[:,:,1,1] = -hz*Z
    hR[:,:,:,1] = hM[:,:,:,1]

    hM = permutedims(hM, (3,2,4,1)) #[m1,ket,m2,bra]
    hL = permutedims(hL, (3,2,4,1)) #[m1,ket,m2,bra]
    hR = permutedims(hR, (3,2,4,1)) #[m1,ket,m2,bra]

    hM = convert(MPOTensor{T}, hM)
    hL = convert(MPOTensor{T}, hL)
    hR = convert(MPOTensor{T}, hR)

    (hL, hM, hR)
end

function ANNNI_PBC_MPO_split(::Type{T}; hz::Real=1.0, delta1::Real=0.0, delta2::Real=0.0, lambda::Real=1.0)::MPO_PBC_uniform_split{T} where {T}
    #Turn off hz for boundary terms, since they're already in the OBC part.
    hL, hM, hR = ANNNI_OBC_MPO(T, hz=0.0, delta1=delta1, delta2=delta2, lambda=lambda)

    #Pick out only the needed terms from the OBC Hamiltonian. Reduces the max. bond dimension to 3.
    hL = hL[:,:,4:5,:]
    hM1 = hM[4:5,:,2:4,:]
    hM2 = hM[2:4,:,1:2,:]
    hR = hR[1:2,:,:,:]

    h_B = MPOTensor{T}[hL, hM1, hM2, hR]

    (ANNNI_OBC_MPO(T, hz=hz, delta1=delta1, delta2=delta2, lambda=lambda), h_B)
end

function ANNNI_Hn_MPO_split(::Type{T}, n::Integer, N::Integer; hz::Real=1.0, delta1::Real=0.0, delta2::Real=0.0, lambda::Real=1.0) where {T}
    (hL, hM, hR), (hb1, hb2, hb3, hb4) = ANNNI_PBC_MPO_split(T, hz=hz, delta1=delta1, delta2=delta2, lambda=lambda)

    Z = [1.0 0.0; 0.0 -1.0]

    hL[1,:,1,:] *= cis(n*2π/N)
    hL[1,:,2,:] *= cis(n*1.5*2π/N)
    hL[1,:,3,:] = -delta2 * cis(n*(1.5)*2π/N) * Z - hz * cis(n*2*2π/N) * I
    hL[1,:,4,:] *= cis(n*2*2π/N)

    get_hM = (j::Integer)->begin
        hM_j = copy(hM)
        hM_j[5,:,2,:] *= cis(n*(j+0.5)*2π/N)
        hM_j[5,:,3,:] = -delta2 * cis(n*(j+0.5)*2π/N) * Z - hz * cis(n*(j+1)*2π/N) * I
        hM_j[5,:,4,:] *= cis(n*(j+1)*2π/N)
        hM_j
    end

    Hn_OBC = MPOTensor{T}[hL, (get_hM(n) for n in 2:N-1)..., hR]

    #Boundary MPO tensor for site N-1
    hb1[1,:,1,:] *= cis(n*N*2π/N)

    #Boundary tensor for site N
    hb2[2,:,1,:] *= cis(n*0.5*2π/N)
    hb2[2,:,2,:] = -delta2 * cis(n*0.5*2π/N) * Z
    hb2[2,:,3,:] *= cis(n*1*2π/N)

    Hn_b = MPOTensor{T}[hb1, hb2, hb3, hb4]

    n, (Hn_OBC, Hn_b)
end

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

function OBF_local_MPO(::Type{T}; lambda1::Real=1.0, lambda2::Real=1.0,
                            lambdaC::Real=0.0)::MPO_open{T} where {T}
    E = Matrix{Float64}(I,2,2)
    X = [0.0 1.0; 1.0 0.0]
    Z = [1.0 0.0; 0.0 -1.0]
    Y = [0.0im -1.0im; 1.0im 0.0im]

    Tin = (lambdaC == 0.0 ? Float64 : ComplexF64)
    d = (lambdaC == 0.0 ? 2 : 3)

    h1 = zeros(Tin, 2,2,1,3)
    h2 = zeros(Tin, 2,2,3,d)
    h3 = zeros(Tin, 2,2,d,1)

    h1[:,:,1,1] = lambda2*Z
    h1[:,:,1,2] = lambda2*X - lambda1*E
    h1[:,:,1,3] = E

    h2[:,:,1,1] = Z
    h2[:,:,2,2] = Z
    h2[:,:,3,1] = -lambda1*E

    h3[:,:,1,1] = X
    h3[:,:,2,1] = Z

    if lambdaC != 0.0
        h2[:,:,3,2] = lambdaC * Y
        h2[:,:,3,3] = -lambdaC * Z

        h3[:,:,3,1] = Y
    end

    h1 = permutedims(h1, (3,2,4,1)) #[m1,ket,m2,bra]
    h2 = permutedims(h2, (3,2,4,1)) #[m1,ket,m2,bra]
    h3 = permutedims(h3, (3,2,4,1)) #[m1,ket,m2,bra]

    h1 = convert(MPOTensor{T}, h1)
    h2 = convert(MPOTensor{T}, h2)
    h3 = convert(MPOTensor{T}, h3)

    MPOTensor{T}[h1, h2, h3]
end

function OBF_OBC_MPO(::Type{T}; lambda1::Real=1.0, lambda2::Real=1.0,
                        lambdaC::Real=0.0)::MPO_open_uniform{T} where {T}
    E = Matrix{Float64}(I,2,2)
    X = [0.0 1.0; 1.0 0.0]
    Z = [1.0 0.0; 0.0 -1.0]
    Y = [0.0im -1.0im; 1.0im 0.0im]

    Tin = (lambdaC == 0.0 ? Float64 : ComplexF64)
    hM = zeros(Tin, 2,2,6,6)

    hM[:,:,1,1] = E
    hM[:,:,2,1] = X
    hM[:,:,3,1] = Z
    hM[:,:,4,2] = lambda2*Z #Z
    hM[:,:,5,3] = Z
    hM[:,:,6,1] = -lambda1*X
    hM[:,:,6,3] = -lambda1*Z
    hM[:,:,6,4] = Z #lambda2*Z
    hM[:,:,6,5] = lambda2*X
    hM[:,:,6,6] = E

    if lambdaC != 0.0
        hM[:,:,4,1] = -lambdaC*Y
        hM[:,:,6,3] += lambdaC*Y
    end

    hL = zeros(Tin, 2,2,1,6)
    hR = zeros(Tin, 2,2,6,1)

    hL[:,:,1,:] = hM[:,:,6,:]
    hR[:,:,:,1] = hM[:,:,:,1]

    hM = permutedims(hM, (3,2,4,1)) #[m1,ket,m2,bra]
    hL = permutedims(hL, (3,2,4,1)) #[m1,ket,m2,bra]
    hR = permutedims(hR, (3,2,4,1)) #[m1,ket,m2,bra]

    hM = convert(MPOTensor{T}, hM)
    hL = convert(MPOTensor{T}, hL)
    hR = convert(MPOTensor{T}, hR)

    (hL, hM, hR)
end

function OBF_PBC_MPO_split(::Type{T}; lambda1::Real=1.0, lambda2::Real=1.0,
                                lambdaC::Real=0.0)::MPO_PBC_uniform_split{T} where {T}
    hL, hM, hR = OBF_OBC_MPO(T, lambda1=lambda1, lambda2=lambda2, lambdaC=lambdaC)

    #Pick out only the needed terms from the OBC Hamiltonian. Reduces the max. bond dimension to 3.
    hL = hL[:,:,4:6,:]
    hM1 = hM[4:6,:,2:5,:]
    hM2 = hM[2:5,:,1:3,:]
    hR = hR[1:3,:,:,:]

    h_B = MPOTensor{T}[hL, hM1, hM2, hR]

    (OBF_OBC_MPO(T, lambda1=lambda1, lambda2=lambda2, lambdaC=lambdaC), h_B)
end

function OBF_Hn_MPO_split(::Type{T}, n::Integer, N::Integer; lambda1::Real=1.0, lambda2::Real=1.0,
    ZZXoff::Number=5/4, XZZoff::Number=3/4) where {T}
    (hL, hM, hR), (hb1, hb2, hb3, hb4) = OBF_PBC_MPO_split(T, lambda1=lambda1, lambda2=lambda2)

    Z = [1.0 0.0; 0.0 -1.0]

    hL[1,:,1,:] *= cis(n*2π/N)
    hL[1,:,3,:] *= cis(n*1.5*2π/N)

    #hL[1,:,4,:] *= cis(n*(1 + ZZXoff)*2π/N) #ZZX
    hL[1,:,5,:] *= cis(n*(1 + XZZoff)*2π/N) #XZZ

    get_hM = (j::Integer)->begin
        hM_j = copy(hM)
        hM_j[6,:,1,:] *= cis(n*j*2π/N)
        hM_j[6,:,3,:] *= cis(n*(j+0.5)*2π/N)
        #hM_j[6,:,4,:] *= cis(n*(j+ZZXoff)*2π/N) #ZZX
        hM_j[4,:,2,:] *= cis(n*(j-1+ZZXoff)*2π/N) #ZZX
        hM_j[6,:,5,:] *= cis(n*(j+XZZoff)*2π/N) #XZZ
        hM_j
    end

    Hn_OBC = MPOTensor{T}[hL, (get_hM(n) for n in 2:N-1)..., hR]

    #Boundary tensor for site N-1
    #hb1[1,:,1,:] *= cis(n*(N-1 + ZZXoff)*2π/N) #ZXX
    hb1[1,:,2,:] *= cis(n*(N-1 + XZZoff)*2π/N) #XXZ (N-1)

    #Boundary tensor for site N
    hb2[3,:,2,:] *= cis(n*0.5*2π/N)
    #hb2[3,:,3,:] *= cis(n*ZZXoff*2π/N) #ZZX
    hb2[1,:,1,:] *= cis(n*(N-1 + ZZXoff)*2π/N) #ZZX (N-1)
    hb2[3,:,4,:] *= cis(n*XZZoff*2π/N) #XXZ (N)

    #Boundary tensor for site 1
    hb3[3,:,2,:] *= cis(n*(ZZXoff)*2π/N) #ZZX (N)

    Hn_b = MPOTensor{T}[hb1, hb2, hb3, hb4]

    n, (Hn_OBC, Hn_b)
end

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

function weylops(p::Integer)
    om = cis(2π / p)
    U = diagm(ComplexF64[om^j for j in 0:p-1])
    V = diagm(ones(p - 1), 1)
    V[end, 1] = 1
    U, V, om
end

function potts3_OBC_MPO(::Type{T}; h::Real=1.0) where {T}
    E = Matrix{ComplexF64}(I,3,3)
    U, V, om = weylops(3)

    hM = zeros(ComplexF64, 3,3,4,4)

    hM[:,:,1,1] = E
    hM[:,:,2,1] = U
    hM[:,:,3,1] = U'
    hM[:,:,4,1] = -h/2 * (V+V')
    hM[:,:,4,2] = -0.5*U'
    hM[:,:,4,3] = -0.5*U
    hM[:,:,4,4] = E

    hL = zeros(ComplexF64, 3,3,1,4)
    hR = zeros(ComplexF64, 3,3,4,1)

    hL[:,:,1,:] = hM[:,:,4,:]
    hR[:,:,:,1] = hM[:,:,:,1]

    hM = permutedims(hM, (3,2,4,1)) #[m1,ket,m2,bra]
    hL = permutedims(hL, (3,2,4,1)) #[m1,ket,m2,bra]
    hR = permutedims(hR, (3,2,4,1)) #[m1,ket,m2,bra]

    hM = convert(MPOTensor{T}, hM)
    hL = convert(MPOTensor{T}, hL)
    hR = convert(MPOTensor{T}, hR)

    (hL, hM, hR)
end

function potts3_local_MPO(::Type{T}; h::Real=1.0) where {T}
    hL, hM, hR = potts3_OBC_MPO(T, h=h)

    MPOTensor{T}[hL[:,:,1:3,:], hR[1:3,:,:,:]]
end

function potts3_PBC_MPO_split(::Type{T}; h::Real=1.0)::MPO_PBC_uniform_split{T} where {T}
    hL, hM, hR = potts3_OBC_MPO(T, h=h)

    h_B = MPOTensor{T}[hL[:,:,2:3,:], hR[2:3,:,:,:]]

    ((hL, hM, hR), h_B)
end

function potts3_Hn_MPO_split(::Type{T}, n::Integer, N::Integer; h::Real=1.0) where {T}
    (hL, hM, hR), (hb1, hb2) = potts3_PBC_MPO_split(T, h=h)

    hL[1,:,1,:] *= cis(n*2π/N)
    hL[1,:,2,:] *= cis(n*1.5*2π/N)
    hL[1,:,3,:] *= cis(n*1.5*2π/N)

    hR[4,:,1,:] *= cis(n*N*2π/N)

    get_hM = (j::Integer)->begin
        hM_j = copy(hM)
        hM_j[4,:,1,:] *= cis(n*j*2π/N)
        hM_j[4,:,2,:] *= cis(n*(j+0.5)*2π/N)
        hM_j[4,:,3,:] *= cis(n*(j+0.5)*2π/N)
        hM_j
    end

    Hn_OBC = MPOTensor{T}[hL, (get_hM(n) for n in 2:N-1)..., hR]

    hb1 *= cis(n*0.5*2π/N)

    Hn_b = MPOTensor{T}[hb1, hb2]

    n, (Hn_OBC, Hn_b)
end

max_stars_repo_name: FuTen/puMPS

max_issues_repo_name: FuTen/puMPS

max_forks_repo_name: FuTen/puMPS


text:
# Some basic functions to interact with the data.

function gencov(d::Integer)
  A = rand(d, d)
  eye(d) + (A + A')/2
end

function plotclouds(cloudA, cloudB)
  scatter(cloudA[1,:], cloudA[2,:], label="-1")
  scatter!(cloudB[1,:], cloudB[2,:], label="1")
end

function errorrate(classA::MvNormal,
                   classB::MvNormal,
                   w::AbstractVector, n::Integer)
  testA = rand(classA, n)
  testB = rand(classB, n)
  X = [[testA testB]' ones(2*n)]
  Y = [-1.*ones(n); ones(n)]
  sum(sign.(X*w).≠Y)/n
end

function drawborder(f)
  xs = linspace(xlims()..., 100)
  ys = linspace(ylims()..., 100)
  contour!(xs, ys, f, levels=[0], colorbar=false, ls=:dash)
end

function plotdualitygap(m::Integer,
                        μ::Real,
                        numstepsarray::AbstractVector{T}) where T<:Integer
  dualgaps = [m/μ^i for i=0:length(numstepsarray)]
  numstepssumsarray = reduce(
    (arr, n) -> [arr; arr[end]+n] ,[0], numstepsarray)
  dualgaps2 = repeat(dualgaps, inner=2)[1:end-1]
  numstepssumsarray2 = repeat(numstepssumsarray, inner=2)[2:end]
  plot(numstepssumsarray2, dualgaps2, yscale=:log10,
       xlabel="Newton iterations", ylabel="duality gap")
end


text:
#=
   Likelihhod is chi^2 comparing observations of the radial velocity of the star to a model,
   where the model computes the velocity of the star as the linear superposition of the Keplerian orbit
   induced by each planet, i.e., neglecting mutual planet-planet interactions
   Priors based on SAMSI reference priors from Ford & Gregory 2006
=#

include("utils_internal.jl")                 # Manually add code to support autodiff of mod2pi

module RvModelCircularIndepNoise

include("rv_model_circ_param.jl")       # Code to convert parameter vector to named parameters
include("rv_model_circ.jl")             # Code to compute RV model for circular orbits
include("stat_model.jl")                # Common code to compute target density, logtarget, logprior_planets
include("noise_model_indep.jl")      # Code to compute target density, loglikelihood, logprior_noise

end # module RvModelCircularIndepNoise

module RvModelCircularGPNoise

include("rv_model_circ_param.jl")       # Code to convert parameter vector to named parameters
include("rv_model_circ.jl")             # Code to compute RV model for circular orbits
include("stat_model.jl")                # Common code to compute target density, logtarget, logprior_planets
include("noise_model_gp.jl")             # Code to compute target density, loglikelihood, logprior_noise

end # module RvModelCircularGPNoise

module RvModelKeplerianIndepNoise

include("rv_model_keplerian_param.jl")  # Code to convert parameter vector to named parameters
include("rv_model_keplerian.jl")        # Code to compute RV model for Keplerian orbits
include("stat_model.jl")                # Common code to compute target density, logtarget, logprior_planets
include("noise_model_indep.jl")      # Code to compute target density, loglikelihood, logprior_noise

end # module RvModelKeplerianIndepNoise

module RvModelKeplerianGPNoise

include("rv_model_keplerian_param.jl")  # Code to convert parameter vector to named parameters
include("rv_model_keplerian.jl")        # Code to compute RV model for Keplerian orbits
include("stat_model.jl")                # Common code to compute target density, logtarget, logprior_planets
include("noise_model_gp.jl")             # Code to compute target density, loglikelihood, logprior_noise

end # module RvModelKeplerianGPNoise


text:
export FactorGaussian

include(joinpath(@__DIR__, "../prob_operations.jl"))
include(joinpath(@__DIR__, "../graph_operations.jl"))
include(joinpath(@__DIR__, "../util.jl"))

mutable struct FactorGaussian
    """
    Composite Gaussian factor node.

        f(y, x, a, u, z) = 𝓝(y | a*x + u, z)

    This is a Gaussian node, but slightly more general: it contains a
    multiplication node with a transition coefficient and an addition node with
    a control coefficient. By default, these are clamped such that a standard
    Gaussian is employed.

    Examples
    ≡≡≡≡≡≡≡≡≡≡

    # States
    julia> x_1 = VarGaussian(:x_1);
    julia> x_2 = VarGaussian(:x_2);

    # Define Gaussian factor for state transition
    julia> f = FactorGaussian(:f, out=:x_2, mean=:x_1)

    """

    # Identifiers of edges/nodes in factor graph
    id::Symbol
    variables::Dict{String, Union{Float64, Symbol}}
    time::Union{Integer,Nothing}

    # Reaction parameters
    threshold::Float64
    fired::Bool

    # Additional properties
    verbose::Bool

    function FactorGaussian(id::Symbol;
                            out::Union{Float64, Symbol}=0.0,
                            mean::Union{Float64, Symbol}=0.0,
                            transition::Union{Float64, Symbol}=1.0,
                            control::Union{Float64, Symbol}=0.0,
                            precision::Union{Float64, Symbol}=1.0,
                            threshold::Float64=0.0,
                            time::Union{Integer,Nothing}=0,
                            fired::Bool=false,
                            verbose::Bool=false)

        # Keep track of recognition distributions
        variables = Dict{String, Union{Float64, Symbol}}("out" => out,
                                                         "mean" => mean,
                                                         "transition" => transition,
                                                         "control" => control,
                                                         "precision" => precision)

        # Construct instance
        return new(id, variables, time, threshold, fired, verbose)
    end
end

function message(graph::MetaGraph, node::FactorGaussian, var_id::Symbol)
    "Compute a message to an edge."

    # Extract means and variances of beliefs
    Ex, Vx = meanvar(graph, node.variables["out"])
    Em, Vm = meanvar(graph, node.variables["mean"])
    Ea, Va = meanvar(graph, node.variables["transition"])
    Eu, Vu = meanvar(graph, node.variables["control"])
    Eγ, Vγ = meanvar(graph, node.variables["precision"])

    # Find which variable belongs to var_id
    var_name = key_from_value(node.variables, var_id)

    if var_name == "out"

        message = Normal(Ea*Em, sqrt(inv(Eγ)))

    elseif var_name == "mean"

        message = Normal(Ex*Ea / (Va + Ea^2), sqrt(inv(Eγ*(Va + Ea^2))))

    elseif var_name == "transition"

        message = Normal(Ex*Em / (Vm + Em^2), sqrt(inv(Eγ*(Vm + Em^2))))

    elseif var_name == "control"

        error("Exception: not implemented yet.")

    elseif var_name == "precision"

        inv_scale = (Vx + Ex^2 - 2*Ex*Ea*Em + (Em^2 + Vm)*(Ea^2 + Va))/2
        message = Gamma(3/2., inv(inv_scale))

    else
        throw("Exception: variable unknown.")
    end

    return message
end

function act!(graph::MetaGraph, node::FactorGaussian, var_id::Symbol)
    "Send out message to connected edge."

    # Compute message
    message_f2v = message(graph, node, var_id)

    # Find edge id
    edge_id = Edge(graph[var_id, :id], graph[node.id, :id])

    # Pass message
    set_prop!(graph, edge_id, :message_factor2var, message_f2v)
end

function react!(graph::MetaGraph, node::FactorGaussian)
    "React to incoming messages from edges."

    # Keep track of whether node fired
    node.fired = false

    # Set of variable ids
    var_ids = Set{Symbol}([x for x in values(node.variables) if isa(x, Symbol)])

    # Loop over variables to check for incoming messages
    for var_id in var_ids

        # Find edge id
        edge = Edge(graph[var_id, :id], graph[node.id, :id])

        # If change in FE is large enough, record variable as updated
        if get_prop(graph, edge, :Δfree_energy) >= node.threshold

            # Loop over other variables
            for other_var in setdiff(var_ids, [var_id])

                # Pass message to other variable
                act!(graph, node, other_var)
            end

            # Record that node has fired
            node.fired = true
        end
    end
end

max_stars_repo_name: wmkouw/schedulefreeMP

max_issues_repo_name: wmkouw/schedulefreeMP

max_forks_repo_name: wmkouw/schedulefreeMP


text:
module JungleHelperSwingCassetteBlock

using ..Ahorn, Maple

@mapdef Entity "JungleHelper/SwingCassetteBlock" SwingCassetteBlock(x::Integer, y::Integer, width::Integer=Maple.defaultBlockWidth, height::Integer=Maple.defaultBlockHeight, index::Integer=0, tempo::Number=1.0)

const colorNames = Dict{String, Int}(
    "Blue" => 0,
    "Rose" => 1,
    "Bright Sun" => 2,
    "Malachite" => 3
)

const placements = Ahorn.PlacementDict(
    "Swing Cassette Block ($index - $color) (Jungle Helper)" => Ahorn.EntityPlacement(
        SwingCassetteBlock,
        "rectangle",
        Dict{String, Any}(
            "index" => index,
        )
    ) for (color, index) in colorNames
)

Ahorn.editingOptions(entity::SwingCassetteBlock) = Dict{String, Any}(
    "index" => colorNames
)

Ahorn.minimumSize(entity::SwingCassetteBlock) = 16, 16
Ahorn.resizable(entity::SwingCassetteBlock) = true, true

Ahorn.selection(entity::SwingCassetteBlock) = Ahorn.getEntityRectangle(entity)

const colors = Dict{Int, Ahorn.colorTupleType}(
    1 => (240, 73, 190, 255) ./ 255,
	2 => (252, 220, 58, 255) ./ 255,
	3 => (56, 224, 78, 255) ./ 255,
)

const defaultColor = (73, 170, 240, 255) ./ 255
const borderMultiplier = (0.9, 0.9, 0.9, 1)

const frame = "objects/cassetteblock/solid"

function getCassetteBlockRectangles(room::Maple.Room)
    entities = filter(e -> e.name == "cassetteBlock" || e.name == "JungleHelper/SwingCassetteBlock", room.entities)
    rects = Dict{Int, Array{Ahorn.Rectangle, 1}}()

    for e in entities
        index = get(e.data, "index", 0)
        rectList = get!(rects, index) do
            Ahorn.Rectangle[]
        end

        push!(rectList, Ahorn.Rectangle(
            Int(get(e.data, "x", 0)),
            Int(get(e.data, "y", 0)),
            Int(get(e.data, "width", 8)),
            Int(get(e.data, "height", 8))
        ))
    end

    return rects
end

# Is there a casette block we should connect to at the offset?
function notAdjacent(entity::SwingCassetteBlock, ox, oy, rects)
    x, y = Ahorn.position(entity)
    rect = Ahorn.Rectangle(x + ox + 4, y + oy + 4, 1, 1)

    for r in rects
        if Ahorn.checkCollision(r, rect)
            return false
        end
    end

    return true
end

function drawCassetteBlock(ctx::Ahorn.Cairo.CairoContext, entity::SwingCassetteBlock, room::Maple.Room)
    cassetteBlockRectangles = getCassetteBlockRectangles(room)

    x, y = Ahorn.position(entity)

    width = Int(get(entity.data, "width", 32))
    height = Int(get(entity.data, "height", 32))

    tileWidth = ceil(Int, width / 8)
    tileHeight = ceil(Int, height / 8)

    index = Int(get(entity.data, "index", 0))
    color = get(colors, index, defaultColor)

    rect = Ahorn.Rectangle(x, y, width, height)
    rects = get(cassetteBlockRectangles, index, Ahorn.Rectangle[])

    if !(rect in rects)
        push!(rects, rect)
    end

    for x in 1:tileWidth, y in 1:tileHeight
        drawX, drawY = (x - 1) * 8, (y - 1) * 8

        closedLeft = !notAdjacent(entity, drawX - 8, drawY, rects)
        closedRight = !notAdjacent(entity, drawX + 8, drawY, rects)
        closedUp = !notAdjacent(entity, drawX, drawY - 8, rects)
        closedDown = !notAdjacent(entity, drawX, drawY + 8, rects)
        completelyClosed = closedLeft && closedRight && closedUp && closedDown

        if completelyClosed
            if notAdjacent(entity, drawX + 8, drawY - 8, rects)
                Ahorn.drawImage(ctx, frame, drawX, drawY, 24, 0, 8, 8, tint=color)

            elseif notAdjacent(entity, drawX - 8, drawY - 8, rects)
                Ahorn.drawImage(ctx, frame, drawX, drawY, 24, 8, 8, 8, tint=color)

            elseif notAdjacent(entity, drawX + 8, drawY + 8, rects)
                Ahorn.drawImage(ctx, frame, drawX, drawY, 24, 16, 8, 8, tint=color)

            elseif notAdjacent(entity, drawX - 8, drawY + 8, rects)
                Ahorn.drawImage(ctx, frame, drawX, drawY, 24, 24, 8, 8, tint=color)

            else
                Ahorn.drawImage(ctx, frame, drawX, drawY, 8, 8, 8, 8, tint=color)
            end

        else
            if closedLeft && closedRight && !closedUp && closedDown
                Ahorn.drawImage(ctx, frame, drawX, drawY, 8, 0, 8, 8, tint=color)

            elseif closedLeft && closedRight && closedUp && !closedDown
                Ahorn.drawImage(ctx, frame, drawX, drawY, 8, 16, 8, 8, tint=color)

            elseif closedLeft && !closedRight && closedUp && closedDown
                Ahorn.drawImage(ctx, frame, drawX, drawY, 16, 8, 8, 8, tint=color)

            elseif !closedLeft && closedRight && closedUp && closedDown
                Ahorn.drawImage(ctx, frame, drawX, drawY, 0, 8, 8, 8, tint=color)

            elseif closedLeft && !closedRight && !closedUp && closedDown
                Ahorn.drawImage(ctx, frame, drawX, drawY, 16, 0, 8, 8, tint=color)

            elseif !closedLeft && closedRight && !closedUp && closedDown
                Ahorn.drawImage(ctx, frame, drawX, drawY, 0, 0, 8, 8, tint=color)

            elseif !closedLeft && closedRight && closedUp && !closedDown
                Ahorn.drawImage(ctx, frame, drawX, drawY, 0, 16, 8, 8, tint=color)

            elseif closedLeft && !closedRight && closedUp && !closedDown
                Ahorn.drawImage(ctx, frame, drawX, drawY, 16, 16, 8, 8, tint=color)
            end
        end
    end
end

Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::SwingCassetteBlock, room::Maple.Room) = drawCassetteBlock(ctx, entity, room)

end

max_stars_repo_name: Jackal-Celeste/JungleHelper

max_issues_repo_name: Jackal-Celeste/JungleHelper

max_forks_repo_name: Jackal-Celeste/JungleHelper


text:
using Luxor

struct TableCell
    content::Drawable
    valign::Symbol
    halign::Symbol
end

function TableCell(
    content;
    valign=nothing,
    halign=nothing,
    )
    valign = isnothing(valign) ? :middle : valign
    halign = isnothing(halign) ? :center : halign
    TableCell(
        content,
        valign,
        halign,
    )
end

tablecell(::Nothing; kwargs...) = nothing
tablecell(c::TableCell; kwargs...) = c # @todo kwargs needs to override here
tablecell(c; kwargs...) = TableCell(c; kwargs...)
getwidth(c::TableCell) = getwidth(c.content)
getheight(c::TableCell) = getheight(c.content)

function draw(c::TableCell, pos, cell_size)
    gsave()
    cell_height = cell_size[2]
    if c.valign == :top
        translate_y = cell_height÷2 - getheight(c)÷2
        translate(pos - (0, translate_y))
    else
        translate(pos)
    end
    draw(c.content)
    grestore()
end


text:
# query.jl

module QueryTests

using Test
using Dates
using QDates

# Name functions
mtk = QDates.QDate(2017,1,1) #先勝
ksrg = QDates.QDate(2017,2,1) #友引
yyi = QDates.QDate(2017,3,1) #先負
udk = QDates.QDate(2017,4,1) #仏滅
stk = QDates.QDate(2017,5,1) #大安
lstk = QDates.QDate(2017,5,true,1) #大安
mndk = QDates.QDate(2017,6,1) #赤口
fdk = QDates.QDate(2017,7,1) #先勝
hdk = QDates.QDate(2017,8,1) #友引
ngtk = QDates.QDate(2017,9,1) #先負
kmndk = QDates.QDate(2017,10,1) #仏滅
smtk = QDates.QDate(2017,11,1) #大安
sws = QDates.QDate(2017,12,1) #赤口

monthnames = ["睦月","如月","弥生","卯月","皐月","閏皐月","水無月",
              "文月","葉月","長月","神無月","霜月","師走"]
daysofweek = [QDates.先勝,QDates.友引,QDates.先負,QDates.仏滅,QDates.大安,QDates.大安,QDates.赤口,
              QDates.先勝,QDates.友引,QDates.先負,QDates.仏滅,QDates.大安,QDates.赤口]
dows = ["先勝","友引","先負","仏滅","大安","大安","赤口",
        "先勝","友引","先負","仏滅","大安","赤口"]
daysinmonth = [29,30,29,30,29,29,30,29,30,29,30,30,30]

@testset "Name functions" begin
    for (i,qdt) in enumerate([mtk,ksrg,yyi,udk,stk,lstk,mndk,fdk,hdk,ngtk,kmndk,smtk,sws])
        let m = i > 5 ? (i - 1) : i
            @test QDates.month(qdt) == m
            @test QDates.monthname(qdt) == monthnames[i]
            @test QDates.monthname(m) == monthnames[i == 6 ? 5 : i]
            @test QDates.monthabbr(qdt) == monthnames[i]
            @test QDates.monthabbr(m) == monthnames[i == 6 ? 5 : i]
            @test QDates.dayofweek(qdt) == daysofweek[i]
            @test QDates.dayname(qdt) == QDates.dayname(QDates.dayofweek(qdt)) == dows[i]
            @test QDates.dayabbr(qdt) == QDates.dayabbr(QDates.dayofweek(qdt)) == dows[i]
            @test QDates.daysinmonth(qdt) == daysinmonth[i]
        end
    end
    @test QDates.is先勝(mtk)
end
@testset "leap years" begin
    @test QDates.isleapyear(QDates.QDate(2010)) == false
    @test QDates.isleapyear(QDates.QDate(2012)) == true
    @test QDates.isleapyear(QDates.QDate(2016)) == false
    @test QDates.isleapyear(QDates.QDate(2017)) == true

    @test QDates.daysinyear(QDates.QDate(2010)) == 354
    @test QDates.daysinyear(QDates.QDate(2017)) == 384

    @test QDates.monthsinyear(QDates.QDate(2010)) == 12
    @test QDates.monthsinyear(QDates.QDate(2017)) == 13
end
@testset "dayofyear" begin
    @test QDates.dayofyear(2017,1,1) == 1
    @test QDates.dayofyear(2017,5,1) == 119
    @test QDates.dayofyear(2017,5,true,1) == 148
    @test QDates.dayofyear(QDates.QDate(2017,1,1)) == 1
    @test QDates.dayofyear(QDates.QDate(2017,5,1)) == 119
    @test QDates.dayofyear(QDates.QDate(2017,5,true,1)) == 148
end
@testset "privilleged 2100 (broken)" begin
    # (til v0.1.0) Available just 2100/01/01 - 2100/12/01
    @test_broken QDates.daysinyear(QDates.QDate(2100)) == 326
    @test_broken QDates.monthsinyear(QDates.QDate(2100)) == 11
end
@testset "test for years 2100-2200" begin
    # Available 2100/01/01 - 2200/12/29
    @test QDates.daysinyear(QDates.QDate(2100)) == 354
    @test QDates.monthsinyear(QDates.QDate(2100)) == 12
    @test QDates.daysinyear(QDates.QDate(2110)) == 354
    @test QDates.monthsinyear(QDates.QDate(2110)) == 12
    @test QDates.daysinyear(QDates.QDate(2120)) == 384
    @test QDates.monthsinyear(QDates.QDate(2120)) == 13
    @test QDates.daysinyear(QDates.QDate(2130)) == 355
    @test QDates.monthsinyear(QDates.QDate(2130)) == 12
    @test QDates.daysinyear(QDates.QDate(2140)) == 355
    @test QDates.monthsinyear(QDates.QDate(2140)) == 12
    @test QDates.daysinyear(QDates.QDate(2150)) == 384
    @test QDates.monthsinyear(QDates.QDate(2150)) == 13
    @test QDates.daysinyear(QDates.QDate(2160)) == 354
    @test QDates.monthsinyear(QDates.QDate(2160)) == 12
    @test QDates.daysinyear(QDates.QDate(2170)) == 354
    @test QDates.monthsinyear(QDates.QDate(2170)) == 12
    @test QDates.daysinyear(QDates.QDate(2180)) == 384
    @test QDates.monthsinyear(QDates.QDate(2180)) == 13
    @test QDates.daysinyear(QDates.QDate(2190)) == 354
    @test QDates.monthsinyear(QDates.QDate(2190)) == 12
    @test QDates.daysinyear(QDates.QDate(2200)) == 354
    @test QDates.monthsinyear(QDates.QDate(2200)) == 12
end

end


text:
module JacobiSVD
using LinearAlgebra

export gesvj!, gejsv!, jsvd!

include("lapack.jl")

"""
    jsvd!(A; full::Bool = false, precondition::Bool = true)

Compute the singular value decomposition (SVD) of `A` and return an `SVD` object.

Uses a one-sided Jacobi SVD algorithm, or (normally) a preconditioned variant.
These schemes often produce more accurate small singular values than the
standard algorithms. Preconditioning often improves performance.
Usage is similar to [`LinearAlgebra.svd!`](@ref).
"""
function jsvd!(A::StridedMatrix{T}; full::Bool = false,
               precondition::Bool = true,
               ) where T<:LinearAlgebra.BlasFloat
    m,n = size(A)
    squat = m < n
    if precondition
        if squat
            A = Matrix(A')
            xfull = true
        else
            xfull = full
        end
        U, S, V, scales = gejsv!('F', xfull ? 'F' : 'U', 'V', A)
        S .= (scales[2]/scales[1]) .* S
        if squat
            U, V = V, U
            if !full
                V = V[1:n,1:m]
            end
        else
            if !full
                U = U[1:m,1:n]
            end
        end
    else
        if full
            throw(ArgumentError("full U is not available from gesvj!"))
        end
        U, S, V, scale, ranka = gesvj!('G', 'U', 'V', A)
        if ranka < n
            # CHECKME: are trailing entries in S zeroed out?
            @warn "matrix is rank deficient; only $ranka left vectors are valid"
        end
        S .= scale .* S
    end
    # CHECKME: should we materialize the transpose?
    LinearAlgebra.SVD(U,S,V')
end

end # module


text:
module SE

using DataFrames
using Random
using XLSX
using StructArrays
using StatsBase
using CSV
using Main.JOH
using JuMP
using JSON

""" create a variety of SSIT methods. Accept a parameter to multiply each time
limit by. """
function make_SSIT_methods(m=60; n_threads=6)
    [
        JOH.Matheur.SSIT.make_SSIT_method(
			[.001, .005, .01, .02, .05],
			[m*5, m*5, m*5, m*5, m*5],
			"even time", n_threads)
        JOH.Matheur.SSIT.make_SSIT_method(
			[.001, .001, .001, .001, .001],
			[m*5,m*5,m*5,m*5,m*5],
			"one tolerance", n_threads)
        JOH.Matheur.SSIT.make_SSIT_method(
			[.001, .001, .005, .005, .005],
			[m*5, m*5, m*5, m*5, m*5],
			"tight tolerances", n_threads)
        JOH.Matheur.SSIT.make_SSIT_method(
			[.001, .005, .01, .02, .05],
			[m*2, m*4, m*5, m*6, m*8],
            "increasing time", n_threads)
        JOH.Matheur.SSIT.make_SSIT_method(
			[.001, .005, .01, .02, .05],
			[m*8, m*6, m*5, m*4, m*2],
            "decreasing time", n_threads)
    ]
end

struct MethodProblemResult
	method_name
	problem_id
	lowest_gap
	highest_reached_tolerance
	total_time
	last_phase_time
	cplex_obj
	true_obj
	infeasibility
	n_phases
end

MethodProblemResult(method, problem, solution, last_row) = MethodProblemResult(
	method.name,
	problem.id.id,
	last_row[!, :gap],
	last_row[!, :rtol],
	last_row[!, :elapsed_time],
	last_row[!, :solve_time],
	last_row[!, :objective],
	solution.objective,
	solution.infeasibility,
	last_row[!, :index])

mutable struct ExperimentResults{T}
	problem_ids::Vector{T}
	SSIT_phases::Vector{DataFrame}
	method_problem_results::Vector{MethodProblemResult}
end

ExperimentResults() = ExperimentResults([], [], [])

function flatten_ssit(df::DataFrame, tolerances)
	times = []
	gaps = []
	objectives = []

	n_rows = length(df[!, 1])

	for i in 1:n_rows
		push!(times, df[i, :][:elapsed_time])
		push!(gaps, df[i, :][:gap])
		push!(objectives, df[i, :][:objective])
	end

	flat_df = DataFrame()
	for i in 1:n_rows
		flat_df[!, Symbol("phase $i time")] = [times[i]]
		flat_df[!, Symbol("phase $i gap")] = [gaps[i]]
		flat_df[!, Symbol("phase $i obj")] = [objectives[i]]
	end
	flat_df[!, :termination] = [last(df)[:term_stat]]
	flat_df[!, :lowest_gap] = [last(df)[:gap]]
	flat_df
end

function include_aux_data(df::DataFrame, method, problem_id)
	df[!, :method] .= method.name
	for field in fieldnames(typeof(problem_id))
		val = getfield(problem_id, field)
		df[!, Symbol("problem_$(field)")] = [val]
	end
	df
end

function include_sol_data(df, ssit_phases, model)
	lp = last(ssit_phases, 1)
	try
		df[!, :objective] = [objective_value(model)]
	catch e
		df[!, :objective] = [-1]
	end
	df
end

_get_id(p) = try p["id"]  catch LoadError p.id end
function summarize_ssit(ssit_phases::DataFrame, method, problem, model)
	df = include_aux_data(flatten_ssit(ssit_phases, method.tolerances),
		method, _get_id(problem))
	df = include_sol_data(df, ssit_phases, model)
	df
end

function generate_comparison_data(
		method::JOH.Matheur.SSIT.SSIT_method,
		problems::Vector{T},
		mips_model; results_dir="results") where T

	results = []
	for problem in problems
		model = mips_model(problem)
		ssit_phases = JOH.Matheur.evaluate(model, method)
		result_df = summarize_ssit(ssit_phases, method, problem, model)
		CSV.write("$(results_dir)/$(problem.id.id).csv", result_df)

		push!(results, result_df)
	end

	results
end

function generate_comparison_data2(
		method::JOH.Matheur.SSIT.SSIT_method,
		problems::Vector{T},
		; results_dir="results") where T

	results = []
	for problem in problems
		ssit_phases = JOH.Matheur.evaluate(problem.model, method)
		result_df = summarize_ssit(ssit_phases, method, problem, problem.model)
		CSV.write("$(results_dir)/$(_get_id(problem.id)).csv", result_df)

		push!(results, result_df)
		problem.model = nothing
	end

	results
end

function save_model(m, m_path, s_path)
	write_to_file(m, m_path)
	try
		open(s_path, "w") do f
			print(f, json(value.(all_variables(m))))
		end
	catch e
		if !(isa(e, JuMP.OptimizeNotCalled))
			rethrow()
		end
	end
end

function read_solution(s_path::String)
	try
		open(s_path, "r") do f
			JSON.parse(read(f,String))
		end
	catch e
		false
	end
end

function log_ssit_run(m::JuMP.Model, method, res_dir::String, optimizer,
		getdettime)
	JOH.Matheur.set_threads!(m, method.num_threads)

	for i in 1:length(method.tolerances)
		#create a directory to store this phase's information
		phase_dir = joinpath(res_dir, "$(i)")
		mkpath(phase_dir)

		#generate paths to data files
		m_path, s_path, r_path = map(n->joinpath(phase_dir, n),
			["model.mps", "start_sol.json", "results.json"])

		#save the model and starting solution
		save_model(m, m_path, s_path) #TODO: save the method and phase as well

		#delete the current model, then replace from the saved record
		#this is to make starting from the saved files deterministic
		m = nothing
		m = read_from_file(m_path)
		solution = read_solution(s_path)
		set_optimizer(m, optimizer)

		#update the model according to the SSIT phase parameters
		JOH.Matheur.set_tolerance!(m, method.tolerances[i])
		JOH.Matheur.set_time!(m, method.times[i])
		JOH.Matheur.set_threads!(m, method.num_threads)

		if solution != false
			set_start_value.(all_variables(m), solution)
		end

		# run the optimization, and record the elapsed time
		start_time = time()
		optimize!(m)
		end_time = time()
		elapsed_time = end_time - start_time

		#make sure julia deletes the C optimizers memory
		GC.gc() #this doesn't always happen automatically

		row = JOH.Matheur.get_DF_row(m, elapsed_time=elapsed_time, index=i,
			getdettime=getdettime)

		GC.gc() #this doesn't always happen automatically

		open(r_path, "w") do f
			print(f, json(row))
		end

		if termination_status(m) == MOI.OPTIMAL || termination_status(m) ==
				MOI.INFEASIBLE
			break
		end
	end
end

function log_method_results(
		method::JOH.Matheur.SSIT.SSIT_method,
		problems::Vector{T},
		mips_model, res_dir, optimizer, getdettime) where T <: JOH.Problem

	rm(res_dir, force=true, recursive=true)
	mkpath(res_dir)

	for problem in problems
		problem_dir = joinpath(res_dir, "$(problem.id.id)")
		mkdir(problem_dir)
		log_ssit_run(mips_model(problem), method, problem_dir, optimizer,
			getdettime)
	end
end

function ba_rep(ba)
	join([b == 1 ? "1" : "0" for b in ba], "")
end

end

max_stars_repo_name: Dibillilia/JuliaOptHeuristics

max_issues_repo_name: Dibillilia/JuliaOptHeuristics

max_forks_repo_name: Dibillilia/JuliaOptHeuristics


text:
"""
`PDBML <: FileFormat`

Protein Data Bank Markup Language (PDBML), a representation of PDB data in XML format.
"""
struct PDBML <: FileFormat end

function _get_text(elem, name)
    sub = find_element(elem, name)
    if sub !== nothing
        return(content(sub))
    end
    throw(ErrorException("There is not $name for $elem"))
end

function _get_ins_code(elem)
    sub = find_element(elem, "pdbx_PDB_ins_code")
    if sub !== nothing
        return(content(sub))
    else
        return("")
    end
end

function _get_atom_iterator(document::LightXML.XMLDocument)
    pdbroot = root(document)
    child_elements(get_elements_by_tagname(pdbroot, "atom_siteCategory")[1])
end

"Used for parsing a PDB file into `Vector{PDBResidue}`"
function _generate_residues(residue_dict::OrderedDict{PDBResidueIdentifier, Vector{PDBAtom}},
                            occupancyfilter::Bool=false)
    if occupancyfilter
        return( PDBResidue[ PDBResidue(k, bestoccupancy(v)) for (k,v) in residue_dict ] )
    else
        return( PDBResidue[ PDBResidue(k, v) for (k,v) in residue_dict ] )
    end
end

"""
`parse(pdbml, ::Type{PDBML}; chain=All, model=All, group=All, atomname=All, onlyheavy=false, label=true, occupancyfilter=false)`

Reads a `LightXML.XMLDocument` representing a pdb file.
Returns a list of `PDBResidue`s (view `MIToS.PDB.PDBResidues`).
Setting `chain`, `model`, `group`, `atomname` and `onlyheavy` values can be used to select
of a subset of all residues. If not set, all residues are returned. If the keyword
argument `label` (default: `true`) is `false`,the **auth_** attributes will be use instead
of the **label_** attributes for `chain`, `atom` and residue `name` fields. The **auth_**
attributes are alternatives provided by an author in order to match the
identification/values used in the publication that describes the structure. If the
keyword argument `occupancyfilter` (default: `false`) is `true`, only the atoms with the
best occupancy are returned.
"""
function Base.parse(pdbml::LightXML.XMLDocument, ::Type{PDBML};
                    chain::Union{String,Type{All}} = All,
                    model::Union{String,Type{All}} = All,
                    group::Union{String,Type{All}} = All,
                    atomname::Union{String,Type{All}} = All,
                    onlyheavy::Bool=false,
                    label::Bool=true,
                    occupancyfilter::Bool=false)

    residues = Vector{PDBResidue}()

    prefix = label ? "label" : "auth"
    chain_attribute = string(prefix, "_asym_id")
    atom_attribute = string(prefix, "_atom_id")
    comp_attribute = string(prefix, "_comp_id")

    residue_id = PDBResidueIdentifier("", "", "", "", "", "")
    atoms = _get_atom_iterator(pdbml)
    for atom in atoms

        atom_name = _get_text(atom, atom_attribute)
        if !_is(atom_name, atomname)
            continue
        end

        element = _get_text(atom, "type_symbol")
        if onlyheavy && element=="H"
            continue
        end

        atom_group = _get_text(atom, "group_PDB")
        if !_is(atom_group, group)
            continue
        end

        atom_chain = _get_text(atom, chain_attribute)
        if !_is(atom_chain, chain)
            continue
        end

        atom_model = _get_text(atom, "pdbx_PDB_model_num")
        if !_is(atom_model, model)
            continue
        end

        PDBe_number = _get_text(atom, "label_seq_id")

        #  Residue_No  _atom_site.auth_seq_id
        #  Ins_Code    _atom_site.pdbx_PDB_ins_code
        PDB_number = string(_get_text(atom, "auth_seq_id"), _get_ins_code(atom))
        name = _get_text(atom, comp_attribute)

        if (residue_id.PDBe_number != PDBe_number) ||
                (residue_id.number != PDB_number) ||
                (residue_id.name != name) ||
                (residue_id.chain != atom_chain) ||
                (residue_id.group != atom_group) ||
                (residue_id.model != atom_model)

            n_res = length(residues)
            if occupancyfilter && n_res > 0
                residues[n_res].atoms = bestoccupancy(residues[n_res].atoms)
            end

            residue_id = PDBResidueIdentifier(PDBe_number,
                                              PDB_number,
                                              name,
                                              atom_group,
                                              atom_model,
                                              atom_chain)
            push!(residues, PDBResidue(residue_id, Vector{PDBAtom}()))
        end

        x = parse(Float64, _get_text(atom, "Cartn_x"))
        y = parse(Float64, _get_text(atom, "Cartn_y"))
        z = parse(Float64, _get_text(atom, "Cartn_z"))
        occupancy = parse(Float64, _get_text(atom, "occupancy"))
        B = _get_text(atom, "B_iso_or_equiv")

        push!(residues[end].atoms, PDBAtom(Coordinates(x,y,z),
                                           atom_name,
                                           element,
                                           occupancy,
                                           B))
    end

    if occupancyfilter
        residues[end].atoms = bestoccupancy(residues[end].atoms)
    end

    residues
end

# function Base.parse(pdbml::LightXML.XMLDocument, ::Type{PDBML};
#                     chain::Union{String,Type{All}} = All,
#                     model::Union{String,Type{All}} = All,
#                     group::Union{String,Type{All}} = All,
#                     atomname::Union{String,Type{All}} = All,
#                     onlyheavy::Bool=false,
#                     label::Bool=true,
#                     occupancyfilter::Bool=false)
#
#     residue_dict = OrderedDict{PDBResidueIdentifier, Vector{PDBAtom}}()
#
#     prefix = label ? "label" : "auth"
#     chain_attribute = string(prefix, "_asym_id")
#     atom_attribute = string(prefix, "_atom_id")
#     comp_attribute = string(prefix, "_comp_id")
#
#     atoms = _get_atom_iterator(pdbml)
#     for atom in atoms
#
#         atom_group = _get_text(atom, "group_PDB")
#         atom_model = _get_text(atom, "pdbx_PDB_model_num")
#         atom_chain = _get_text(atom, chain_attribute)
#         atom_name = _get_text(atom, atom_attribute)
#         element = _get_text(atom, "type_symbol")
#
#         if  _is(atom_group, group) && _is(atom_chain,chain) &&
#                 _is(atom_model, model) && _is(atom_name,atomname) &&
#                 (!onlyheavy || element!="H")
#
#             PDBe_number = _get_text(atom, "label_seq_id")
#
#             #  Residue_No  _atom_site.auth_seq_id
#             #  Ins_Code    _atom_site.pdbx_PDB_ins_code
#             PDB_number = string(_get_text(atom, "auth_seq_id"), _get_ins_code(atom))
#             name = _get_text(atom, comp_attribute)
#             x = parse(Float64, _get_text(atom, "Cartn_x"))
#             y = parse(Float64, _get_text(atom, "Cartn_y"))
#             z = parse(Float64, _get_text(atom, "Cartn_z"))
#             occupancy = parse(Float64, _get_text(atom, "occupancy"))
#             B = _get_text(atom, "B_iso_or_equiv")
#
#             residue_id = PDBResidueIdentifier(PDBe_number, PDB_number, name, atom_group, atom_model, atom_chain)
#             atom_data  = PDBAtom(Coordinates(x,y,z), atom_name, element, occupancy, B)
#
#             value = get!(residue_dict, residue_id, PDBAtom[])
#             push!(value, atom_data)
#         end
#     end
#     _generate_residues(residue_dict, occupancyfilter)
# end

# Download PDB
# ============

function _inputnameforgzip(outfile)
    if endswith(outfile, ".gz")
        return(outfile)
    end
    string(outfile, ".gz")
end

_file_extension(format::Type{PDBML}) = ".xml.gz"
_file_extension(format::Type{PDBFile}) = ".pdb.gz"

"""
It downloads a gzipped PDB file from PDB database.
It requires a four character `pdbcode`.
Its default `format` is `PDBML` (PDB XML) and It uses the
`baseurl` "http://www.rcsb.org/pdb/files/".
`filename` is the path/name of the output file.
This function calls `MIToS.Utils.download_file` that calls `HTTP.open`.
You can use keyword arguments from `HTTP.request` (e.g. `redirect`).
Use the `headers` keyword argument to pass a `Dict{String, String}`
with the header information.
"""
function downloadpdb(pdbcode::String;
                     format::Type{T} = PDBML,
                     filename::String= uppercase(pdbcode)*_file_extension(format),
                     baseurl::String = "http://www.rcsb.org/pdb/files/",
                     headers::Dict{String, String} = Dict{String, String}(),
                     kargs...) where T <: FileFormat
    if check_pdbcode(pdbcode)
        pdbfilename = uppercase(pdbcode) * _file_extension(format)
        filename = _inputnameforgzip(filename)
        sepchar = endswith(baseurl,"/") ? "" : "/";
        download_file(string(baseurl,sepchar,pdbfilename), filename; headers=headers, kargs...)
    else
        throw(ErrorException("$pdbcode is not a correct PDB code"))
    end
    filename
end

# RESTful PDB interface
# =====================

function _graphql_query(pdbcode::String)
    """
    {
      entry(entry_id: "$pdbcode") {
        entry {
          id
        }
        rcsb_entry_info {
          experimental_method
          assembly_count
          resolution_combined
        }
        rcsb_accession_info {
          initial_release_date
        }
        polymer_entities {
          rcsb_polymer_entity_container_identifiers {
            entity_id
            auth_asym_ids
          }
          entity_poly {
            rcsb_entity_polymer_type
          }
        }
      }
    }
    """
end

function _pdbheader(pdbcode::String; kargs...)
    pdbcode = uppercase(pdbcode)
    if check_pdbcode(pdbcode)
        String(
            HTTP.request(
                "GET",
                "https://data.rcsb.org/graphql";
                query = ["query" => _graphql_query(pdbcode)],
                kargs...
                ).body
            )
    else
        throw(ErrorException("$pdbcode is not a correct PDB code"))
    end
end

"""
It downloads a JSON file containing the PDB header information.
"""
function downloadpdbheader(pdbcode::String; filename::String=tempname(), kargs...)
    open(filename, "w") do fh
        write(fh, _pdbheader(pdbcode; kargs...))
    end
    filename
end

"""
Access general information about a PDB entry (e.g., Header information) using the
GraphQL interface of the PDB database. It parses the JSON answer into a Dict.
"""
function getpdbdescription(pdbcode::String; kargs...)
    JSON.parse(_pdbheader(pdbcode; kargs...))["data"]["entry"]
end


text:
module ExternalModule

function foo()
    println("ExternalModule.foo()")
end

end

max_stars_repo_name: erwanlecarpentier/JuliaSandbox

max_issues_repo_name: erwanlecarpentier/JuliaSandbox

max_forks_repo_name: erwanlecarpentier/JuliaSandbox


text:
using UnivariateFunctions
using Test

# Run tests

println("Test of UnivariateFunctions package.")
@time @test include("1_basic_tests.jl")
@time @test include("2_date_tests.jl")
@time @test include("3_piecewise_tests.jl")
@time @test include("4_schumaker_test.jl")
@time @test include("5_interpolation_test.jl")
@time @test include("6_test_regressions.jl")