text:
# This code is part of QuantumCircuits.
#
# (C) Copyright Rafał Pracht 2022.
#
# This code is licensed under the Apache License, Version 2.0. You may
# obtain a copy of this license in the LICENSE.txt file in the root directory
# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
#
# Any modifications or derivative works of this code must retain this
# copyright notice, and modified files need to carry a notice indicating
# that they have been altered from the originals.

module Circuit

using QuantumCircuits.QCircuits.QBase
using QuantumCircuits.QCircuits.Gates
using QuantumCircuits.QCircuits.Gates: ParamT, appendparams!, appendRandParams!
using QuantumCircuits.QCircuits.ComplexGates
using QuantumCircuits.QCircuits.ComplexGates: U4_params
using QuantumCircuits.QCircuits.OtherGates
using QuantumCircuits.QCircuits.Instructions
using QuantumCircuits.QCircuits.Registers
using QuantumCircuits.QCircuits.Qiskit
using QuantumCircuits.QCircuits.Graph
using QuantumCircuits.QCircuits.Math
using QuantumCircuits.QCircuits.Qiskit: IS_QISKIT, NoQiskitError

using MacroTools

using LinearAlgebra

using CBOOCall: @cbooify

import QuantumCircuits.QCircuits.QBase: add!, tomatrix, setparameters!, simplify,
       standardGateError, decompose, measure!, bindparameters!
import Base.show

export QCircuit, getparameters, getRandParameters, toString, setClassicalRegister!, toPythonQiskitCode

"Nothing function"
const nop = () -> nothing

@doc raw"""
    QCircuit - Quantum Circuit structure
The fundamental element of quantum computing is the quantum circuit. A quantum circuit is a computational routine consisting of coherent quantum operations on quantum data, such as qubits. It is an ordered sequence of quantum gates, measurements and resets, which may be conditioned on real-time classical computation. A set of quantum gates is said to be universal if any unitary transformation of the quantum data can be efficiently approximated arbitrarily well as as sequence of gates in the set. Any quantum program can be represented by a sequence of quantum circuits and classical near-time computation.

In QuantumCircuits, this core element is represented by the QCircuit struct.
Below is an example of a quantum circuit that makes a three-qubit GHZ state defined as:
```math
|\psi\rangle = \frac{|000\rangle + |111\rangle}{\sqrt{2}}
```

```julia
using QuantumCircuits
qc = QCircuit(3)
qc.h(0)
qc.cx(0, 1)
qc.cx(1, 2)
```

```
      ┌───┐
q0_0: ┤ H ├──■───────
      └───┘┌─┴─┐
q0_1: ─────┤ X ├──■──
           └───┘┌─┴─┐
q0_2: ──────────┤ X ├
                └───┘
c0: 3/═══════════════

```

"""
mutable struct QCircuit <: QuantumCircuit
    qubits::Int
    qRegisters::Vector{QuantumRegister}
    cRegisters::Vector{ClassicalRegister}
    dcg::DirectedGraph{Int}
    has_code::Bool
    code::Vector{QuantumGate}
    vqubits::Vector{Qubit}
    vcbits::Vector{Cbit}
    measures::Vector{Pair{Qubit, Cbit}}
    measures_matrix::Matrix{Float64}

    function QCircuit(qRegs::Vector{QuantumRegister}, cRegs::Vector{ClassicalRegister}; inline_optimization=true)
        n = sum([length(i) for i in qRegs])

        qc = new(n, qRegs, cRegs, DirectedGraph{Int}(n, inline_optimization), false, QuantumGate[], Qubit[], ClassicalRegister[], Pair{Qubit, Cbit}[], eye(2^n))

        # assign the qregister to circuit
        for r in qRegs
            append!(qc.vqubits, r.bits)
        end
        for (i, q) in enumerate(qc.vqubits)
            setid(q, i - 1)
        end

        # and classical one
        for r in cRegs
            append!(qc.vcbits, r.bits)
        end
        for (i, q) in enumerate(qc.vcbits)
            setid(q, i - 1)
        end

        return qc
    end
end
QCircuit(n::Integer; inline_optimization=true) = QCircuit([QuantumRegister(n)], [ClassicalRegister(n)], inline_optimization=inline_optimization)
QCircuit(qReg::QuantumRegister, cReg::ClassicalRegister; inline_optimization=true) = QCircuit([qReg], [cReg], inline_optimization=inline_optimization)
QCircuit(reg::QuantumRegister; inline_optimization=true) = QCircuit([reg], ClassicalRegister[], inline_optimization=inline_optimization)
QCircuit(regs::Vector{QuantumRegister}; inline_optimization=true) = QCircuit(regs, [ClassicalRegister(sum([length(i) for i in regs]))], inline_optimization=inline_optimization)
QCircuit(regs::Vector{QuantumRegister}, cReg::ClassicalRegister; inline_optimization=true) = QCircuit(regs, [cReg], inline_optimization=inline_optimization)
function QCircuit(qc::QCircuit)
    qregs = [QuantumRegister(length(r), r.name) for r in qc.qRegisters]
    cregs = [ClassicalRegister(length(r), r.name) for r in qc.cRegisters]

    QCircuit(qregs, cregs)
end

###################################################################################
@cbooify QCircuit (x, sx, y, z, h, cx, s, sdg, t, tdg, u, u3, rx, ry, rz, rzx, u4, barrier, measure)

"Add function macro"
macro addfunction(name, gate)
    eval(quote
        $name(qc::QCircuit, args...) = add!(qc, $gate, args...)
    end)
end
macro addfunction1param(name, gate)
    eval(quote
        $name(qc::QCircuit, q, θ=ParameterT(rand()*2π)) = add!(qc, $gate, q, θ)
    end)
end
macro addfunction3param(name, gate)
    eval(quote
        $name(qc::QCircuit, q, θ=ParameterT(rand()*2π), ϕ=ParameterT(rand()*2π), λ=ParameterT(rand()*2π)) = add!(qc, $gate, q, θ, ϕ, λ)
    end)
end

@addfunction(x, X)
@addfunction(sx, Sx)
@addfunction(y, Y)
@addfunction(z, Z)
@addfunction(h, H)
@addfunction(cx, CX)
@addfunction(s, S)
@addfunction(sdg, Sd)
@addfunction(t, T)
@addfunction(tdg, Td)
@addfunction(u, U)

@addfunction1param(rx, Rx)
@addfunction1param(ry, Ry)
@addfunction1param(rz, Rz)

@addfunction3param(u3, U3)

rzx(qc::QCircuit, q1, q2, θ=ParameterT(rand()*2π)) = add!(qc, Rzx, q1, q2, θ)
u4(qc::QCircuit, q1, q2, params=[ParameterT(rand() * 2π) for i in 1:U4_params]) = add!(qc, U4, q1, q2, params)

barrier(qc::QCircuit) = add!(qc, Barrier(qc.vqubits))
measure(qc::QCircuit, q, c) = measure!(qc, q, c)

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

function getCode(c::QCircuit)
    if c.has_code
        return c.code
    else
        c.code = to_vector(c.dcg)
        c.has_code = true

        return c.code
    end
end

Base.:(==)(c1::QCircuit, c2::QCircuit) = c1.qubits == c2.qubits && getCode(c1) == getCode(c2) && c1.measures_matrix == c2.measures_matrix
Base.hash(c::QCircuit, h::UInt) = hash((c.qubits, getCode(code), c.measures_matrix), h)

"Add gate to circuit"
function add!(qc::QCircuit, gate::QuantumGate)
    add!(qc.dcg, gate)
    qc.has_code = false
    nothing
end
function add!(qc::QCircuit, gate::Type{<:QuantumGate}, qubit::Qubit)
    @assert qubit in qc.vqubits "Qubits is out of circuit qubits range."

    add!(qc.dcg, gate(qubit))
    qc.has_code = false
    nothing
end
function add!(qc::QCircuit, gate::Type{<:QuantumGate}, args::Tuple)
    add!(qc.dcg, gate(args...))
    qc.has_code = false
    nothing
end
function add!(qc::QCircuit, gate::Type{<:QuantumGate}, qubit::Integer)
    @assert qubit < qc.qubits "Qubits $qubit is out of circuit qubits range $(qc.qubits)."

    add!(qc, gate, qc.vqubits[qubit + 1])

    nothing
end
function add!(qc::QCircuit, gate::Type{<:QuantumGate}, qubits::AbstractVector)
    for qubit in qubits
        add!(qc, gate, qubit)
    end
end
function add!(qc::QCircuit, gate::Type{<:QuantumGate}, qubit1::Qubit, qubit2::Qubit)
    @assert qubit1 in qc.vqubits "Qubits is out of circuit qubits range."
    @assert qubit2 in qc.vqubits "Qubits is out of circuit qubits range."

    add!(qc.dcg, gate(qubit1, qubit2))
    qc.has_code = false
    nothing
end
function add!(qc::QCircuit, gate::Type{<:QuantumGate}, qubit1::Integer, qubit2::Integer)
    @assert qubit1 < qc.qubits "Qubits $qubit1 is out of circuit qubits range $(qc.qubits)."
    @assert qubit2 < qc.qubits "Qubits $qubit2 is out of circuit qubits range $(qc.qubits)."

    add!(qc, gate, qc.vqubits[qubit1 + 1], qc.vqubits[qubit2 + 1])

    nothing
end
function add!(qc::QCircuit, gate::Type{<:QuantumGate}, qubit1::Qubit, qubit2::Qubit, params)
    @assert qubit1 in qc.vqubits "Qubits is out of circuit qubits range."
    @assert qubit2 in qc.vqubits "Qubits is out of circuit qubits range."

    if isnothing(params)
        add!(qc.dcg, gate(qubit1, qubit2))
    else
        add!(qc.dcg, gate(qubit1, qubit2, params))
    end
    qc.has_code = false
    nothing
end
function add!(qc::QCircuit, gate::Type{<:QuantumGate}, qubit1::Integer, qubit2::Integer, params)
    @assert qubit1 < qc.qubits "Qubits $qubit1 is out of circuit qubits range $(qc.qubits)."
    @assert qubit2 < qc.qubits "Qubits $qubit2 is out of circuit qubits range $(qc.qubits)."

    add!(qc, gate, qc.vqubits[qubit1 + 1], qc.vqubits[qubit2 + 1], params)

    nothing
end

function add!(qc::QCircuit, gate::Type{<:QuantumGate}, qubit::Qubit, θ)
    @assert qubit in qc.vqubits "Qubits is out of circuit qubits range."

    add!(qc.dcg, gate(qubit, θ))
    qc.has_code = false
    nothing
end
function add!(qc::QCircuit, gate::Type{<:QuantumGate}, qubit::Integer, θ)
    @assert qubit < qc.qubits "Qubits $qubit is out of circuit qubits range $(qc.qubits)."

    add!(qc, gate, qc.vqubits[qubit + 1], θ)

    nothing
end

function add!(qc::QCircuit, gate::Type{<:QuantumGate}, qubits::Vector, θ)
    for qubit in qubits
        add!(qc, gate, qubit, θ)
    end
end

function add!(qc::QCircuit, gate::Type{<:QuantumGate}, qubit::Qubit, θ, ϕ, λ)
    @assert qubit in qc.vqubits "Qubits is out of circuit qubits range."

    add!(qc.dcg, gate(qubit, θ, ϕ, λ))
    qc.has_code = false
    nothing
end
function add!(qc::QCircuit, gate::Type{<:QuantumGate}, qubit::Integer, θ, ϕ, λ)
    @assert qubit < qc.qubits "Qubits $qubit is out of circuit qubits range $(qc.qubits)."

    add!(qc, gate, qc.vqubits[qubit + 1], θ, ϕ, λ)

    nothing
end
function add!(qc::QCircuit, gate::Type{<:QuantumGate}, qubits::AbstractVector, θ, ϕ, λ)
    for qubit in qubits
        add!(qc, gate, qubit, ParameterT(θ), ParameterT(ϕ), ParameterT(λ))
    end
end
function add!(qc::QCircuit, gates::Vector{<:QuantumGate})
    for g in gates
        add!(qc, g)
    end
    nothing
end

"Add measures to circuit"
function measure!(qc::QCircuit, qubit::Qubit, cbit::Cbit, setMatrix::Bool=true)
    push!(qc.measures, qubit => cbit)

    if setMatrix
        setMeasureMatrix!(qc)
    end
end

"Add measures to circuit"
function measure!(qc::QCircuit, qubit::Integer, cbit::Integer, setMatrix::Bool=true)
    @assert qubit < qc.qubits "Qubits $qubit is out of circuit qubits range $(qc.qubits)."
    @assert cbit < length(qc.vcbits) "Qubits $qubit is out of circuit qubits range $(length(qc.vcbits))."

    measure!(qc, qc.vqubits[qubit + 1], qc.vcbits[cbit + 1], setMatrix)
end
function measure!(qc::QCircuit, qubits::AbstractVector{<:Integer}, cbits::AbstractVector{<:Integer})
    @assert length(qubits) == length(cbits) "The length of qubits and classical bits should be equal."

    for (i, j) in zip(eachindex(qubits), eachindex(cbits))
        measure!(qc, qubits[i], cbits[j], false)
    end

    setMeasureMatrix!(qc)
end

function setMeasureMatrix!(qc::QCircuit)
    measured = [(i-1, v) for (i, v) in enumerate(sort!([getid(q)  for (q, c) in qc.measures]))]
    measured_qubits = length(measured)
    mes_matrix = zeros(2^measured_qubits, 2^qc.qubits)
    for i in 0:(2^qc.qubits-1)
        bin = digits(i, base=2, pad=qc.qubits)

        # calculate new index
        ni = 0
        for (idx, m) in measured
            if bin[m + 1] == 1
                ni += 2^idx
            end
        end
        mes_matrix[ni+1, i+1] = 1.0
    end

    qc.measures_matrix = mes_matrix
end

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

"Add the classical register to circuit."
function setClassicalRegister!(qc::QCircuit, cr::ClassicalRegister)
    @assert qc.cRegisters == ClassicalRegister[] "Unable to set register."
    @assert qc.vcbits == Cbit[] "Internal error :P"

    push!(qc.cRegisters, cr)
    append!(qc.vcbits, cr.bits)
    for (i, q) in enumerate(qc.vcbits)
        setid(q, i - 1)
    end
end

"Function to convert QCircuit to Qiskit circuit"
function toQiskit(circuit::QCircuit)
    if !IS_QISKIT
        throw(NoQiskitError("Unable to convert to Qiskit, qiskit is not installed."))
    end

    if isempty(circuit.cRegisters)
        qc = QiskitCircuit(circuit.qRegisters, nothing)
    else
        qc = QiskitCircuit(circuit.qRegisters, circuit.cRegisters)
    end

    for gate in getCode(circuit)
        add!(qc, gate)
    end

    for (q, c) in circuit.measures
        measure!(qc, q, c)
    end

    return qc
end

"Function to convert QCircuit to Python code using Qiskit library."
function toPythonQiskitCode(circuit::QCircuit)
    code = ""

    args = ""
    for (i, r) in enumerate(circuit.qRegisters)
        if isnothing(r.name)
            code *= "qr$i = QuantumRegister($(length(r.bits)))\n"
        else
            code *= "qr$i = QuantumRegister($(length(r.bits)), \"$(r.name)\")\n"
        end
        if i > 1
            args *= ", "
        end
        args *= "qr$i"
    end

    if !isempty(circuit.cRegisters)
        args *= ", "
        for (i, r) in enumerate(circuit.cRegisters)
            if isnothing(r.name)
                code *= "cr$i = ClassicalRegister($(length(r.bits)))\n"
            else
                code *= "cr$i = ClassicalRegister($(length(r.bits)), \"$(r.name)\")\n"
            end
            if i > 1
                args *= ", "
            end
            args *= "cr$i"
        end
    end

    code *= "qc = QuantumCircuit($args)\n"

    for gate in getCode(circuit)
        code *= getPythonCode("qc", gate)
    end

    for (q, c) in circuit.measures
        code *= "qc.measure($q, $(getid(c)))\n"
    end

    return code
end

"Show method"
function Base.show(io::IO, qc::QCircuit)
    if IS_QISKIT
        print(io, toQiskit(qc))
    else
        print(io, "QCircuit")
    end
end

"Function create new parameters vector"
function getparameters(qc::QCircuit)
    params = ParamT[]
    for gate in getCode(qc)
        appendparams!(params, gate)
    end

    return params
end

"Function create new random parameters vector"
function getRandParameters(qc::QCircuit)
    params = ParamT[]
    for gate in getCode(qc)
        appendRandParams!(params, gate)
    end

    return params
end

"Save the parameters to the quantum circuits."
function setparameters!(qc::QCircuit, params)
    i = 1
    for gate in getCode(qc)
        l = length(gate)
        if l > 0
            @views setparameters!(gate, params[i:(i+l-1)])
            i = i + l
        end
    end

    @assert i - 1 == length(params) "There is different number of parameters in circuit $(i-1) and in vector $(length(params))."
end

"Function bind the parameters, so now there only values."
function bindparameters!(qc::QCircuit)
    for gate in getCode(qc)
        bindparameters!(gate)
    end
end

function tomatrix(qc::QCircuit, params=nothing)
    i = 1
    unitary = LinearAlgebra.I * (1.0 + 0im)
    for gate in getCode(qc)
        l = length(gate)
        if !isnothing(params) && l > 0
            @views gate_unitary = tomatrix(qc.qubits, gate, params[i:(i+l-1)])
            # TODO gate_unitary = tomatrix(qc.qubits, gate, params[i:(i+l-1)])
            i = i + l
        else
            gate_unitary = tomatrix(qc.qubits, gate)
        end

        #println("GU: $gate_unitary")
        unitary = gate_unitary * unitary
        #println("U: $unitary")
    end
    @assert i - 1 == length(params) "There is different number of parameters in circuit $(i-1) and in vector $(length(params))."

    return unitary
end

function standardGateError(qc::QCircuit, params=nothing)
    err = 0.0
    i = 1
    for gate in getCode(qc)
        l = length(gate)
        if !isnothing(params) && l > 0
            @views gerr = standardGateError(gate, params[i:(i+l-1)])
            i = i + l
        else
            gerr = standardGateError(gate)
        end

        err +=  gerr
    end

    #return err / (100 * length(qc.code))
    return err
end

function Base.inv(qc::QCircuit)
    newqc = QCircuit(qc)

    for gate in reverse(getCode(qc))
        newgate = inv(gate)
        add!(newqc, typeof(newgate), getArgs(newgate)...)
    end

    return newqc
end

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

function decompose(qc::QCircuit)
    newqc = QCircuit(qc)

    for gate in getCode(qc)
        newgates = decompose(gate)
        for newgate in newgates
            add!(newqc, typeof(newgate), getArgs(newgate)...)
        end
    end

    return newqc
end

function simplify(qc::QCircuit)
    newqc = QCircuit(qc)

    for gate in getCode(qc)
        newgates = simplify(gate)
        for newgate in newgates
            add!(newqc, typeof(newgate), getArgs(newgate)...)
        end
    end

    return newqc
end

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

function Base.append!(qc::QCircuit, apqc::QCircuit)
    @assert length(qc.measures) == 0 "Unable append to measured circuit."
    @assert length(apqc.measures) == 0 "Unable append measured circuit."
    @assert length(qc.vqubits) == length(apqc.vqubits) "Unable append, the qubits are different"

    for (q, aq) in zip(qc.vqubits, apqc.vqubits)
        @assert q == aq "There are difference in qubits."
    end

    for c in getCode(apqc)
        add!(qc, typeof(c), getArgs(c)...)
    end
end

end  # module Circuit


text:
##############################################################
#
#  GENETIC TEMPERING ON A GRAPHICS PROCESSING UNIT (GT-GPU)
#
##############################################################
# Made by Thomas E. Baker (2018)
# See accompanying license with this program
# This code is native to the julia programming language (v0.6.1)
include("GPUpart.jl")

const L = 16 #lattice size on one side
const latsize = L*L #total number of sites...can be extended to non-square lattices
const J = 1 #1 is ferromagnetic

const Nloop = 1_000 #number of MC loops
const GPUloops = 1_000 #Number of instances of lattices on the GPU...an upper limit depends on the GPU's memory size

#adding a multiplier doesn't add much overhead, better to have higher M and lower Nloop on a laptop for heat reasons
#number of cycles for hysteresis loop
const delayloops = 5 #number of loops to check while thermalizing in generating hysteresis loop
const delaymult = 1 #store data for delaymult*Nloop iterations
#number of cycles in the final sweep
const relaxloops = 10 #delay on taking final sweep data
const maxmult = 10 #total Monte Carlo iterations: maxmult*Nloop iterations

const var = 1 #"1" = energy

const go_up = true #increase temperature or decrease?
const ordered = true #start from an ordered or a random state?
#####
const Tmin = 0.5 #2.25 #minimum temperature
const Tmax = 5.0 #2.35 #maximum temperature
const DeltaT = 0.1 #0.001 #temperature step
const numT = round(Int32,(Tmax-Tmin)/DeltaT)+1

#scan over a specific window...will add these points to the run schedule without duplicates
const Tc = true
const Tcmin = 2.2
const Tcmax = 2.3
const DeltaTc = 0.01
const numTc = round(Int32,(Tcmax-Tcmin)/DeltaTc)+1

result_buff=0

#information about the number of calcualtions for each lattice...can be removed
numcalcs =  Nloop*(2*delaymult+maxmult)*GPUloops #delay > Nloop ? Nloop*maxmult*GPUloops  : GPUloops*(Nloop*maxmult-delay)
println("single-thread equivalent iterations (per temperature): ",(@sprintf "%.3E" numcalcs))
println("Ideally, about ",floor(Int64,numcalcs/latsize)," lattice sites will be visited over ",GPUloops," threads and ",3," sweeps (each lattice site is visited about ",floor(Int64,numcalcs/latsize/GPUloops)," times per thread)")

#1536MB 1.2288E10
#4GB 3.2e+10
#8GB 6.4e+10
percentused = (2*4*GPUloops*Nloop)/1.2288E10
println("approximate upper limit of 1536MB GPU: ",(2*4*GPUloops*Nloop)/(8E6)," MB of 1536 MB (using ",GPUloops*Nloop," samples or ",100*percentused,"% of the total)")

#define functions for averages and standard deviation
avg(x::Array{Float64,1}...) = Float64[sum(x[p])/size(x[p],1) for p = 1:size(x,1)]
stDevsq(avgx::Array{Float64,1},x::Array{Float64,1}...) = Float64[sum(a->(x[p][a]-avgx[p])^2,1:size(x[p],1))/(size(x[p],1)*(size(x[p],1)-1)) for p = 1:size(x,1)]

#determines initial wavefunction psi across all GPU threads
if ordered
  psi = ones(Int8,GPUloops*latsize)
  else
  psi = Int8[(rand() > 0.5 ? Int8(1) : Int8(-1)) for j = 1:GPUloops*latsize]
end

device, ctx, queue = cl.create_compute_context() #creates a "context" for the GPU
psi_buff = cl.Buffer(Int8, ctx, (:rw, :copy), hostbuf=psi) #buffer wavefunction to GPU (can be written to and read from)
const intvec = Int32[Nloop,latsize,GPUloops,L,J]
intvec_buff = cl.Buffer(Int32, ctx, (:r, :copy), hostbuf=intvec)

#printed out data columns
println("#T       <E>     +/-E      <M>     +/-M      <|M|>   +/-|M|    C_v    +/-C_v   chi_|M|   +/-chi     B      +/-B")

#generates list of temperatures to run through
incrT = (go_up ? (Tmin:DeltaT:Tmax) : (Tmax:-DeltaT:Tmin))
incrTc = (go_up ? (Tcmin:DeltaTc:Tcmax) : (Tcmax:-DeltaTc:Tcmin))
runs = Float64[T for T = incrT]
Tc ? runs = sort(unique(vcat(runs,[T for T = incrTc]))) : 0.
go_up ? 0. : runs = reverse(runs)
const nruns = size(runs,1)
const interval = -nruns:nruns

data = zeros(Float64,Tc ? numTc : nruns,13) #output data
data[:,1] = Tc ? [T for T = (go_up ? (Tcmin:DeltaTc:Tcmax) : (Tcmax:-DeltaTc:Tcmin))] : runs
datafore = zeros(Float64,nruns,13) #output data going up in temperature
databack = zeros(Float64,nruns,13) #output data going down in temperature
datafore[:,1] = databack[:,1] = runs

savedpsi=Array{Int8,2}[zeros(Int8,GPUloops,latsize) for p = 1:nruns]#Array{Array{Int8,2}}(numT)
savedE=Array{Float64,1}[zeros(Float64,GPUloops) for p = 1:nruns]#Array{Array{Float64,1}}(numT)

geneticpsi=deepcopy(savedpsi) #saves wavefunctions for the final sweep

function dataline(dataArray::Array{Float64,2},m::Int64,t2::Float64,t1::Float64,gputime::Float64)
    println("T = ",dataArray[m,1],", loop time = ",t2-t1,", total GPU time = ",gputime)
    print((@sprintf "%.4f" dataArray[m,1])," ")
    for k = 2:2:size(dataArray,2)
      print((@sprintf "%.6f" dataArray[m,k])," ")
      print((@sprintf "%.6f" sqrt(dataArray[m,k+1]))," ")
    end
    println()
end

function savedata(dataArray::Array{Float64,2},m::Int64,t2::Float64,t1::Float64,gputime::Float64,values::Array{Float64,1},errors::Array{Float64,1})
  sizedatArray = size(dataArray,2)
  dataArray[m,2:2:sizedatArray-2] = deepcopy(values[1:size(values,1)-1])/latsize
  dataArray[m,3:2:sizedatArray-2] = deepcopy(errors[1:size(errors,1)-1])/latsize^2
  dataArray[m,sizedatArray-1:sizedatArray] = [values[size(values,1)],errors[size(errors,1)]]
  dataline(dataArray,m,t2,t1,gputime)
end

floatconvert(fcount::Float64,x::Array{Float32,1}) = [Float64(x[p])/fcount for p = 1:size(x,1)]

#start sweeps through temperature
counterloop2 = 0
breaksweep = false
@time for sweep = 1:2 #sweep loop: forward, backward, final loop
  breaksweep ? break : 0
for qq = interval
  if qq == 0
    if sweep == 2
      breaksweep = true
      break
    end
    continue
  end
  m = qq <= 0 ? nruns - abs(qq+1) : nruns + 1 - qq #site index
  temp = runs[m]
  if sweep > 1 && Tc
    (Tcmin-1E-10 <= temp <= Tcmax+1E-10) ? 0. : continue
  end
  sweep > 1 ? counterloop2 += 1 : 0
  tgpu = 0. #time spent computing on GPU
  t1 = time() #times the outerloop

  #initializing average values
  JavgE,JavgEsq,JavgM,JavgMsq,JavgMfour,JavgAbsM,JstdE,JstdM,JstdAbsM,sumsaveCount = 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.
  singleHalf = -Float32(temp*0.5) #stores -T/2 for future reference
  function GPUcall()
    randpos = rand(Int32(1):Int32(latsize),Nloop*GPUloops) #define random positions
    randNum = log.(rand(Float32,Nloop*GPUloops))*singleHalf #define random numbers
    randpos_buff = cl.Buffer(Int32,ctx, (:r, :copy), hostbuf=randpos) #buffer data to GPU
    randNum_buff = cl.Buffer(Float32,ctx, (:r, :copy), hostbuf=randNum)

    res_buff = cl.Buffer(Float32, ctx, :w, GPUloops*7) #"magic number" of 7 here is for the number of observables wanted (see OpenCL part)
    p = cl.Program(ctx, source=MCloop_kernel) |> cl.build! #generates opencl kernel
    k = cl.Kernel(p, "MCloop")
    t3 = time()
    queue(k, GPUloops, nothing, psi_buff, intvec_buff, randpos_buff, randNum_buff, res_buff) #actual run of the GPU
    t4 = time()
    tgpu += t4-t3
    return res_buff
  end
  sumsaveCount = zeros(GPUloops)
  #thermalization
  mult = sweep == 1 ? delayloops : relaxloops
  for pqr = 1:mult
    result_buff = GPUcall()
  end

  #sets all states to a state near the observable
  if sweep == 1
    if qq == -nruns && !ordered
      s = deepcopy(cl.read(queue, result_buff)) #transfer resulting data from the GPU to the CPU
      s = reshape(s,GPUloops,7) #reshapes from a vector into a GPUloops x 7 matrix
      findval = sum(s[:,var])/GPUloops
      foundval = sortperm([abs(s[b,var]-findval) for b = 1:GPUloops])

      psi = deepcopy(cl.read(queue, psi_buff))
      psi = reshape(psi,GPUloops,latsize)
      copypsi = deepcopy(psi[foundval[1],:])
      for w = 1:GPUloops
        psi[w,:] = deepcopy(copypsi)
      end
      psi_buff = cl.Buffer(Int8, ctx, (:rw, :copy), hostbuf=psi)
    end
  else
    psi_buff = cl.Buffer(Int8, ctx, (:rw, :copy), hostbuf=reshape(geneticpsi[m],GPUloops*latsize))
  end

  #actual sweeps
  mult = sweep == 1 ? delaymult : maxmult
  for pqr = 1:mult
    result_buff = GPUcall()
    s = deepcopy(cl.read(queue, result_buff)) #transfer resulting data from the GPU to the CPU
    s = reshape(s,GPUloops,7) #reshapes from a vector into a GPUloops x 7 matrix

    sumsaveCount += Float64[s[w,7]==0 ? Float64(1) : Float64(s[w,7]) for w = 1:GPUloops] #saves the number of counts from the GPU

    JavgE += s[:,1] #adds averages to the stored average values on the CPU
    JavgEsq += s[:,2]
    JavgM += s[:,3]
    JavgMsq += s[:,4]
    JavgAbsM += s[:,5]
    JavgMfour += s[:,6]
  end

  #stores values for each thread
  threadE,threadEsq,threadM,threadMsq,threadAbsM,threadMfour = Float64[],Float64[],Float64[],Float64[],Float64[],Float64[],Float64[]
  threadCv,threadChiM,threadB = Float64[],Float64[],Float64[]
  floatconvert(fcount::Float64,x::Array{Float32,1}) = Float64[Float64(x[p])/fcount for p = 1:size(x,1)]
  for b = 1:GPUloops
    fcount = Float64(sumsaveCount[b])#floating point count
    En,Esq,avgM,Msq,AbsM,Mfour = floatconvert(fcount,Float32[JavgE[b],JavgEsq[b],JavgM[b],JavgMsq[b],JavgAbsM[b],JavgMfour[b]])
    push!(threadE,En) #pushed value to vector
    push!(threadEsq,Esq)
    push!(threadM,avgM) #didn't need average of M below, so only the vector is taken
    push!(threadMsq,Msq)
    push!(threadAbsM,AbsM)
    push!(threadMfour,Mfour)
    push!(threadCv,(Esq-En^2)/temp^2)#calculate specific heat for each thread
    push!(threadChiM,(Msq-AbsM^2)/temp)#calculate magnetic susceptibility for each thread
    push!(threadB,1-(Mfour/Msq^2)/3)#calculate binder cumulant for each thread
  end

  #calculate averages of values between all threads
  averages = avg(threadE,threadM,threadAbsM,threadMfour,threadEsq,threadMsq,threadCv,threadChiM,threadB)
  JavgE,JavgM,JavgAbsM,JavgMfour,JavgEsq,JavgMsq,specificHeat,susceptibility,binder = deepcopy(averages)
  #calculates errors for a given thread (standard deviation of the mean)
  errors = stDevsq(averages,threadE,threadM,threadAbsM,threadMfour,threadEsq,threadMsq,threadCv,threadChiM,threadB)
  JstdE,JstdM,JstdAbsM,JstdMfour,JstdEsq,JstdMsq,errorCv,errorChiM,errorB = deepcopy(errors)

  #loads select states into geneticpsi for the final sweep
  if sweep == 1
    gcounter = 1
    if (qq < 0)
      savedpsi[m] = deepcopy(reshape(cl.read(queue, psi_buff),GPUloops,latsize))
      savedE[m] = deepcopy(threadE)
    else
        #generate a random set of wavefunctions inside the hysteresis loop
        minval = min(databack[m,var+1]-sqrt(databack[m,var+2]*GPUloops),datafore[m,var+1]-sqrt(datafore[m,var+2]*GPUloops))
        maxval = max(databack[m,var+1]+sqrt(databack[m,var+2]*GPUloops),datafore[m,var+1]+sqrt(datafore[m,var+2]*GPUloops))
        psi = deepcopy(cl.read(queue, psi_buff))
        psi = reshape(psi,GPUloops,latsize)
        while true
          w = rand(1:GPUloops)
          randnumber = round(rand()) == 1.0
          if randnumber
            if ((minval <= savedE[m][w]/latsize <= maxval) | (maxval <= savedE[m][w]/latsize <= minval))
              geneticpsi[m][gcounter,:] = deepcopy(savedpsi[m][w,:])
              gcounter += 1
            end
          else
            if ((minval <= threadE[w]/latsize <= maxval) | (maxval <= threadE[w]/latsize <= minval))
              geneticpsi[m][gcounter,:] = deepcopy(psi[w,:])
              gcounter += 1
            end
          end
          gcounter > GPUloops ? break : 0
        end
    end
  end

  valtosave = [JavgE,JavgM,JavgAbsM,specificHeat,susceptibility,binder]
  errtosave = [JstdE,JstdM,JstdAbsM,errorCv,errorChiM,errorB]
  #record data
  t2 = time()
  if sweep == 1
    qq < 0 ? savedata(datafore,m,t2,t1,tgpu,valtosave,errtosave) : savedata(databack,m,t2,t1,tgpu,valtosave,errtosave)
  else
    savedata(data,counterloop2,t2,t1,tgpu,valtosave,errtosave)
  end
end #end sweep loop

end #end temperature scan loop
#write data to files
function savedata(name::String,dataArray::Array{Float64,2})
  f = open(name,"w")
  write(f,"#T        <E>      +/-E       <M>      +/-M      <|M|>     +/-|M|    C_v      +/-C_v    chi_|M|    +/-chi     B       +/-B\n")
  for i = 1:size(dataArray,1)
    write(f,join([(@sprintf "%.4f" dataArray[i,1]),"  "]))
    for k = 2:2:size(dataArray,2)-1
      write(f,join([(@sprintf "%.6f" dataArray[i,k]),"  ",(@sprintf "%.6f" sqrt(dataArray[i,k+1]))," "]))
    end
    write(f,"\n")
  end
  close(f)
end

savedata("dataforward_L=$(L).txt",datafore)
savedata("databackward_L=$(L).txt",databack)
savedata("data_L=$(L).txt",data)

ext: jl

max_stars_repo_name: bakerte/GeneticTemperingMC

max_issues_repo_name: bakerte/GeneticTemperingMC

max_forks_repo_name: bakerte/GeneticTemperingMC


text:
using DataStructures: Queue, enqueue!, dequeue!

function makenodelist(adjacencylist)
    nNodes = length(adjacencylist)
    nodes = []
    for i = 1:nNodes push!(nodes,Node(i)) end # Initialize array of nodes
    i = 1
    for i_node in adjacencylist
        for neighbour in i_node
            push!(nodes[i].neighbours, nodes[neighbour])
        end
        i = i+1
    end
    return nodes
end

function printnodelist(nodelist)
    println("Skriver ut nodeliste med printnodelist")
    for node in nodelist
        neighbourlist = [neighbours.id for neighbours in node.neighbours]
        println("id: $(node.id), neighbours: $neighbourlist")
    end
end

function main(;n=5, degree=3)
    println("Kjører makenodelist")
    adjacencylist = generateadjacencylist(n, degree)
    @info "adjacencylist" adjacencylist
    nodelist = makenodelist(adjacencylist)
    printnodelist(nodelist)
end

generateadjacencylist(n, degree) = [rand(1:n, degree) for id = 1:n]

# Det kan være nyttig å kjøre mange tester for å se om programmet man har laget
# krasjer for ulike instanser
function runmanytests(;maxlistsize=15)
    # Kjører n tester for hver listestørrelse – én for hver grad
    for n = 1:maxlistsize
        for degree = 1:n
            adjacencylist = generateadjacencylist(n, degree)
            @info "Listelendge $n" n, degree #, adjacencylist
            makenodelist(adjacencylist)
        end
    end
end

function bfs!(nodes, start)
    for v in nodes
        v.color = "WHITE"
        v.distance = typemax(Int)
        #v.predecessor = Nothing
    end
    start.color = "GRAY"
    start.distance = 0
    #start.predecessor = Nothing
    Q = Queue{Node}()
    enqueue!(Q,start)
    while length(Q) > 0
        u = dequeue!(Q)
        for v in u.neighbours
            if v.color == "WHITE"
                v.color == "GRAY"
                v.distance = u.distance + 1
                v.predecessor = u
                enqueue!(Q,v)
            end
        end
        u.color = "BLACK"
        # node object as argument and checks if this is the node we're looking for
    end
    return nodes
end

main()
#runmanytests()

max_stars_repo_name: evenlwanvik-student/algdat

max_issues_repo_name: evenlwanvik-student/algdat

max_forks_repo_name: evenlwanvik-student/algdat


text:
#####
##### Gradient AD implementation using Tracker
#####

import .Tracker

struct TrackerGradientLogDensity{L} <: ADGradientWrapper
    ℓ::L
end

"""
    ADgradient(:Tracker, ℓ)
    ADgradient(Val(:Tracker), ℓ)

Gradient using algorithmic/automatic differentiation via Tracker.

This package has been deprecated in favor of Zygote, but we keep the interface available.
"""
ADgradient(::Val{:Tracker}, ℓ) = TrackerGradientLogDensity(ℓ)

Base.show(io::IO, ∇ℓ::TrackerGradientLogDensity) = print(io, "Tracker AD wrapper for ", ∇ℓ.ℓ)

function logdensity_and_gradient(∇ℓ::TrackerGradientLogDensity, x::AbstractVector{T}) where {T}
    @unpack ℓ = ∇ℓ
    y, back = Tracker.forward(x -> logdensity(ℓ, x), x)
    yval = Tracker.data(y)
    # work around https://github.com/FluxML/Flux.jl/issues/497
    z = T <: Real ? zero(T) : 0.0
    S = typeof(z + 0.0)
    S(yval)::S, (S.(first(Tracker.data.(back(1)))))::Vector{S}
end


text:
""" Tonti Diagrams

This file includes a framework for constructing and evaluating Tonti diagrams.
Tonti diagrams are stored as ACSets, and have an imperative interface for
describing physical variables and the relationships between them. This tooling
also lets a Tonti diagram be converted to a vectorfield, allowing for
simulation of physical systems through any traditional vectorfield solver.
"""

module TontiDiagrams

using Catlab.Present
using Catlab.Theories
using Catlab.CategoricalAlgebra
using Catlab.CategoricalAlgebra.FinSets
using Catlab.Graphs
import Catlab.Graphs: Graph
using Catlab.Graphics
using Catlab.Programs
using Catlab.WiringDiagrams
using CombinatorialSpaces
using CombinatorialSpaces: ⋆

export TheoryTontiDiagram, TontiDiagram, Space,
        add_variable!, add_variables!,
        add_derivative!, add_derivatives!,
        add_time_dep!,
        add_laplacian!, add_transition!, add_bc!,
        vectorfield, Open, gen_form

""" ACSet definition for a Tonti diagram.

This diagram is visualized in the q.uiver framework
[here](https://q.uiver.app/?q=WzAsMTAsWzEsMiwiTyJdLFsyLDIsIkkiXSxbMSwxLCJWIl0sWzEsMywiVCJdLFswLDIsIkJDIl0sWzAsMCwiVEQiXSxbMCwzLCJGdW5jIl0sWzMsMCwiQ29tcCJdLFszLDEsIkRpbWVuc2lvbiJdLFszLDIsIkxhYmVsIl0sWzQsMiwiYmN2Il0sWzUsMiwiaW50ZWciLDAseyJvZmZzZXQiOi0xfV0sWzAsMiwib3YiLDFdLFswLDMsIm90IiwxXSxbMSwyLCJpdiIsMV0sWzEsMywiaXQiLDFdLFs1LDIsImRlcml2IiwyLHsib2Zmc2V0IjoxfV0sWzQsNiwiYmNmdW5jIiwyLHsic3R5bGUiOnsiYm9keSI6eyJuYW1lIjoiZGFzaGVkIn19fV0sWzMsNiwidGZ1bmMiLDAseyJzdHlsZSI6eyJib2R5Ijp7Im5hbWUiOiJkYXNoZWQifX19XSxbMiw5LCJzeW1ib2wiLDEseyJzdHlsZSI6eyJib2R5Ijp7Im5hbWUiOiJkYXNoZWQifX19XSxbMiw4LCJkaW1lbnNpb24iLDEseyJzdHlsZSI6eyJib2R5Ijp7Im5hbWUiOiJkYXNoZWQifX19XSxbMiw3LCJjb21wbGV4IiwxLHsic3R5bGUiOnsiYm9keSI6eyJuYW1lIjoiZGFzaGVkIn19fV1d).
This Tonti diagram definition is currently very close to the simulation level
and generalized away from the DEC tooling. Future work will be done to provide
a closer integration between Tonti diagrams and DEC.

See Catlab.jl documentation for a description of the @present syntax.
"""
@present TheoryTontiDiagram(FreeSchema) begin
  Func::Data
  Comp::Data
  Dimension::Data
  Label::Data

  V::Ob
  I::Ob
  O::Ob
  T::Ob
  BC::Ob
  TD::Ob

  iv::Hom(I,V)
  it::Hom(I,T)
  ov::Hom(O,V)
  ot::Hom(O,T)
  bcv::Hom(BC,V)
  deriv::Hom(TD,V)
  integ::Hom(TD,V)

  tfunc::Attr(T,Func)
  bcfunc::Attr(BC,Func)

  complex::Attr(V,Comp)
  dimension::Attr(V,Dimension)
  symbol::Attr(V, Label)
  tlabel::Attr(T, Label)
  bclabel::Attr(BC, Label)
end

""" Space

Structure for storing the computed values for the primal/dual complex, the
boundary operators, the hodge stars, the laplacian operators, and the lie
operators. This is the operator which provides the DEC connection for the Tonti
diagra tooling.

Accessing the boundary, hodge, and laplacian operators is slightly complicated
by using arrays in a 1-indexed system. The hodge and boundary operators are
accessed with two indices, where the first determines the complex (1 for primal
and 2 for dual) and the second determines the dimension (1 for 0-forms, 2 for
1-forms, 3 for 2-forms, etc.). Thus, the boundary operator from dual 1-forms to
dual 2-forms is given as:

```julia
sp = Space(s)
sp.boundary[2,2]
```

The other goal of this structure is to cache the computed operators, yet
optimally this caching should be perfomed during the `vectorfield` operation.
"""
struct Space
  s::EmbeddedDeltaSet2D
  sd::EmbeddedDeltaDualComplex2D
  boundary
  hodge
  laplacian
  lie
end

""" Space(s::EmbeddedDeltaSet3D)

Caclulates all of the values stored in the `Space` object for a given complex
`s`.
"""
function Space(s::EmbeddedDeltaSet2D{O, P}) where {O, P}
  sd = EmbeddedDeltaDualComplex2D{O, eltype(P), P}(s)
  subdivide_duals!(sd, Barycenter())

  boundary = Array{typeof(d(0,s)), 2}(undef, 2,2)
  boundary[[1,2,3,4]] .= [d(0,s), dual_derivative(0,sd), d(1,s), dual_derivative(1,sd)]

  hodge = Array{typeof(⋆(0,sd)), 2}(undef, 2,3)
  hodge[[1,2,3,4,5,6]] .= [⋆(0,sd), inv(⋆(2,sd)), ⋆(1,sd), -1 .* inv(⋆(1,sd)), ⋆(2,sd), inv(⋆(0,sd))]

  # TODO:
  # Add Laplacian and Lie operators here
  laplacian  = Array{typeof(hodge[1,2]*boundary[1,1]),1}(undef,3)
  laplacian[1] = hodge[2,3]*boundary[2,2]*hodge[1,2]*boundary[1,1]

  laplacian[2] = hodge[2,2]*boundary[2,1]*hodge[1,3]*boundary[1,2] .+
                 boundary[1,1]*hodge[2,3]*boundary[2,2]*hodge[1,2]
  laplacian[3] = boundary[1,2]*hodge[2,2]*boundary[2,1]*hodge[1,3]

  lie_op = [hodge[2,3]*boundary[2,2]*hodge[1,2], hodge[2,2]*boundary[2,1]*hodge[1,3]]
  lie = [(∂s,s,v)->(∂s .= boundary[1,2]*∧(Tuple{0,1},sd,s,bound[2,1]*v)),
         (∂u,u,v)->(∂u .= (star_1*∧(Tuple{1,0},sd,v,inv_star_0*cobound_1_2*u) .+ cobound_0_1*star_2*∧(Tuple{1,1},sd,v,inv_star_1*u))),

        ]
  Space(s, sd, boundary, hodge, laplacian, nothing)
end

# TODO:
# Make more petri-style simulation ACSet
#
# Make more Tonti-aware ACSet w/ support of multiple complexes

# These functions are necessary for defining a TontiDiagram data structure from
# the ACSet
const AbstractTontiDiagram = AbstractACSetType(TheoryTontiDiagram)
const TontiDiagram = ACSetType(TheoryTontiDiagram,index=[:iv, :it, :ov, :ot, :bcv],
                                                 unique_index=[:symbol])

# Define an interface for an OpenTontiDiagram which allows Tonti diagrams to
# be composed over their corners
const OpenTontiDiagramOb, OpenTontiDiagram = OpenACSetTypes(TontiDiagram,
                                                            :V)

""" Open(td::TontiDiagram, states::Symbol)

Generates an OpenTontiDiagram with cospan legs on variables defined by the
symbols included in `states`. This OpenTontiDiagram can then be composed with
other OpenTontiDiagrams over a pattern given by an undirected wiring diagram.

```julia
OpenTontiDiagram(td, :x, :v)
```
"""
Open(td, states...) = OpenTontiDiagram{Function, Bool, Int64, Symbol}(td,
                        map(v->FinFunction([incident(td, v, :symbol)], nparts(td,:V)), states)...)

""" TontiDiagram()

Initialize an empty TontiDiagram object.
"""
TontiDiagram() = TontiDiagram{Function, Bool, Int64, Symbol}()

""" vectorfield(td::AbstractTontiDiagram, sp::Space)

Generates a Julia function which calculates the vectorfield of the Tonti
diagram. The state of the system is defined by a single vector which is a
flattening of all state variables of the system. Thus, this function returns
both the indices of each variable in the state-vector along with the
vectorfield function itself.

The resulting function has a signature of the form `f!(du, u, p, t)` and can be
passed to the DifferentialEquations.jl solver package.
"""
function vectorfield(td, sp::Space)
  fv_mem, ft_mem = init_mem(td, sp)
  dg = Graph(td)

  order = topological_sort(dg)

  state_vars = filter(x->length(incident(td,x,:ov))==0, 1:nparts(td, :V))

  input_vars = Dict{Symbol, Tuple{Int64,Int64}}()
  cur_head = 1
  for i in state_vars
    v_size = length(fv_mem[i])
    input_vars[td[i,:symbol]] = (cur_head,cur_head+v_size-1)
    cur_head += v_size
  end

  data_buffer = Dict{Type, Tuple{Array, Array}}(Float64=>(fv_mem, ft_mem))
  function system(du, u, t, p)
    if !(eltype(u) ∈ keys(data_buffer))
      data_buffer[eltype(u)] = init_mem(td, sp, type=eltype(u))
    end

    v_mem, t_mem = data_buffer[eltype(u)]

    for cur in order
      # Check if current is a variable or transition
      if cur > nparts(td, :V)
        cur -= nparts(td, :V)
        inputs = td[incident(td, cur, :it), :iv]
        outputs = incident(td, cur, :ot)
        td[cur,:tfunc](t_mem[outputs]..., v_mem[inputs]...)
      else
        inputs = incident(td, cur, :ov)
        if(length(inputs) == 0)
          # This means this is a state variable
          data_source = input_vars[td[cur,:symbol]]
          v_mem[cur] .= u[data_source[1]:data_source[2]]
        else
          v_mem[cur] .= 0
          for i in inputs
            v_mem[cur] .+= t_mem[i]
          end
        end
        bcs = incident(td, cur, :bcv)
        for bc in bcs
          td[bc, :bcfunc](v_mem[cur])
        end
      end
    end
    # If a state variable does not have a derivative defined, we keep it out
    # (we'll want to move these to the parameter argument instead)
    du .= 0
    for i in state_vars
      state_range = input_vars[td[i,:symbol]]
      out_var = td[incident(td, i, :integ), :deriv]
      if length(out_var) != 0
        du[state_range[1]:state_range[2]] .= v_mem[out_var[1]]
      end
    end
  end
  input_vars, system
end

""" add_variable!(td:TontiDiagram, symbol::Symbol, dimension::Int64, complex::Bool)

Adds a variable to the TontiDiagram system which can later be referenced by its
`symbol`. This constructor requires the dimensionality of the variable (0 ->
point, 1 -> line, etc.) and the complex it is defined on
(true -> primal/straight, false -> dual/twisted).

Definiting a system with a variable `v` defined on the primal lines would be
constructed as:
```@example
td = TontiDiagram()
add_variable!(td, :v, 1, true)
```
"""
function add_variable!(td, symbol::Symbol, dimension::Int64, complex::Bool)
  add_part!(td, :V,symbol=symbol, dimension=dimension, complex=complex)
end

""" add_variables!(td:TontiDiagram, vars::Tuple{Symbol, Int64, Bool}...)

Adds multiple variables to the TontiDiagram system which can later be
referenced by their symbols. This constructor follows the same pattern as
`add_variable!` with each variable specified as a tuple of:
```julia
(symbol, dimension, complex)
```

Defining a system with a variable `v` defined on the primal lines, a variable
`p̃` defined on the dual surfaces, and a variable `C` defined on the primal
surfaces would be constructed as:
```@example
td = TontiDiagram()
add_variables!(td, (:v, 1, true), (:p̃, 1, true), (:C, 2, true))
```
"""
function add_variables!(td, vars::Tuple{Symbol, Int64, Bool}...)
  for v in vars
    add_variable!(td, v...)
  end
end

""" add_transition!(td, dom_sym::Vector{Symbol}, func!, codom_sym::Vector{Symbol})

Adds a transition function from variables `dom_sym` to variables `codom_sym`
with its transition defined by `func!`. `func!` is expected to have the signature
`func(codom_sym..., dom_sym...)` and is expected to modify the values of the
`codom_sym` variables.

Defining a transition from variables `x` and `y` to `m` that calculates the
magnitude of the values in `x` and `y` as vector coordinates would be defined as:

```@example
td = TontiDiagram()
add_variables!(td, (:x, 0, 1), (:y, 0, 1). (:m, 0, 1))
add_transition!(td, [:x, :y], (m,x,y)->(m .= sqrt.(x .* y)), [:m])
```
"""
function add_transition!(td, dom_sym::Array{Symbol,1}, func!, codom_sym::Array{Symbol,1}; label=Symbol(""))
  dom = [findfirst(v->v == s, td[:symbol]) for s in dom_sym]
  codom = [findfirst(v->v == s, td[:symbol]) for s in codom_sym]
  t = add_part!(td, :T, tfunc=func!, tlabel=label)
  add_parts!(td, :I, length(dom), iv=dom, it=t)
  add_parts!(td, :O, length(codom), ov=codom, ot=t)
end

""" add_derivative!(td::TontiDiagram, sp::Space, dom_sym::Symbol,
codom_sym::Symbol)

Adds a derivative transition from variable `dom_sym` to variable `codom_sym`
using the boundary operators from `sp`. This function determines which boundary
operator to use and inserts an appropriate transition between the two
variables.

Defining a spatial derivative relationship between the primal 0-form `x` and
the primal 1-form `Δx` cis given as follows:
```julia
add_derivative(td, sp, :x, :Δx)
```
"""
function add_derivative!(td, sp, dom_sym, codom_sym)
  dom = findfirst(v->v==dom_sym, td[:symbol])
  codom = findfirst(v->v==codom_sym, td[:symbol])

  # TODO:
  # Add tests for proper dimensions, complexes, etc.
  # This will later be replaced as we pre-initialize all boundary operators
  bound = sp.boundary[(td[dom,:complex] ? 1 : 2), td[dom,:dimension]+1]
  func(x,y) = (x.=bound*y)
  add_transition!(td, [dom_sym],func,[codom_sym]; label=:d)
end

""" add_derivatives!(td::TontiDiagram, sp::Space, vars:Pair{Symbol, Symbol}...)

Adds multiple derivative transition between pairs of variables, using the same
syntax as in `add_derivative!`.

Example usage:
```julia
add_derivatives!(td, sp, (:x,:Δx), (:y, :Δy))
```
"""
function add_derivatives!(td, sp, vars::Pair{Symbol, Symbol}...)
  for v in vars
    add_derivative!(td, sp, v[1],v[2])
  end
end

""" add_time_dep!(td::TontiDiagram, deriv_sym::Symbol, integ_sym::Symbol)

Adds a time derivative relationship between the variables `deriv_sym` and
`integ_sym` (where `deriv_sym` is the time derivative of `integ_sym`). These
relationships are used to determine the state-variables of the system.
"""
function add_time_dep!(td, deriv_sym::Symbol, integ_sym::Symbol)
  deriv = findfirst(v->v==deriv_sym, td[:symbol])
  integ = findfirst(v->v==integ_sym, td[:symbol])

  add_part!(td, :TD, integ=integ, deriv=deriv)
end

""" add_bc!(td::TontiDiagram, var_sys::Symbol, func::Function)

Adds a "boundary condition" to the variable `var_sys` by applying `func` to the
values of this variable during simulation. This function is the last one
evaluated on the data of the variable `var_sys`, and so can be used to enforce
any relevant boundary conditions.

TODO:
Add time dependency of boundary condition function to allow for time-varying BCs.
"""
function add_bc!(td, var_sym, func; label=:BC)
  var = findfirst(v->v==var_sym, td[:symbol])
  add_part!(td, :BC, bcfunc=func, bcv=var, bclabel=label)
end

# Note: This function can be made more efficient if combined with existing
# transformations.
# e.g. Advection-diffusion can be merged after the initial wedge
# product/coboundary operator
#
# Currently only defined on primal complices (can this be applied to dual
# complices?)

""" add_laplacian!(td::TontiDiagrams, sp::Space, dom_sym::Symbol, codom_sym::Symbol; coef::Float64)

Adds a transition which defines `codom_sym` as the laplacian of `dom_sym` with
a constant scaling factor of `coef`.
"""
function add_laplacian!(td, sp, dom_sym, codom_sym; coef=1.0)
  sd = sp.sd
  dom = findfirst(v->v==dom_sym, td[:symbol])
  codom = findfirst(v->v==codom_sym, td[:symbol])

  lap_op = sp.laplacian[td[dom,:dimension]+1] # laplace_beltrami(Val{td[dom,:dimension]},sd)
  func(x,y) = (x .= coef * (lap_op*y))
  add_transition!(td, [dom_sym], func, [codom_sym]; label=:Δ)
end

function init_mem(td, s::EmbeddedDeltaSet1D)
  # Fill out this function
end

function init_mem(td, sp::Space; type=Float64)
  s = sp.s
  primal_size = [nv(s), ne(s), ntriangles(s)]
  dual_size   = [ntriangles(s), ne(s), nv(s)]

  t_mem = Array{Array{type,1},1}()
  v_mem = Array{Array{type,1},1}()

  for i in 1:nparts(td, :O)
    var = td[i,:ov]
    push!(t_mem, zeros(type,
                      td[var,:complex] ? primal_size[td[var,:dimension]+1] :
                      dual_size[td[var,:dimension]+1]))
  end

  for v in 1:nparts(td,:V)
    push!(v_mem, zeros(type,
                      td[v,:complex] ? primal_size[td[v,:dimension]+1] :
                      dual_size[td[v,:dimension]+1]))
  end

  v_mem, t_mem
end

function Graph(td; extended=false)
  g = Graph()
  add_vertices!(g, nparts(td, :V) + nparts(td, :T))
  nvars = nparts(td, :V)
  for i in 1:nparts(td, :I)
    add_edge!(g, td[i,:iv], td[i,:it] + nvars)
  end
  for o in 1:nparts(td, :O)
    add_edge!(g, td[o,:ot] + nvars, td[o,:ov])
  end
  if extended
    for bc in 1:nparts(td, :BC)
      add_edge!(g, td[bc, :bcv], td[bc, :bcv])
    end
    for t in 1:nparts(td, :TD)
      add_edge!(g, td[t, :deriv], td[t, :integ])
    end
  end
  g
end

""" Construct a 0-form based on a scalar function

This operator accepts a scalar function and evaulates it at each point on the
simplex, returning a 0-form.
"""
function gen_form(s::EmbeddedDeltaSet2D, f::Function)
  map(f, point(s))
end

function propertygraph(td::TontiDiagram;
                       prog::AbstractString="neato", graph_attrs::AbstractDict=Dict(),
                       node_attrs::AbstractDict=Dict(), edge_attrs::AbstractDict=Dict(:len=>"1.5"),
                       node_labels::Bool=true, edge_labels::Bool=false)

  nvars = nparts(td, :V)
  isVar(v) = (v <= nvars)
  g = Graph(td, extended = true)

  dim_colors = ["green", "yellow", "red", "blue"]

	node_labeler(v) = begin
    if isVar(v)
      Dict(:label=>"$(td[v, :symbol])",
           :style=>td[v,:complex] ? "filled" : "filled, dotted",
           :width=>"0.75", :height=>"0.75", :fixedsize=>"false",
           :shape=>"circle",
           :fillcolor=>dim_colors[td[v,:dimension]+1],
           )
    else
      Dict(:label=>"$(td[v - nvars, :tlabel])",
           :width=>"0.75", :height=>"0.75", :fixedsize=>"false",
           :shape=>"square"
           )
    end
  end

  edge_labeler(e) = begin
    if e <= nparts(td, :I) + nparts(td, :O)
      Dict(:color=>"black")
    else
      e -= nparts(td, :I) + nparts(td, :O)
      if e <= nparts(td, :BC)
        Dict(:color=>"black",
             :label=>"$(td[e, :bclabel])")
      else
        e -= nparts(td, :BC)
        Dict(:color=>"black", :style=>"dotted")
      end
    end
  end

  PropertyGraph{Any}(g, node_labeler, edge_labeler;
                        prog = prog,
                        graph = merge!(Dict(:rankdir => "TB"), graph_attrs),
                        node = merge!(Graphics.GraphvizGraphs.default_node_attrs(node_labels), node_attrs),
                        edge = merge!(Dict(:arrowsize => "0.5"), edge_attrs),
                        )
end

function Graphics.to_graphviz(td::TontiDiagram; kwargs...)
  propertygraph(td; kwargs...) |> to_graphviz
end

end


text:
using AbstractOperators, RecursiveArrayTools
using ProximalOperators

A1 = [ 1.0  -2.0   3.0;
       2.0  -1.0   0.0;
      -1.0   0.0   4.0;
      -1.0  -1.0  -1.0]
A2 = [-4.0  5.0;
	  -1.0  3.0;
	  -3.0  2.0;
	   1.0  3.0]
A = [A1 A2]

opA1 = MatrixOp(A1)
opA2 = MatrixOp(A2)
opA = HCAT(opA1, opA2)

b = [1.0, 2.0, 3.0, 4.0]

m, n = size(A)

f = Translate(SqrNormL2(), -b)
lam = 0.1*norm(A'*b, Inf)
g = SeparableSum(NormL1(lam), NormL1(lam))

x_star = ArrayPartition([-3.877278911564627e-01, 0, 0], [2.174149659863943e-02, 6.168435374149660e-01])

# Nonfast/Nonadaptive

x0 = zero(x_star)
@time it, x, sol = ProximalAlgorithms.FBS(x0; fq=f, Aq=opA, g=g, gamma=1.0/opnorm(A)^2)
@test maximum(abs,x .- x_star) <= 1e-4
@test it < 150
println(sol)

# Nonfast/Adaptive

x0 = zero(x_star)
@time it, x, sol = ProximalAlgorithms.FBS(x0; fq=f, Aq=opA, g=g, adaptive=true)
@test maximum(abs,x .- x_star) <= 1e-4
@test it < 300
println(sol)

# Fast/Nonadaptive

x0 = zero(x_star)
@time it, x, sol = ProximalAlgorithms.FBS(x0; fq=f, Aq=opA, g=g, gamma=1.0/opnorm(A)^2, fast=true)
@test maximum(abs, x .- x_star) <= 1e-4
@test it < 100
println(sol)

# Fast/Adaptive

x0 = zero(x_star)
@time it, x, sol = ProximalAlgorithms.FBS(x0; fq=f, Aq=opA, g=g, adaptive=true, fast=true)
@test maximum(abs,x .- x_star) <= 1e-4
@test it < 200
println(sol)

# ZeroFPR/Adaptive

x0 = zero(x_star)
@time it, x, sol = ProximalAlgorithms.ZeroFPR(x0; fq=f, Aq=opA, g=g, adaptive=true)
@test maximum(abs,x .- x_star) <= 1e-4
@test it < 15
println(sol)

# PANOC/Adaptive

x0 = zero(x_star)
@time it, x, sol = ProximalAlgorithms.PANOC(x0; fq=f, Aq=opA, g=g, adaptive=true)
@test maximum(abs,x .- x_star) <= 1e-4
@test it < 20
println(sol)


text:
using ArgParse
using Random
using JLD
using NPZ
using UUIDs

include("./FESOM.jl")

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table s begin
        "--temperature", "-T"
            help = "Temperature of system."
            arg_type = Float64
            default = 0.0
        "--number_of_iterations", "-N"
            help = "Number of generations before genetic algorithm quits."
            arg_type = Int64
            default = 10000
        "--stop_minimum_fitness"
            help = "Regularization constant."
            arg_type = Float64
            default = 1.0e-8
        "--save_file_dir"
            help = "Directory to save results in."
            arg_type = String
            default = "./fesomresults"
        "--seed"
            help = "Seed to pass to random number generator."
            arg_type = Int64
            default = 1
        "--number_of_blas_threads"
            help = "Number of BLAS threads to utilize if > 0 otherwise uses Julia default value."
            arg_type = Int64
            default = 0
        "--track_stats"
            help = "Track minimum fitness."
            action = :store_true
        "qmc_data"
            help = "*.npz or *.jld file containing quantum Monte Carlo data with columns: IMAGINARY_TIME, INTERMEDIATE_SCATTERING_FUNCTION, ERROR"
            arg_type = String
            required = true
        "initial_dsf"
            help = "*.npz or *.jld file containing initial guess for dynamic structure factor and frequency with columns: DSF, FREQUENCY"
            arg_type = String
            required = true
    end

    return parse_args(s)
end

function main()
    start = time();
    parsed_args = parse_commandline()

    # create data directory
    try
        mkpath(parsed_args["save_file_dir"]);
    catch
        nothing
    end
    save_dir = parsed_args["save_file_dir"];

    _ext = splitext(parsed_args["qmc_data"])[2];
    if _ext == ".npz"
        qmcdata=NPZ.npzread(parsed_args["qmc_data"]);
    elseif _ext == ".jld"
        qmcdata = load(parsed_args["qmc_data"]);
    else
        throw(AssertionError("qmc_data must be *.jld or *.npz"));
    end
    imaginary_time = qmcdata["tau"];
    isf = qmcdata["isf"];
    isf_error = qmcdata["error"];

    _ext = splitext(parsed_args["initial_dsf"])[2];
    if _ext == ".npz"
        dsfdata=NPZ.npzread(parsed_args["initial_dsf"]);
    elseif _ext == ".jld"
        dsfdata = load(parsed_args["initial_dsf"]);
    else
        throw(AssertionError("initial_dsf must be *.jld or *.npz"));
    end
    dsf = dsfdata["dsf"];
    frequency_bins = dsfdata["frequency"];

    initial_dsf = copy(dsf);

    u4 = uuid4();
    rng = MersenneTwister(parsed_args["seed"]);
    beta = 1/parsed_args["temperature"];
    dsf .*= (1 .+ exp.(-beta .* frequency_bins));
    rng,dsf_result,fitness, minimum_fitness = FESOM.fesom(rng,dsf,isf,isf_error,
                                                          frequency_bins,imaginary_time,
                                                          temperature = parsed_args["temperature"],
                                                          number_of_iterations = parsed_args["number_of_iterations"],
                                                          stop_minimum_fitness = parsed_args["stop_minimum_fitness"],
                                                          track_stats = parsed_args["track_stats"],
                                                          number_of_blas_threads = parsed_args["number_of_blas_threads"])
    elapsed = time() - start;
    seed = lpad(parsed_args["seed"],4,'0');
    filename = "$(save_dir)/fesom_results_$(seed)_$u4.jld";
    println("Saving results to $filename");
    save(filename,
         "u4",u4,
         "dsf",dsf_result,
         "frequency",frequency_bins,
         "fitness",fitness,
         "elapsed_time",elapsed);
    filename = "$(save_dir)/fesom_stats_$(seed)_$u4.jld";

    if parsed_args["track_stats"]
        println("Saving stats to $filename");
        save(filename,
             "u4",u4,
             "minimum_fitness",minimum_fitness);
    end

    filename = "$(save_dir)/fesom_params_$(seed)_$u4.jld";
    println("Saving parameters to $filename");
    save(filename,
         "u4",u4,
         "rng",rng,
         "imaginary_time",imaginary_time,
         "isf",isf,
         "isf_error",isf_error,
         "dsf",dsf_result,
         "frequency",frequency_bins,
         "initial_dsf",initial_dsf,
         "number_of_iterations",parsed_args["number_of_iterations"],
         "temperature",parsed_args["temperature"],
         "stop_minimum_fitness",parsed_args["stop_minimum_fitness"],
         "seed",parsed_args["seed"])
    nothing
end

main()


text:
############################
# Circuits
############################

type Circuit{N} <: AbstractQuCircuit{N}
    gates::Array{GateUnit,1}
end

Circuit(num::Integer,gates::Array{GateUnit,1}) = Circuit{num}(gates)
Circuit(num::Integer)=Circuit{num}(Array(GateUnit,0))

# TODO
# function show(io::IO,circuit::Circuit)
# end

############################
# Circuit Constructor
############################

max(a::Int) = a

# NOTE
# Should we open AbstractGateUnit as API for users?

function addgate!{N,M}(circuit::AbstractQuCircuit{N},gate::AbstractGateUnit{M})
  push!(circuit.gates,gate)
  sort!(circuit,gates,alg=QuickSort,by=x->x.pos[1])
end

# For gate units
function addgate!{T,N,M}(circuit::AbstractQuCircuit{N},gate::Gate{T,M},pos::Vector{Int})
    # Bounds check
    @assert length(pos)==M+1 "number of qubits do not match"
    @assert maximum(pos[2:end])<=N "bit's id out of range, maximum is $N"

    addgate!(circuit,GateUnit(gate,pos))
end

addgate!{T,N,M}(circuit::AbstractQuCircuit{N},gate::Gate{T,M},pos::Int...) = addgate!(circuit,gate,collect(pos))

# for controled gates
function addgate!{T,N,M}(circuit::AbstractQuCircuit{N},gate::Gate{T,M},pos::Vector{Int},ctrl::Vector{Int})
    # Bounds check
    @assert length(pos)+length(ctrl) == M+1
    @assert max(maximum(pos[2:end]),maximum(ctrl)) <= N "bit's id out of range, maximum is $N"

    addgate!(circuit,CtrlGateUnit(gate,pos,ctrl))
end

function removegate!{N}(circuit::AbstractQuCircuit{N},pos::Vector{Int})
    for i in eachindnex(circuits.gates)
        if circuits.gates.pos == pos
            deleteat!(circuits.gates,i)
        end
    end
end


text:
__precompile__()

module Netpbm

using Images, FileIO, ColorTypes, FixedPointNumbers
typealias AbstractGray{T} Color{T, 1}

import FileIO: load, save

const is_little_endian = ENDIAN_BOM == 0x04030201
const ufixedtype = Dict(10=>UFixed10, 12=>UFixed12, 14=>UFixed14, 16=>UFixed16)

function load(f::Union{File{format"PBMBinary"},File{format"PGMBinary"},File{format"PPMBinary"}})
    open(f) do s
        skipmagic(s)
        load(s)
    end
end

function load(s::Stream{format"PBMBinary"})
    io = stream(s)
    w, h = parse_netpbm_size(io)
    dat = BitArray(w, h)
    nbytes_per_row = ceil(Int, w/8)
    for irow = 1:h, j = 1:nbytes_per_row
        tmp = read(io, UInt8)
        offset = (j-1)*8
        for k = 1:min(8, w-offset)
            dat[offset+k, irow] = (tmp>>>(8-k))&0x01
        end
    end
    Image(dat, Dict("spatialorder" => ["x", "y"], "pixelspacing" => [1,1]))
end

function load(s::Stream{format"PGMBinary"})
    io = stream(s)
    w, h = parse_netpbm_size(io)
    maxval = parse_netpbm_maxval(io)
    local dat
    if maxval <= 255
        dat = reinterpret(Gray{U8}, read(io, U8, w, h), (w, h))
    elseif maxval <= typemax(UInt16)
        datraw = Array(UInt16, w, h)
        if !is_little_endian
            for indx = 1:w*h
                datraw[indx] = read(io, UInt16)
            end
        else
            for indx = 1:w*h
                datraw[indx] = bswap(read(io, UInt16))
            end
        end
        # Determine the appropriate UFixed type
        T = ufixedtype[ceil(Int, log2(maxval)/2)<<1]
        dat = reinterpret(Gray{T}, datraw, (w, h))
    else
        error("Image file may be corrupt. Are there really more than 16 bits in this image?")
    end
    T = eltype(dat)
    Image(dat, Dict("colorspace" => "Gray", "spatialorder" => ["x", "y"], "pixelspacing" => [1,1]))
end

function load(s::Stream{format"PPMBinary"})
    io = stream(s)
    w, h = parse_netpbm_size(io)
    maxval = parse_netpbm_maxval(io)
    local dat
    if maxval <= 255
        datraw = read(io, U8, 3, w, h)
        dat = reinterpret(RGB{U8}, datraw, (w, h))
    elseif maxval <= typemax(UInt16)
        # read first as UInt16 so the loop is type-stable, then convert to UFixed
        datraw = Array(UInt16, 3, w, h)
        # there is no endian standard, but netpbm is big-endian
        if !is_little_endian
            for indx = 1:3*w*h
                datraw[indx] = read(io, UInt16)
            end
        else
            for indx = 1:3*w*h
                datraw[indx] = bswap(read(io, UInt16))
            end
        end
        # Determine the appropriate UFixed type
        T = ufixedtype[ceil(Int, log2(maxval)/2)<<1]
        dat = reinterpret(RGB{T}, datraw, (w, h))
    else
        error("Image file may be corrupt. Are there really more than 16 bits in this image?")
    end
    T = eltype(dat)
    Image(dat, Dict("spatialorder" => ["x", "y"], "pixelspacing" => [1,1]))
end

function save(filename::File{format"PGMBinary"}, img; kwargs...)
    open(filename, "w") do s
        io = stream(s)
        write(io, "P5\n")
        write(io, "# pgm file written by Julia\n")
        save(s, img; kwargs...)
    end
end

function save(filename::File{format"PPMBinary"}, img; kwargs...)
    open(filename, "w") do s
        io = stream(s)
        write(io, "P6\n")
        write(io, "# ppm file written by Julia\n")
        save(s, img; kwargs...)
    end
end

pnmmax{T<:AbstractFloat}(::Type{T}) = 255
pnmmax{T<:UFixed}(::Type{T}) = reinterpret(FixedPointNumbers.rawtype(T), one(T))
pnmmax{T<:Unsigned}(::Type{T}) = typemax(T)

function save{T<:Gray}(s::Stream{format"PGMBinary"}, img::AbstractArray{T}; mapi = mapinfo(img))
    w, h = widthheight(img)
    TE = eltype(T)
    mx = pnmmax(TE)
    write(s, "$w $h\n$mx\n")
    p = permutation_horizontal(img)
    writepermuted(s, img, mapi, p)
end

function save{T<:Number}(s::Stream{format"PGMBinary"}, img::AbstractArray{T}; mapi = mapinfo(img))
    io = stream(s)
    w, h = widthheight(img)
    cs = colorspace(img)
    cs == "Gray" || error("colorspace $cs not supported")
    mx = pnmmax(T)
    write(io, "$w $h\n$mx\n")
    p = permutation_horizontal(img)
    writepermuted(io, img, mapi, p)
end

function save{T<:Color}(s::Stream{format"PPMBinary"}, img::AbstractArray{T}; mapi = mapinfo(img))
    w, h = widthheight(img)
    TE = eltype(eltype(mapi))
    mx = pnmmax(TE)
    write(s, "$w $h\n$mx\n")
    p = permutation_horizontal(img)
    writepermuted(s, img, mapi, p; gray2color = T <: AbstractGray)
end

function save{T}(s::Stream{format"PPMBinary"}, img::AbstractArray{T}; mapi = mapinfo(img))
    io = stream(s)
    w, h = widthheight(img)
    cs = colorspace(img)
    in(cs, ("RGB", "Gray")) || error("colorspace $cs not supported")
    mx = pnmmax(T)
    write(io, "$w $h\n$mx\n")
    p = permutation_horizontal(img)
    writepermuted(io, img, mapi, p; gray2color = cs == "Gray")
end

# Permute to a color, horizontal, vertical, ... storage order (with time always last)
function permutation_horizontal(img)
    cd = colordim(img)
    td = timedim(img)
    p = spatialpermutation(["x", "y"], img)
    if cd != 0
        p[p .>= cd] += 1
        insert!(p, 1, cd)
    end
    if td != 0
        push!(p, td)
    end
    p
end

permutedims_horizontal(img) = permutedims(img, permutation_horizontal(img))

# Write values in permuted order
let method_cache = Dict()
global writepermuted
# Delete the following once img[i,j] returns the Color for an ImageCmap
writepermuted(stream, img::ImageCmap, mapi::MapInfo, perm; gray2color::Bool = false) =
    writepermuted(stream, convert(Image, img), mapi, perm; gray2color=gray2color)

function writepermuted(stream, img, mapi::MapInfo, perm; gray2color::Bool = false)
    cd = colordim(img)
    key = (perm, cd, gray2color)
    if !haskey(method_cache, key)
        swapfunc = is_little_endian ? :mybswap : :identity
        loopsyms = [symbol(string("i_",d)) for d = 1:ndims(img)]
        body = gray2color ? quote
                g = $swapfunc(map(mapi, img[$(loopsyms...)]))
                write(stream, g)
                write(stream, g)
                write(stream, g)
            end : quote
                write(stream, $swapfunc(map(mapi, img[$(loopsyms...)])))
            end
        loopargs = [:($(loopsyms[d]) = 1:size(img, $d)) for d = 1:ndims(img)]
        loopexpr = Expr(:for, Expr(:block, loopargs[perm[end:-1:1]]...), body)
        f = @eval begin
            local _writefunc_
            function _writefunc_(stream, img, mapi)
                $loopexpr
            end
        end
    else
        f = method_cache[key]
    end
    f(stream, img, mapi)
    nothing
end
end

function parse_netpbm_size(stream::IO)
    szline = strip(readline(stream))
    while isempty(szline) || szline[1] == '#'
        szline = strip(readline(stream))
    end
    parseints(szline, 2)
end

function parse_netpbm_maxval(stream::IO)
    skipchars(stream, isspace, linecomment='#')
    maxvalline = strip(readline(stream))
    parse(Int, maxvalline)
end

function parseints(line, n)
    ret = Array(Int, n)
    pos = 1
    for i = 1:n
        pos2 = search(line, ' ', pos)
        if pos2 == 0
            pos2 = length(line)+1
        end
        ret[i] = parse(Int, line[pos:pos2-1])
        pos = pos2+1
        if pos > length(line) && i < n
            error("Line terminated without finding all ", n, " integers")
        end
    end
    tuple(ret...)

end

function Base.write{T<:UFixed}(io::IO, c::AbstractRGB{T})
    write(io, reinterpret(red(c)))
    write(io, reinterpret(green(c)))
    write(io, reinterpret(blue(c)))
end

function Base.write(io::IO, c::RGB24)
    write(io, red(c))
    write(io, green(c))
    write(io, blue(c))
end

Base.write(io::IO, c::Gray) = write(io, reinterpret(gray(c)))
Base.write(io::IO, c::UFixed) = write(io, reinterpret(c))

mybswap(i::Integer) = bswap(i)
mybswap(i::UFixed) = bswap(reinterpret(i))
mybswap(c::RGB24) = c
mybswap{T<:UFixed}(c::AbstractRGB{T}) = RGB{T}(T(bswap(reinterpret(red(c))),0),
                                               T(bswap(reinterpret(green(c))),0),
                                               T(bswap(reinterpret(blue(c))),0))

# Netpbm mapinfo client. Converts to RGB and uses UFixed.
mapinfo{T<:Unsigned}(img::AbstractArray{T}) = MapNone{T}()
mapinfo{T<:UFixed}(img::AbstractArray{T}) = MapNone{T}()
mapinfo{T<:AbstractFloat}(img::AbstractArray{T}) = MapNone{U8}()
for ACV in (Color, AbstractRGB)
    for CV in subtypes(ACV)
        (length(CV.parameters) == 1 && !(CV.abstract)) || continue
        CVnew = CV<:AbstractGray ? Gray : RGB
        @eval mapinfo{T<:UFixed}(img::AbstractArray{$CV{T}}) = MapNone{$CVnew{T}}()
        @eval mapinfo{CV<:$CV}(img::AbstractArray{CV}) = MapNone{$CVnew{U8}}()
    end
end
mapinfo(img::AbstractArray{RGB24}) = MapNone{RGB{U8}}()

end # module


text:
using DiffOpt
using MathOptInterface
using Test
using OSQP
using Ipopt
using Clp
using Random
using SCS
using LinearAlgebra
using DelimitedFiles
using GLPK

const MOI = MathOptInterface;
const MOIU = MathOptInterface.Utilities;
const MOIT = MathOptInterface.Test

const ATOL = 1e-4
const RTOL = 1e-4

@testset "Examples" begin
    include(joinpath(@__DIR__, "../examples/solve-LP.jl"))
    include(joinpath(@__DIR__, "../examples/solve-QP.jl"))
end

@testset "Generate random problems" begin
    include("gen_random_problem.jl")
end

@testset "MOI_wrapper" begin
    include("MOI_wrapper.jl")
end

@testset "Utility Methods" begin
    include("utils.jl")
end

@testset "Solver Interface" begin
    include("solver_interface.jl")
end


text:
function setupCommBuffers!(domain::Domain, edgeNodes::Int)
    function cache_align_real(n::Int64)
        n = Int32(n)
        (n + CACHE_COHERENCE_PAD_REAL%Int32 - one(Int32)) & ~(CACHE_COHERENCE_PAD_REAL%Int32 - one(Int32))
    end
    @unpack_Domain domain
    # allocate a buffer large enough for nodal ghost data
    maxEdgeSize = max(domain.sizeX, max(domain.sizeY, domain.sizeZ))+1

    m_maxPlaneSize = cache_align_real(maxEdgeSize*maxEdgeSize)
    m_maxEdgeSize = cache_align_real(maxEdgeSize)
    maxPlaneSize = m_maxPlaneSize
    maxEdgeSize = m_maxEdgeSize

    # assume communication to 6 neighbors by default
    m_rowMin = (m_rowLoc == 0)        ? 0 : 1
    m_rowMax = (m_rowLoc == m_tp-1)     ? 0 : 1
    m_colMin = (m_colLoc == 0)        ? 0 : 1
    m_colMax = (m_colLoc == m_tp-1)     ? 0 : 1
    m_planeMin = (m_planeLoc == 0)    ? 0 : 1
    m_planeMax = (m_planeLoc == m_tp-1) ? 0 : 1

    # account for face communication
    comBufSize =(
        (m_rowMin + m_rowMax + m_colMin + m_colMax + m_planeMin + m_planeMax) *
        m_maxPlaneSize * MAX_FIELDS_PER_MPI_COMM
    )

    # account for edge communication
    comBufSize += (
        ((m_rowMin & m_colMin) + (m_rowMin & m_planeMin) + (m_colMin & m_planeMin) +
        (m_rowMax & m_colMax) + (m_rowMax & m_planeMax) + (m_colMax & m_planeMax) +
        (m_rowMax & m_colMin) + (m_rowMin & m_planeMax) + (m_colMin & m_planeMax) +
        (m_rowMin & m_colMax) + (m_rowMax & m_planeMin) + (m_colMax & m_planeMin)) *
        m_maxEdgeSize * MAX_FIELDS_PER_MPI_COMM
    )

    # account for corner communication
    # factor of 16 is so each buffer has its own cache line
    comBufSize += (((m_rowMin & m_colMin & m_planeMin) +
            (m_rowMin & m_colMin & m_planeMax) +
            (m_rowMin & m_colMax & m_planeMin) +
            (m_rowMin & m_colMax & m_planeMax) +
            (m_rowMax & m_colMin & m_planeMin) +
            (m_rowMax & m_colMin & m_planeMax) +
            (m_rowMax & m_colMax & m_planeMin) +
            (m_rowMax & m_colMax & m_planeMax)) * CACHE_COHERENCE_PAD_REAL
            )

    commDataSend = Vector{Float64}(undef, comBufSize)
    commDataRecv = Vector{Float64}(undef, comBufSize)
    # prevent floating point exceptions
    fill!(commDataSend, 0)
    fill!(commDataRecv, 0)

    # Boundary nodesets
    if (m_colLoc == 0)
        resize!(symmX, edgeNodes*edgeNodes)
    end
    if (m_rowLoc == 0)
        resize!(symmY, edgeNodes*edgeNodes)
    end
    if (m_planeLoc == 0)
        resize!(symmZ, edgeNodes*edgeNodes)
    end
    @pack_Domain! domain
end

getMyRank(comm::MPI.Comm) = MPI.Comm_rank(comm)
getMyRank(::Nothing) = 0

getNumRanks(comm::MPI.Comm) = MPI.Comm_size(comm)
getNumRanks(::Nothing) = 1

getWtime(::MPI.Comm) = MPI.Wtime()
getWtime(::Nothing) = time()

comm_max(data::Float64, comm::MPI.Comm) = MPI.Allreduce(data, MPI.MAX, comm)
comm_max(data::Float64, ::Nothing) = data

comm_min(data::Float64, comm::MPI.Comm) = MPI.Allreduce(data, MPI.MIN, comm)
comm_min(data::Float64, ::Nothing) = data

comm_barrier(comm::MPI.Comm) = MPI.Barrier(comm)
comm_barrier(::Nothing) = nothing


text:
using ExploratoryDataAnalysis
using DataFrames
using CategoricalArrays
# using Plots

df = DataFrame(t = [fill(0, 500); fill(1, 500)], x1 = collect(1:1000), x2 = randn(1000))

df.tc = categorical(df.t)

eda(df, :t)

eda(df, :x1)

f = [9800 9925 11825; 190200 190075 188175]'
infovalue(f)
mutualinfo(f)


text:
# TODO - Check deg2rad
function calculateTheta(configData)
    theta = deg2rad(Float64(configData["core"]["aoa"]))::Float64
    return theta
end

function compute_cl_cd_cm(globaldata, configData, shapeindices)

    rho_inf::Float64 = configData["core"]["rho_inf"]::Float64
    Mach::Float64 = configData["core"]["mach"]::Float64
    pr_inf::Float64 = configData["core"]["pr_inf"]::Float64
    shapes::Int64 = configData["core"]["shapes"]::Float64
    theta = calculateTheta(configData)

    temp = 0.5*rho_inf*Mach*Mach

    H = zeros(Float64, shapes)
    V = zeros(Float64, shapes)
    pitch_mom = zeros(Float64, shapes)

    Cl = zeros(Float64, shapes)
    Cd = zeros(Float64, shapes)
    Cm = zeros(Float64, shapes)

    open("cp_file", "w") do the_file
        for itm in shapeindices
            # print(wallindices)
            left = globaldata[itm].left
            right = globaldata[itm].right
            lx = globaldata[left].x
            ly = globaldata[left].y
            rx = globaldata[right].x
            ry = globaldata[right].y
            mx = globaldata[itm].x
            my = globaldata[itm].y

            ds1 = hypot(mx - lx, my - ly)
            ds2 = hypot(rx - mx, ry - my)

            ds = 0.5*(ds1 + ds2)
            nx = globaldata[itm].nx
            ny = globaldata[itm].ny
            cp = globaldata[itm].prim[4] - pr_inf
            cp = -cp/temp

            flag_2 = globaldata[itm].flag_2

            print(the_file, flag_2, " ", mx, " ",  cp, "\n")
            H[flag_2] += cp * nx * ds
            V[flag_2] += cp * ny * ds

            pitch_mom[flag_2] += (-cp * ny * ds * (mx - 0.25)) + (cp * nx * ds * my)
        end
    end

    Cl = V*cos(theta) - H*sin(theta)
    Cd = H*cos(theta) + V*sin(theta)
    Cm = pitch_mom

end

max_stars_repo_name: TestSubjector/HybridMeshfreeSolver

max_issues_repo_name: TestSubjector/HybridMeshfreeSolver

max_forks_repo_name: TestSubjector/HybridMeshfreeSolver


text:
function test11()
  filename = "test11.asc"
  # exectuablepath = null string will not run LTspice.exe.  Test parsing only.
  sim = LTspiceSimulation(filename,executablepath="")
  @test LTspice.does_circuitfilearray_file_match(sim)
  show(IOBuffer(),sim)

  asteps = [1.0,2.0]
  bsteps = [10.0,15.0,20.0,25.0]
  csteps = [100.0,200.0,300.0]

  @test stepvalues(sim) == (asteps,bsteps,csteps)
  @test measurementnames(sim) == ("sum","sump1000")
  @test stepnames(sim) == ("a","b","c")
  @test logpath(sim) != ""

  @static if Sys.iswindows()
      @test circuitpath(sim) == filename
  end

  @test parametervalues(sim) == []
  @test measurementvalues(sim)[1,1,1,1] == 111.0
  @test executablepath(sim) == ""
  @test haskey(sim,"sum") == true  # measurments in stepped files are not a Dict
  @test length(sim) == 48

  for i in eachindex(asteps)
      for j in eachindex(bsteps)
          for k in eachindex(csteps)
              @test measurementvalues(sim)[i,j,k,1] == asteps[i]+bsteps[j]+csteps[k]
              @test sim["sum"][i,j,k] == asteps[i]+bsteps[j]+csteps[k]
              @test measurementvalues(sim)[i,j,k,2] == asteps[i]+bsteps[j]+csteps[k]+1000
              @test sim["sump1000"][i,j,k] == asteps[i]+bsteps[j]+csteps[k]+1000
          end
      end
  end

  for line in perlineiterator(sim)
      @test (line[1]+line[2]+line[3] == line[4])
      @test (line[4] + 1000 == line[5])
  end

  for line in perlineiterator(sim,resultnames=["sum"],steporder=["b","a","c"])
      @test (line[1]+line[2]+line[3] == line[4])
  end
  show(IOBuffer(),sim)
end
test11()


text:
module ShiftedArrays

import Base: checkbounds, getindex, setindex!, parent, size
export ShiftedArray, ShiftedVector, shifts, default
export CircShiftedArray, CircShiftedVector
export lag, lead
export reduce_vec, mapreduce_vec

include("shiftedarray.jl")
include("circshiftedarray.jl")
include("lag.jl")
include("circshift.jl")

end


text:
"""
    parsimonious_flux_balance_analysis(
        model::MetabolicModel,
        optimizer;
        modifications = [],
        qp_modifications = [],
        relax_bounds=[1.0, 0.999999, 0.99999, 0.9999, 0.999, 0.99],
    )

Run parsimonious flux balance analysis (pFBA) on the `model`. In short, pFBA
runs two consecutive optimization problems. The first is traditional FBA:
```
max cᵀx = μ
s.t. S x = b
     xₗ ≤ x ≤ xᵤ
```
And the second is a quadratic optimization problem:
```
min Σᵢ xᵢ²
s.t. S x = b
     xₗ ≤ x ≤ xᵤ
     μ = μ⁰
```
Where the optimal solution of the FBA problem, μ⁰, has been added as an
additional constraint. See "Lewis, Nathan E, Hixson, Kim K, Conrad, Tom M,
Lerman, Joshua A, Charusanti, Pep, Polpitiya, Ashoka D, Adkins, Joshua N,
Schramm, Gunnar, Purvine, Samuel O, Lopez‐Ferrer, Daniel, Weitz, Karl K, Eils,
Roland, König, Rainer, Smith, Richard D, Palsson, Bernhard Ø, (2010) Omic data
from evolved E. coli are consistent with computed optimal growth from
genome‐scale models. Molecular Systems Biology, 6. 390. doi:
accession:10.1038/msb.2010.47" for more details.

pFBA gets the model optimum by standard FBA (using
[`flux_balance_analysis`](@ref) with `optimizer` and `modifications`), then
finds a minimal total flux through the model that still satisfies the (slightly
relaxed) optimum. This is done using a quadratic problem optimizer. If the
original optimizer does not support quadratic optimization, it can be changed
using the callback in `qp_modifications`, which are applied after the FBA.

Thhe optimum relaxation sequence can be specified in `relax` parameter, it
defaults to multiplicative range of `[1.0, 0.999999, ..., 0.99]` of the
original bound.

Returns an optimized model that contains the pFBA solution; or `nothing` if the
optimization failed.

# Example
```
optimizer = Gurobi.Optimizer
atts = Dict("OutputFlag" => 0)
model = load_model(StandardModel, "iJO1366.json")
biomass = findfirst(model.reactions, "BIOMASS_Ec_iJO1366_WT_53p95M")
sol = pfba(model, biomass, Gurobi.optimizer)
```
"""
function parsimonious_flux_balance_analysis(
    model::MetabolicModel,
    optimizer;
    modifications = [],
    qp_modifications = [],
    relax_bounds = [1.0, 0.999999, 0.99999, 0.9999, 0.999, 0.99],
)
    # Run FBA
    opt_model = flux_balance_analysis(model, optimizer; modifications = modifications)
    is_solved(opt_model) || return nothing # FBA failed

    # get the objective
    Z = objective_value(opt_model)
    original_objective = COBREXA.JuMP.objective_function(opt_model)

    # prepare the model for pFBA
    for mod in qp_modifications
        mod(model, opt_model)
    end

    # add the minimization constraint for total flux
    v = opt_model[:x] # fluxes
    @objective(opt_model, Min, sum(dot(v, v)))

    for rb in relax_bounds
        lb, ub = objective_bounds(rb)(Z)
        @_models_log @info "pFBA step relaxed to [$lb,$ub]"
        @constraint(opt_model, pfba_constraint, lb <= original_objective <= ub)

        optimize!(opt_model)
        is_solved(opt_model) && break

        COBREXA.JuMP.delete(opt_model, pfba_constraint)
        COBREXA.JuMP.unregister(opt_model, :pfba_constraint)
    end

    is_solved(opt_model) || return nothing # pFBA failed

    return opt_model
end

"""
    parsimonious_flux_balance_analysis_vec(args...; kwargs...)

Perform parsimonious flux balance analysis on `model` using `optimizer`.
Returns a vector of fluxes in the same order as the reactions in `model`.
Arguments are forwarded to [`parsimonious_flux_balance_analysis`](@ref) internally.
"""
function parsimonious_flux_balance_analysis_vec(args...; kwargs...)
    opt_model = parsimonious_flux_balance_analysis(args...; kwargs...)

    isnothing(opt_model) && return nothing

    return value.(opt_model[:x])
end

"""
    parsimonious_flux_balance_analysis_dict(model::MetabolicModel, args...; kwargs...)

Perform parsimonious flux balance analysis on `model` using `optimizer`.
Returns a dictionary mapping the reaction IDs to fluxes.
Arguments are forwarded to [`parsimonious_flux_balance_analysis`](@ref) internally.
"""
function parsimonious_flux_balance_analysis_dict(model::MetabolicModel, args...; kwargs...)
    opt_fluxes = parsimonious_flux_balance_analysis_vec(model, args...; kwargs...)

    isnothing(opt_fluxes) && return nothing

    return Dict(zip(reactions(model), opt_fluxes))
end


text:
##### Momentum tendencies

export GeostrophicForcing
export Coriolis
export Gravity
export RayleighSponge
using CLIMAParameters.Planet: Omega

#####
##### First order fluxes
#####

function flux(::Momentum, ::Advect, atmos, args)
    @unpack state = args
    return state.ρu .* (state.ρu / state.ρ)'
end

function flux(::Momentum, ::PressureGradient, atmos, args)
    @unpack state, aux = args
    @unpack ts = args.precomputed
    s = state.ρu * state.ρu'
    pad = SArray{Tuple{size(s)...}}(ntuple(i -> 0, length(s)))
    ref_state = reference_state(atmos)
    if ref_state isa HydrostaticState && ref_state.subtract_off
        return pad + (air_pressure(ts) - aux.ref_state.p) * I
    else
        return pad + air_pressure(ts) * I
    end
end

#####
##### Second order fluxes
#####

struct ViscousStress <: TendencyDef{Flux{SecondOrder}} end

function flux(::Momentum, ::ViscousStress, atmos, args)
    @unpack state = args
    @unpack τ = args.precomputed.turbulence
    pad = SArray{Tuple{size(τ)...}}(ntuple(i -> 0, length(τ)))
    return pad + τ * state.ρ
end

function flux(::Momentum, ::MoistureDiffusion, atmos, args)
    @unpack state, diffusive = args
    @unpack D_t = args.precomputed.turbulence
    d_q_tot = (-D_t) .* diffusive.moisture.∇q_tot
    return d_q_tot .* state.ρu'
end

function flux(::Momentum, ::HyperdiffViscousFlux, atmos, args)
    @unpack state, hyperdiffusive = args
    return state.ρ * hyperdiffusive.hyperdiffusion.ν∇³u_h
end

#####
##### Sources
#####

struct Gravity <: TendencyDef{Source} end

prognostic_vars(::Gravity) = (Momentum(),)

function source(::Momentum, ::Gravity, m, args)
    @unpack state, aux = args
    ref_state = reference_state(m)
    if ref_state isa HydrostaticState && ref_state.subtract_off
        return -(state.ρ - aux.ref_state.ρ) * aux.orientation.∇Φ
    else
        return -state.ρ * aux.orientation.∇Φ
    end
end

struct Coriolis <: TendencyDef{Source} end

prognostic_vars(::Coriolis) = (Momentum(),)

function source(::Momentum, ::Coriolis, m, args)
    @unpack state = args
    FT = eltype(state)
    param_set = parameter_set(m)
    _Omega::FT = Omega(param_set)
    # note: this assumes a SphericalOrientation
    return -SVector(0, 0, 2 * _Omega) × state.ρu
end

struct GeostrophicForcing{FT} <: TendencyDef{Source}
    f_coriolis::FT
    u_geostrophic::FT
    v_geostrophic::FT
end

prognostic_vars(::GeostrophicForcing) = (Momentum(),)

function source(::Momentum, s::GeostrophicForcing, m, args)
    @unpack state, aux = args
    u_geo = SVector(s.u_geostrophic, s.v_geostrophic, 0)
    ẑ = vertical_unit_vector(m, aux)
    fkvector = s.f_coriolis * ẑ
    return -fkvector × (state.ρu .- state.ρ * u_geo)
end

"""
    RayleighSponge{FT} <: TendencyDef{Source}

Rayleigh Damping (Linear Relaxation) for top wall momentum components
Assumes laterally periodic boundary conditions for LES flows. Momentum components
are relaxed to reference values (zero velocities) at the top boundary.
"""
struct RayleighSponge{FT} <: TendencyDef{Source}
    "Maximum domain altitude (m)"
    z_max::FT
    "Altitude at with sponge starts (m)"
    z_sponge::FT
    "Sponge Strength 0 ⩽ α_max ⩽ 1"
    α_max::FT
    "Relaxation velocity components"
    u_relaxation::SVector{3, FT}
    "Sponge exponent"
    γ::FT
end

prognostic_vars(::RayleighSponge) = (Momentum(),)

function source(::Momentum, s::RayleighSponge, m, args)
    @unpack state, aux = args
    z = altitude(m, aux)
    if z >= s.z_sponge
        r = (z - s.z_sponge) / (s.z_max - s.z_sponge)
        β_sponge = s.α_max * sinpi(r / 2)^s.γ
        return -β_sponge * (state.ρu .- state.ρ * s.u_relaxation)
    else
        FT = eltype(state)
        return SVector{3, FT}(0, 0, 0)
    end
end


text:
export whiten_input

function whiten_input(X::Array{Float32,4})
    for j = 1:size(X, 4)
        x = X[:, :, :, j:j]
        X[:, :, :, j:j] = (x .- mean(x))/std(x)
    end
    return X
end


text:
#include("Header.jl")

function SeisHeaderInfo(filename;ntrace=100000)
	#
	#   print Seis header information
	#

	@compat	key = fieldnames(Header)
	nhead = length(key)
	filename_headers = ParseHeaderName(filename)
	stream = open(filename_headers)
	NX = GetNumTraces(filename)
	h = GrabHeader(stream,1)
	println("Displaying information for ", filename," (",NX," traces):")
	min_h = zeros(Float32,length(key))
	max_h = zeros(Float32,length(key))
	mean_h = zeros(Float32,length(key))

	for ikey=1:length(key)
		min_h[ikey] = convert(Float32,getfield(h,key[ikey]))
		max_h[ikey] = convert(Float32,getfield(h,key[ikey]))
		mean_h[ikey] += convert(Float32,getfield(h,key[ikey]))
	end

	itrace = 2
	while itrace <= NX
		nx = NX - itrace + 1
		ntrace = nx > ntrace ? ntrace : nx
		position = 4*nhead*(itrace-1)
		seek(stream,position)
		h1 = read(stream,Header32Bits,nhead*ntrace)
		h1 = reshape(h1,nhead,convert(Int,ntrace))
		for ikey = 1 : length(key)
			keytype = eval(parse("typeof(Seismic.InitSeisHeader().$(string(key[ikey])))"))
			h2 = reinterpret(keytype,vec(h1[ikey,:]))
			a = minimum(h2)
			b = maximum(h2)
			c = mean(h2)
			if (a < min_h[ikey])
				min_h[ikey] = a
			end
			if (b > max_h[ikey])
				max_h[ikey] = b
			end
			mean_h[ikey] += c*ntrace
		end
		itrace += ntrace
	end

	for ikey=1:length(key)
		mean_h[ikey] /= NX
	end
	close(stream)
	println("       Key          Minimum          Maximum             Mean");
	println("=============================================================")
	for ikey=1:length(key)
		@printf("%10s      %11.3f      %11.3f      %11.3f\n",string(key[ikey]),min_h[ikey],max_h[ikey],mean_h[ikey])
	end

end


text:
# Stream management

export
    CuStream, CuDefaultStream, synchronize

mutable struct CuStream
    handle::CUstream
    ctx::CuContext
end

Base.unsafe_convert(::Type{CUstream}, s::CuStream) = s.handle

Base.:(==)(a::CuStream, b::CuStream) = a.handle == b.handle
Base.hash(s::CuStream, h::UInt) = hash(s.handle, h)

@enum_without_prefix CUstream_flags_enum CU_

"""
    CuStream(flags=STREAM_DEFAULT)

Create a CUDA stream.
"""
function CuStream(flags::CUstream_flags=STREAM_DEFAULT)
    handle_ref = Ref{CUstream}()
    cuStreamCreate(handle_ref, flags)

    ctx = CuCurrentContext()
    obj = CuStream(handle_ref[], ctx)
    finalizer(unsafe_destroy!, obj)
    return obj
end

function unsafe_destroy!(s::CuStream)
    if isvalid(s.ctx)
        cuStreamDestroy(s)
    end
end

"""
    CuDefaultStream()

Return the default stream.
"""
@inline CuDefaultStream() = CuStream(convert(CUstream, C_NULL), CuContext(C_NULL))

"""
    synchronize(s::CuStream)

Wait until a stream's tasks are completed.
"""
synchronize(s::CuStream) = cuStreamSynchronize(s)


text:
module LambdaMain

using HTTP
#=
using AWSS3

function lambda_main(event, context)
    println(s3_list_buckets())
end
=#

Base.@ccallable function julia_main(ARGS::Vector{String})::Cint
    #lambda_main(ARGS[1], ARGS[2])
    println(HTTP.get("http://httpbin.org/ip"))
    return 0
end

end # module LambdaMain


text:
module Mazes
using SimpleGraphs, SimplePartitions, SimpleDrawing, Plots

_vtx = Tuple{Int,Int}
_edge = Tuple{_vtx,_vtx}

import SimpleDrawing: draw
import Base: show, size
export Maze, draw_ans, draw, amaze

"""
`M = Maze(r,c)` creates an `r`-by-`c` maze. To see the maze,
use `draw(M)`. The maze is meant to be solved from the upper left
to the lower right. This is indicated by stars at those locations.
To draw the maze without those stars, use `draw(M,false)`.

To draw the solution to the maze, use `draw_ans(M)`. The solutions is
always drawn without stars.
"""
struct Maze
    r::Int  # number of rows
    c::Int  # number of cols
    T::SimpleGraph{_vtx}  # the tree for this maze
    function Maze(nr::Int, nc::Int)
        G = _gen_tree(nr, nc)
        new(nr, nc, G)
    end
end

Maze(n::Int) = Maze(n, n)

function show(io::IO, M::Maze)
    print(io, "Maze($(M.r),$(M.c))")
end

"""
`size(M::Maze)` returns a tuple `(r,c)` where `r` is the number of
rows and `c` is the number of columns in the maze.
"""
size(M::Maze) = (M.r, M.c)

function _gen_tree(nr, nc)
    @assert nr >= 2 && nc >= 2 "Inputs must both be at least 2"
    G = Grid(nr, nc)

    EE = elist(G)  # all possible edges
    # assign random weights to edges
    wt = Dict{_edge,Float64}()
    for e in EE
        wt[e] = rand()
    end

    # the upper left and lower right squares should be leaves
    # these weights should enforce that

    e1 = ((1, 1), (1, 2))
    e2 = ((1, 1), (2, 1))
    if rand() > 0.5
        wt[e1] = 0
        wt[e2] = 1
    else
        wt[e1] = 1
        wt[e2] = 0
    end

    e1 = (nr - 1, nc), (nr, nc)
    e2 = (nr, nc - 1), (nr, nc)
    if rand() > 0.5
        wt[e1] = 0
        wt[e2] = 1
    else
        wt[e1] = 1
        wt[e2] = 0
    end

    wt_list = [wt[e] for e in EE]
    p = sortperm(wt_list)
    E_list = [EE[j] for j in p] # EE sorted by weight

    # copy vertices from G
    T = SimpleGraph{_vtx}()
    for v in G.V
        add!(T, v)
    end
    P = Partition(G.V)

    for e in E_list
        a = e[1]
        b = e[2]
        if !in_same_part(P, a, b)
            add!(T, a, b)
            merge_parts!(P, a, b)
        end
    end

    d = Dict{_vtx,Array{Int,1}}()
    for v in T.V
        x = v[1]
        y = v[2]
        d[v] = [y, -x]
    end
    embed(T, d)
    return T
end

function _puzzle_draw(M::Maze)
    # get the tree's embedding
    xy = M.T.cache[:xy]

    newdraw()
    # draw(M.T)  # DEBUG #

    # draw outside rectangle
    ul = (0.5, -0.5)
    ur = (M.c + 0.5, -0.5)
    ll = (0.5, -M.r - 0.5)
    lr = (M.c + 0.5, -M.r - 0.5)

    draw_segment(ul..., ur..., color = :black)
    draw_segment(ul..., ll..., color = :black)
    draw_segment(ur..., lr..., color = :black)
    draw_segment(ll..., lr..., color = :black)

    G = Grid(M.r, M.c)

    non_edge_list = [e for e in G.E if !has(M.T, e[1], e[2])]

    for e in non_edge_list
        a = e[1]
        b = e[2]
        mid = 0.5 * (xy[a] + xy[b])

        if a[1] == b[1]  # vertical segment
            x = mid[1]
            y1 = mid[2] - 0.5
            y2 = mid[2] + 0.5
            draw_segment(x, y1, x, y2, color = :black)
        else
            x1 = mid[1] - 0.5
            x2 = mid[1] + 0.5
            y = mid[2]
            draw_segment(x1, y, x2, y, color = :black)
        end

    end
end

"""
`draw(M::Maze)` draws the maze `M` on the screen. Circles mark the start
(upper left) and end (lower right). The size is printed above the maze.

Full syntax: `draw(M::Maze, markers::Bool=true, title::Bool=true)`
"""
function draw(M::Maze, markers::Bool = true, title::Bool = true)
    _puzzle_draw(M)

    if markers
        r = 0.125
        draw_circle(1, -1, r, color = :black)
        draw_circle(M.c, -M.r, r, color = :black)
    end

    if title
        plot!(title = "$(M.r)-by-$(M.c)")
    end

    finish()
end

function _ans_draw(M::Maze, s::_vtx, t::_vtx)
    xy = M.T.cache[:xy]
    @assert has(M.T, s) "$s not a square in this maze"
    @assert has(M.T, t) "$t not a square in this maze"
    if s == t
        @warn "Two locations are the same, no path will be drawn"
    end
    P = find_path(M.T, s, t)
    n = length(P)
    for i = 1:n-1
        a = P[i]
        xa, ya = xy[a]
        b = P[i+1]
        xb, yb = xy[b]
        draw_segment(xa, ya, xb, yb, color = :red, style = :dot)
    end
end

"""
`draw_ans(M::Maze)` draws the maze `M` and the solution (as a dotted red line)
from the start (upper left) to the end (lower right).

To draw a different path (solution) use `draw(M,s,t)` where `s` and `t` are
integer 2-tuples. Note that squares in the maze are indexed like a matrix,
so `(2,3)` refers to the cell in row 2, column 3.
"""
function draw_ans(M::Maze, s::_vtx, t::_vtx)
    _puzzle_draw(M)
    _ans_draw(M, s, t)
    finish()
end

draw_ans(M::Maze) = draw_ans(M, (1, 1), size(M))

"""
`amaze(r,c)` creates and draws a maze.
"""
function amaze(r::Int, c::Int)
    M = Maze(r, c)
    draw(M)
end

end # end of module


text:
# This tutorial uses the TMAR Filter from
#
#    @article{doi:10.1175/MWR-D-16-0220.1,
#      author = {Light, Devin and Durran, Dale},
#      title = {Preserving Nonnegativity in Discontinuous Galerkin
#               Approximations to Scalar Transport via Truncation and Mass
#               Aware Rescaling (TMAR)},
#      journal = {Monthly Weather Review},
#      volume = {144},
#      number = {12},
#      pages = {4771-4786},
#      year = {2016},
#      doi = {10.1175/MWR-D-16-0220.1},
#    }
#
# to reproduce the tutorial in section 4b.  It is a shear swirling
# flow deformation of a transported quantity from LeVeque 1996.  The exact
# solution at the final time is the same as the initial condition.

using MPI
using Test
using ClimateMachine
ClimateMachine.init()
using Logging
using ClimateMachine.Mesh.Topologies
using ClimateMachine.Mesh.Grids
using ClimateMachine.Mesh.Filters
using ClimateMachine.DGMethods
using ClimateMachine.DGMethods.NumericalFluxes
using ClimateMachine.MPIStateArrays
using ClimateMachine.ODESolvers
using LinearAlgebra
using Printf
using Dates
using ClimateMachine.GenericCallbacks:
    EveryXWallTimeSeconds, EveryXSimulationSteps
using ClimateMachine.VTK: writevtk, writepvtu

using ClimateMachine
const clima_dir = dirname(dirname(pathof(ClimateMachine)));
include(joinpath(
    clima_dir,
    "test",
    "Numerics",
    "DGMethods",
    "advection_diffusion",
    "advection_diffusion_model.jl",
))

Base.@kwdef struct SwirlingFlow{FT} <: AdvectionDiffusionProblem
    period::FT = 5
end

init_velocity_diffusion!(::SwirlingFlow, aux::Vars, geom::LocalGeometry) =
    nothing

cosbell(τ, q) = τ ≤ 1 ? ((1 + cospi(τ)) / 2)^q : zero(τ)

function initial_condition!(::SwirlingFlow, state, aux, coord, t)
    FT = eltype(state)
    x, y, _ = aux.coord
    x0, y0 = FT(1 // 4), FT(1 // 4)
    τ = 4 * hypot(x - x0, y - y0)
    state.ρ = cosbell(τ, 3)
end;

has_variable_coefficients(::SwirlingFlow) = true
function update_velocity_diffusion!(
    problem::SwirlingFlow,
    ::AdvectionDiffusion,
    state::Vars,
    aux::Vars,
    t::Real,
)
    x, y, _ = aux.coord
    sx, cx = sinpi(x), cospi(x)
    sy, cy = sinpi(y), cospi(y)
    ct = cospi(t / problem.period)

    u = 2 * sx^2 * sy * cy * ct
    v = -2 * sy^2 * sx * cx * ct
    aux.u = SVector(u, v, 0)
end;

function do_output(mpicomm, vtkdir, vtkstep, dg, Q, model, testname)
    ## name of the file that this MPI rank will write
    filename = @sprintf(
        "%s/%s_mpirank%04d_step%04d",
        vtkdir,
        testname,
        MPI.Comm_rank(mpicomm),
        vtkstep
    )

    statenames = flattenednames(vars_state(model, Prognostic(), eltype(Q)))

    writevtk(filename, Q, dg, statenames)

    ## generate the pvtu file for these vtk files
    if MPI.Comm_rank(mpicomm) == 0
        ## name of the pvtu file
        pvtuprefix = @sprintf("%s/%s_step%04d", vtkdir, testname, vtkstep)

        ## name of each of the ranks vtk files
        prefixes = ntuple(MPI.Comm_size(mpicomm)) do i
            @sprintf("%s_mpirank%04d_step%04d", testname, i - 1, vtkstep)
        end

        writepvtu(pvtuprefix, prefixes, statenames, eltype(Q))

        @info "Done writing VTK: $pvtuprefix"
    end
end;

function test_run(
    mpicomm,
    ArrayType,
    topl,
    problem,
    dt,
    N,
    timeend,
    FT,
    vtkdir,
    outputtime,
)
    grid = DiscontinuousSpectralElementGrid(
        topl,
        FloatType = FT,
        DeviceArray = ArrayType,
        polynomialorder = N,
    )

    model = AdvectionDiffusion{2, false, true}(problem)

    dg = DGModel(
        model,
        grid,
        UpwindNumericalFlux(),
        CentralNumericalFluxSecondOrder(),
        CentralNumericalFluxGradient(),
    )

    Q = init_ode_state(dg, FT(0))

    initialsumQ = weightedsum(Q)

    ## We integrate so that the final solution is equal to the initial solution
    Qe = copy(Q)

    rhs! = function (dQdt, Q, ::Nothing, t; increment = false)
        Filters.apply!(Q, :, grid, TMARFilter())
        dg(dQdt, Q, nothing, t; increment = false)
    end

    odesolver = SSPRK33ShuOsher(rhs!, Q; dt = dt, t0 = 0)

    cbTMAR = EveryXSimulationSteps(1) do
        Filters.apply!(Q, :, grid, TMARFilter())
    end

    ## create output directory on first rank of communicator
    if MPI.Comm_rank(mpicomm) == 0
        mkpath(vtkdir)
    end
    MPI.Barrier(mpicomm)

    vtkstep = 0
    ## output initial step
    do_output(mpicomm, vtkdir, vtkstep, dg, Q, model, "nonnegative")
    ## setup the output callback
    cbvtk = EveryXSimulationSteps(floor(outputtime / dt)) do
        vtkstep += 1
        minQ, maxQ = minimum(Q), maximum(Q)
        sumQ = weightedsum(Q)

        sumerror = (initialsumQ - sumQ) / initialsumQ

        @info @sprintf """Step  = %d
          minimum(Q)  = %.16e
          maximum(Q)  = %.16e
          sum error   = %.16e
          """ vtkstep minQ maxQ sumerror

        do_output(mpicomm, vtkdir, vtkstep, dg, Q, model, "nonnegative")
    end

    callbacks = (cbTMAR, cbvtk)
    solve!(Q, odesolver; timeend = timeend, callbacks = callbacks)

    minQ, maxQ = minimum(Q), maximum(Q)
    finalsumQ = weightedsum(Q)
    sumerror = (initialsumQ - finalsumQ) / initialsumQ
    error = euclidean_distance(Q, Qe)

    @test minQ ≥ 0

    @info @sprintf """Finished
    minimum(Q) = %.16e
    maximum(Q) = %.16e
    L2 error   = %.16e
    sum error  = %.16e
    """ minQ maxQ error sumerror
end;

let
    ArrayType = ClimateMachine.array_type()

    mpicomm = MPI.COMM_WORLD

    FT = Float64
    dim = 2
    Ne = 20
    polynomialorder = 4

    problem = SwirlingFlow()

    brickrange = (
        range(FT(0); length = Ne + 1, stop = 1),
        range(FT(0); length = Ne + 1, stop = 1),
        range(FT(0); length = Ne + 1, stop = 1),
    )

    topology = BrickTopology(
        mpicomm,
        brickrange[1:dim],
        boundary = ntuple(d -> (3, 3), dim),
    )

    maxvelocity = 2
    elementsize = 1 / Ne
    dx = elementsize / polynomialorder^2
    CFL = 1
    dt = CFL * dx / maxvelocity

    vtkdir =
        abspath(joinpath(ClimateMachine.Settings.output_dir, "vtk_nonnegative"))
    outputtime = 0.0625
    dt = outputtime / ceil(Int64, outputtime / dt)

    timeend = problem.period

    @info @sprintf """Starting
    FT               = %s
    dim              = %d
    Ne               = %d
    polynomial order = %d
    final time       = %.16e
    time step        = %.16e
    """ FT dim Ne polynomialorder timeend dt

    test_run(
        mpicomm,
        ArrayType,
        topology,
        problem,
        dt,
        polynomialorder,
        timeend,
        FT,
        vtkdir,
        outputtime,
    )
end;

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
module GRF

using ARCHModels
using Cbc
using CSV
using Dates
using DataFrames
using Distributions
using GLM
using HypothesisTests
using HTTP
using Ipopt
using JSON
using JuMP
using Juniper
using LinearAlgebra
using LsqFit
using Random
using Statistics
using StatsPlots
using XLSX

include("modulo01.jl")

# Notebooks interativos
# binder: https://mybinder.org
# https://mybinder.org/v2/gh/ASaragga/GRF.jl/HEAD
#
# nbviewer: https://nbviewer.jupyter.org
# https://github.com/ASaragga/GRF.jl/blob/master/src/Modulo1.1-1.7.ipynb

end # module