text:
using BinDeps
using Compat
import Compat.String

@BinDeps.setup

libnames    = ["libMagickWand", "CORE_RL_wand_"]
suffixes    = ["", "-Q16", "-6.Q16", "-Q8"]
options     = ["", "HDRI"]
extensions  = ["", ".so.2", ".so.4", ".so.5"]
aliases     = vec(libnames .*
                  reshape(suffixes,(1,length(suffixes))) .*
                  reshape(options,(1,1,length(options))) .*
                  reshape(extensions,(1,1,1,length(extensions))))
libwand     = library_dependency("libwand", aliases = aliases)

initfun = """
function init_deps()
    ccall((:MagickWandGenesis,libwand), Void, ())
end
"""

mpath = get(ENV, "MAGICK_HOME", "") # If MAGICK_HOME is defined, add to library search path
if !isempty(mpath)
    provides(Binaries, mpath, libwand)
    provides(Binaries, joinpath(mpath,"lib"), libwand)
end

if is_linux()
    kwargs = Any[(:onload, initfun)]
    provides(AptGet, "libmagickwand4", libwand; kwargs...)
    provides(AptGet, "libmagickwand5", libwand; kwargs...)
    provides(AptGet, "libmagickwand-6.q16-2", libwand; kwargs...)
    provides(Pacman, "imagemagick", libwand; kwargs...)
    provides(Yum, "ImageMagick", libwand; kwargs...)
end

# TODO: remove me when upstream is fixed
is_windows() && push!(BinDeps.defaults, BuildProcess)

if is_windows()
    const OS_ARCH = (WORD_SIZE == 64) ? "x64" : "x86"

    # TODO: checksums: we have gpg
    # Extract the appropriate filename to download
    magick_base = "http://www.imagemagick.org/download/binaries"
    binariesfn  = download(magick_base)
    str         = readstring(binariesfn)
    pattern     = "ImageMagick-6.9.*?-Q16-$(OS_ARCH)-dll.exe"
    m           = match(Regex(pattern), str)
    magick_exe  = convert(String, m.match)

    magick_tmpdir   = BinDeps.downloadsdir(libwand)
    magick_url      = "$(magick_base)/$(magick_exe)"
    magick_libdir   = joinpath(BinDeps.libdir(libwand), OS_ARCH)
    innounp_url     = "https://bintray.com/artifact/download/julialang/generic/innounp.exe"
    initfun         =
"""
function init_deps()
    ENV["MAGICK_CONFIGURE_PATH"]    = \"$(escape_string(magick_libdir))\"
    ENV["MAGICK_CODER_MODULE_PATH"] = \"$(escape_string(magick_libdir))\"
end
init_deps()
"""

    provides(BuildProcess,
        (@build_steps begin
            CreateDirectory(magick_tmpdir)
            CreateDirectory(magick_libdir)
            FileDownloader(magick_url, joinpath(magick_tmpdir, magick_exe))
            FileDownloader(innounp_url, joinpath(magick_tmpdir, "innounp.exe"))
            @build_steps begin
                ChangeDirectory(magick_tmpdir)
                info("Installing ImageMagick library")
                `innounp.exe -q -y -b -e -x -d$(magick_libdir) $(magick_exe)`
            end
        end),
        libwand, os = :Windows, unpacked_dir = magick_libdir, preload = initfun)
end

if is_apple()
    if Pkg.installed("Homebrew") === nothing
        error("Homebrew package not installed, please run Pkg.add(\"Homebrew\")")
    end
    using Homebrew
    initfun =
"""
function init_deps()
    ENV["MAGICK_CONFIGURE_PATH"] = joinpath("$(Homebrew.prefix("imagemagick"))","lib","ImageMagick","config-Q16")
    ENV["MAGICK_CODER_MODULE_PATH"] = joinpath("$(Homebrew.prefix("imagemagick"))", "lib","ImageMagick","modules-Q16","coders")
    ENV["PATH"] = joinpath("$(Homebrew.prefix("imagemagick"))", "bin") * ":" * ENV["PATH"]
    ccall((:MagickWandGenesis,libwand), Void, ())
end
"""
    provides( Homebrew.HB, "imagemagick", libwand, os = :Darwin, preload = initfun, onload="init_deps()")
end

@BinDeps.install Dict([(:libwand, :libwand)])

# Hack-fix for issue #12
# Check to see whether init_deps is present, and if not add it
if isempty(search(readstring(joinpath(dirname(@__FILE__),"deps.jl")), "init_deps"))
    open("deps.jl", "a") do io
        write(io, initfun)
    end
end

# Save the library version; by checking this now, we avoid a runtime dependency on libwand
# See https://github.com/timholy/Images.jl/issues/184#issuecomment-55643225
module CheckVersion
using Compat
include("deps.jl")
p = ccall((:MagickQueryConfigureOption, libwand), Ptr{UInt8}, (Ptr{UInt8},), "LIB_VERSION_NUMBER")
vstr = string("v\"", join(split(unsafe_string(p), ',')[1:3], '.'), "\"")
open(joinpath(dirname(@__FILE__),"versioninfo.jl"), "w") do file
    write(file, "const libversion = $vstr\n")
end
end

is_windows() && pop!(BinDeps.defaults)


text:
# Autogenerated wrapper script for CUDA_full_jll for powerpc64le-linux-gnu
JLLWrappers.@generate_wrapper_header("CUDA_full")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@generate_init_footer()
end  # __init__()

ext: jl


text:
## TODO: roll into InducedMap work?
struct ComposedMarkovMap{T<:Tuple,D<:Domain,R<:Domain} <: AbstractMarkovMap{D,R}
  maps::T
  domain::D
  rangedomain::R
end
(∘)(f::AbstractMarkovMap,g::AbstractMarkovMap) = ComposedMarkovMap(f,g)

function ComposedMarkovMap(maps...)
  ran = rangedomain(maps[1])
  mps = isa(maps[1],ComposedMarkovMap) ? maps[1].maps : (maps[1],)
  for i = 1:length(maps)-1
    @assert rangedomain(maps[i+1]) == domain(maps[i])
    mps = isa(maps[i+1],ComposedMarkovMap) ? (mps...,maps[i+1].maps...) : (mps...,maps[i+1])
  end
  dom = domain(maps[end])

  ComposedMarkovMap(mps,dom,ran)
end

complength(c::ComposedMarkovMap) = length(c.maps)

function (c::ComposedMarkovMap)(x)
  f = c.maps[end](x)
  for i = complength(c)-1:-1:1
    f = c.maps[i](f)
  end
  f
end

function mapP(c::ComposedMarkovMap,x)
  f,dfdx = mapP(c.maps[end],x)
  for i = complength(c)-1:-1:1
    f,dfdxs = mapP(c.maps[i],f)
    dfdx *= dfdxs
  end
  f,dfdx
end
mapD(c::ComposedMarkovMap,x) = mapP(c,x)[2]

function mapinv(c::ComposedMarkovMap,b,x)
  v = mapinv(c.maps[1],b[1],x)
  for i = 2:complength(c)
    v = mapinv(c.maps[i],b[i],v)
  end
  v
end
function mapinvP(c::ComposedMarkovMap,b,x)
  v,dvdx = mapinvP(c.maps[1],b[1],x)
  for i = 2:complength(c)
    v,dvdxs = mapinvP(c.maps[i],b[i],v)
    dvdx *= dvdxs
  end
  v,dvdx
end
mapinvD(c::ComposedMarkovMap,b,x) = mapinvP(c,b,x)[2]

# TODO: map(P,D)(c,b,x)

function getbranchind(m::ComposedMarkovMap,x)
  x ∈ m.domain || error("DomainError: $x ∉ $(m.domain)")
  fx = x
  br = getbranchind(m.maps[end],fx)
  for i = complength(c)-1:-1:1
    fx = m.maps[i+1](fx)
    br = (getbranchind(m.maps[i],fx),br...)
  end
  br
end

nbranches(C::ComposedMarkovMap) = prod(nbranches(mm for mm in C.maps))
eachbranchindex(C::ComposedMarkovMap) = product(eachbranchindex(mm) for mm in C.maps)
branchindtype(C::ComposedMarkovMap) = Array{promote_type([branchindtype(mm) for mm in C.maps]...),1}
nneutral(C::ComposedMarkovMap) = 0 # assume we can't do this
neutralfixedpoints(C::ComposedMarkovMap) = []

#TODO: must be faster??
function transferfunction(x,m::ComposedMarkovMap{Tuple{M}},f) where {M<:AbstractMarkovMap}
  transferfunction(x,m.maps[1],f)
end

struct TransferCall{M<:AbstractMarkovMap,ff}
  m::M
  f::ff
end
# TransferCall(m,f) = TransferCall{typeof(m),typeof(f),T}(m,f)
(t::TransferCall)(x) = transferfunction(x,t.m,t.f)

function transferfunction(x,m::ComposedMarkovMap,f)
  m2 = ComposedMarkovMap(m.maps[2:end],domain(m.maps[end]),rangedomain(m.maps[2]))
  transferfunction(x,m.maps[1],#x->transferfunction(x,m2,f))
    TransferCall(m2,f))
end

# TODO: transferfunction_int

# TODO
  struct ComposedCircleMap{M<:AbstractMarkovMap,D<:Domain,R<:Domain} <: AbstractMarkovMap{D,R}
    maps::NTuple{M}
    domain::D
    rangedomain::R
    end

    # stuff you can do for both
ComposedMap{M,D<:Domain,R<:Domain} = Union{ComposedMarkovMap{M,D,R},ComposedCircleMap{M,D,R}}

for FUN in (:domain,:rangedomain)
  @eval ($FUN)(c::ComposedMarkovMap) = c.$FUN
end


text:
# module Bukdu

module Keyword

function get(dict::Dict, key::Symbol, default::Any)
    haskey(dict, key) ? dict[key] : default
end

end # module Bukdu.Keyword


text:
using Printf

mutable struct GameLayer <: Ranger.AbstractNode
    base::Nodes.NodeData
    transform::Nodes.TransformProperties{Float64}

    # Collection of nodes.
    children::Array{Ranger.AbstractNode,1}

    min::Geometry.Point{Float64}
    max::Geometry.Point{Float64}
    buc_min::Geometry.Point{Float64}
    buc_max::Geometry.Point{Float64}

    host_node::HostNode

    function GameLayer(world::Ranger.World, name::String, parent::Ranger.AbstractNode)
        o = new()

        o.base = Nodes.NodeData(name, parent, world)
        o.children = Array{Ranger.AbstractNode,1}[]
        o.transform = Nodes.TransformProperties{Float64}()

        o.min = Geometry.Point{Float64}()
        o.max = Geometry.Point{Float64}()
        o.buc_min = Geometry.Point{Float64}()
        o.buc_max = Geometry.Point{Float64}()

        o
    end
end

function build(layer::GameLayer, world::Ranger.World)
    hw = Float64(world.view_width) / 2.0
    hh = Float64(world.view_height) / 2.0

    # top-left
    Geometry.set!(layer.min, -hw, -hh)

    # bottom-right
    Geometry.set!(layer.max, hw, hh)

    layer.host_node = HostNode(world, "HostNode", layer)
    Nodes.set_scale!(layer.host_node, 150.0)
    push!(layer.children, layer.host_node);

end

# --------------------------------------------------------
# Timing
# --------------------------------------------------------
function Nodes.update(layer::GameLayer, dt::Float64)
end

# --------------------------------------------------------
# Rendering
# --------------------------------------------------------
function Nodes.draw(layer::GameLayer, context::Rendering.RenderContext)
    if Nodes.is_dirty(layer)
        Rendering.transform!(context, layer.min, layer.max, layer.buc_min, layer.buc_max)
        Nodes.set_dirty!(layer, false)
    end

    # Render background (i.e. a solid gray color)
    Rendering.set_draw_color(context, RangerGame.darkgray)
    Rendering.render_aa_rectangle(context, layer.buc_min, layer.buc_max, Rendering.FILLED)

    # Draw layer name
    Rendering.set_draw_color(context, RangerGame.white)
    Rendering.draw_text(context, 10, 10, layer.base.name, 2, 2, false)
end

# --------------------------------------------------------
# Life cycle events
# --------------------------------------------------------
function Nodes.enter_node(layer::GameLayer, man::Nodes.NodeManager)
    println("enter ", layer)

    # Register node as a timing target in order to receive updates
    Nodes.register_target(man, layer)

    # Register host so it too can get updates
    Nodes.register_target(man, layer.host_node)

    # Register for io events such as keyboard or mouse
    Nodes.register_event_target(man, layer)
end

function Nodes.exit_node(layer::GameLayer, man::Nodes.NodeManager)
    println("exit ", layer)
    Nodes.unregister_target(man, layer)
    Nodes.unregister_target(man, layer.host_node)

    Nodes.unregister_event_target(man, layer)
end

# --------------------------------------------------------
# Events
# --------------------------------------------------------

# --------------------------------------------------------
# Grouping
# --------------------------------------------------------
function Nodes.get_children(node::GameLayer)
    node.children
end

max_stars_repo_name: wdevore/Ranger-Julia-SDL

max_issues_repo_name: wdevore/Ranger-Julia-SDL

max_forks_repo_name: wdevore/Ranger-Julia-SDL


text:
# get the 'current' working directory
println(pwd())

# read the contents of the 'current' working directory - returns
# an array of file / folder names
println(readdir())

# create / open a file for writing (create it first if it does't exist) ...
function filecreate(file::String, text::String="This is some text")
    io = open(file, "w")
    write(io, text)
    close(io)
end
filecreate("testfile.txt", "Hello from somewhere out there")

# read a file
function fileread(file::String)
    io = open(file, "r")
    contents = read(io, String) # if you don't specify a type you get UInt8 'char' collection returned
    close(io)
    println(contents)
end
fileread("testfile.txt")

# add to a file
function fileappend(file::String, text::String)
    io = open(file, "a")
    write(io, text)
    close(io)
end
fileappend("testfile.txt", "\nHello from somewhere else")

# rename a file
function filerename(file1::String, file2::String)
    mv(file1, file2)
end
filerename("testfile.txt", "testfilenew.txt")

# delete a file
rm("testfilenew.txt")

max_stars_repo_name: rpillar/JuliaExamples

max_issues_repo_name: rpillar/JuliaExamples

max_forks_repo_name: rpillar/JuliaExamples


text:
@testset "Testing ConvConversion" begin
    # Test the Conv flux layer

    input_channels = (10, 16, 20, 64)
    output_channels = (10, 16, 20, 64)
    batch_sizes = (16,)
    backend = nGraph.Backend("CPU")

    for (ic, oc, bs) in Iterators.product(input_channels, output_channels, batch_sizes)
        C = CrossCor((3,3), ic => oc, Flux.relu)
        x = rand(Float32, 14, 14, ic, bs)

        E = C(x)
        f = nGraph.compile(backend, C, x)
        Z = parent(f())

        approx = isapprox(Z, E)
        color = approx ? :green : :red
        printstyled("$ic $oc $bs\n"; color = color)
        @test approx
    end
end

@testset begin
    m = Chain(
        Dense(28^2, 32, relu),
        Dense(32, 10)
    )

    # Construct a dummy input.
    x = rand(Float32, 28^2)

    expected = m(x)

    # nGraph compile
    backend = nGraph.Backend()
    f = nGraph.compile(backend, m, x)

    @test isapprox(expected, parent(f()))
end


text:
###
### Constructors
###
###
### Constructor methods for LongSequences.
###
### This file is a part of BioJulia.
### License is MIT: https://github.com/BioJulia/BioSequences.jl/blob/master/LICENSE.md

function LongSequence{A}(len::Integer) where {A<:Alphabet}
    if len < 0
        throw(ArgumentError("len must be non-negative"))
    end
    return LongSequence{A}(Vector{UInt64}(undef, seq_data_len(A, len)), convert(Int, len))
end

LongSequence(::Type{DNA}) = LongDNASeq()
LongSequence(::Type{RNA}) = LongRNASeq()
LongSequence(::Type{AminoAcid}) = LongAminoAcidSeq()
LongSequence(::Type{Char}) = LongCharSeq()

function LongSequence()
    return LongSequence{VoidAlphabet}(Vector{UInt64}(), 0)
end

function LongSequence{A}(seq::LongSequence{A}, part::UnitRange) where A
    return seq[part]
end

function LongSequence{A}(s::Union{String, SubString{String}}) where {A<:Alphabet}
    return LongSequence{A}(s, codetype(A()))
end

# Generic method for String/Substring.
function LongSequence{A}(s::Union{String, SubString{String}}, ::AlphabetCode) where {A<:Alphabet}
    len = length(s)
    seq = LongSequence{A}(len)
    return copyto!(seq, 1, s, 1, len)
end

function LongSequence{A}(s::Union{String, SubString{String}}, ::AsciiAlphabet) where {A<:Alphabet}
    v = GC.@preserve s unsafe_wrap(Vector{UInt8}, pointer(s), ncodeunits(s))
    seq = LongSequence{A}(length(v))
    return encode_chunks!(seq, 1, v, 1, length(v))
end

function LongSequence{A}(
        src::Union{AbstractString,AbstractVector},
        startpos::Integer=1,
        stoppos::Integer=length(src)) where {A<:Alphabet}
    len = stoppos - startpos + 1
    seq = LongSequence{A}(len)
    return copyto!(seq, 1, src, startpos, len)
end

# create a subsequence
function LongSequence(other::LongSequence, part::UnitRange{<:Integer})
    checkbounds(other, part)
    subseq = typeof(other)(length(part))
    copyto!(subseq, 1, other, first(part), length(part))
    return subseq
end

function (::Type{T})(seq::BioSequence) where {T<:LongSequence}
    newseq = T(length(seq))
    @inbounds for i in eachindex(seq)
        newseq[i] = seq[i]
    end
    return newseq
end

function LongSequence{A}(seq::LongSequence{A}) where {A <: Alphabet}
    return LongSequence{A}(copy(seq.data), seq.len)
end

function (::Type{T})(seq::LongSequence{<:NucleicAcidAlphabet{N}}) where
         {N, T<:LongSequence{<:NucleicAcidAlphabet{N}}}
    return T(copy(seq.data), seq.len)
end

# This exists to fix ambiguity errors with Julia 1.0
function LongSequence{A}(seq::LongSequence{<:NucleicAcidAlphabet{N}}) where {N, A <: NucleicAcidAlphabet{N}}
	return LongSequence{A}(copy(seq.data), seq.len)
end

# Concatenate multiple sequences
function LongSequence{A}(chunks::LongSequence{A}...) where {A}
    len = 0
    for chunk in chunks
        len += length(chunk)
    end
    seq = LongSequence{A}(len)
    offset = 1
    for chunk in chunks
        copyto!(seq, offset, chunk, 1, length(chunk))
        offset += length(chunk)
    end
    return seq
end

function Base.repeat(chunk::LongSequence{A}, n::Integer) where {A}
    seq = LongSequence{A}(length(chunk) * n)
    offset = 1
    for i in 1:n
        copyto!(seq, offset, chunk, 1, length(chunk))
        offset += length(chunk)
    end
    return seq
end

# Concatenation and Base.repeat operators
Base.:*(chunk::LongSequence{A}, chunks::LongSequence{A}...) where {A} =
    LongSequence{A}(chunk, chunks...)
Base.:^(chunk::LongSequence, n::Integer) = repeat(chunk, n)

function Base.similar(seq::LongSequence{A}, len::Integer = length(seq)) where {A}
    return LongSequence{A}(len)
end


text:
#!/cvmfs/soft.computecanada.ca/easybuild/software/2020/avx2/Core/julia/1.6.0/bin/julia
#SBATCH --mail-user=mmcleod2@ualberta.ca
#SBATCH --mail-type=ALL
<<<<<<< HEAD
#SBATCH -o TwodGridWorld_GPI.out # Standard output
#SBATCH -e TwodGridWorld_GPI.err # Standard error
#SBATCH --time=01:30:00
#SBATCH --nodes=1
=======
#SBATCH -o MC_GPI.out # Standard output
#SBATCH -e MC_GPI.err # Standard error
#SBATCH --time=03:00:00
#SBATCH --nodes=3
>>>>>>> mountain_car_dev
#SBATCH --ntasks-per-node=48
#SBATCH --mem=0
#SBATCH --account=def-amw8

using Pkg
Pkg.activate(".")

include(joinpath(ENV["SLURM_SUBMIT_DIR"], "parallel/parallel_config.jl"))
reproduce_config_experiment("configs/twodgridworld_hailmary/twodgridworld_gpi.toml", save_path="~/scratch/curiosity")

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
using FFTW
using Plots

a = 1        # Thermal diffusivity constant
L = 100      # Length of domain
N = 1000     # Number of discretization points
dx = L/N
x = -L/2:dx:L/2-dx # Define x domain

# Define discrete wavenumbers
kappa = (2pi/L)*(-N/2:N/2-1)
kappa = fftshift(kappa) # Re-order fft wavenumbers

# Initial condition
u0 = zeros(N)
u0[(L/2 - L/10)/dx:(L/2 + L/10)/dx] = 1

# Simulate in Fourier frequency domain
# t = 0:0.1:10
# [t,uhat]=ode45(@(t,uhat)rhsHeat(t,uhat,kappa,a),t,fft(u0));

# for k = 1:length(t) # iFFT to return to spatial domain
#     u(k,:) = ifft(uhat(k,:))
# end

# Plot solution in time
figure, waterfall((u(1:10:end,:)));
figure, imagesc(flipud(u));

%% FIGURES (PRODUCTION)
figure
CC = colormap(jet(100));
dt = 0.1;
for k = 1:100
    u(k,:) = ifft(uhat(k,:));
    if(mod(k-1,10)==0)
        plot(x,u(k,:),'Color',CC(k,:),'LineWidth',1.5)
        hold on, grid on, drawnow
    end
end
% xlabel('Spatial variable, x')
% ylabel('Temperature, u(x,t)')
axis([-50 50 -.1 1.1])
set(gca,'LineWidth',1.2,'FontSize',12);
set(gcf,'Position',[100 100 550 220]);
set(gcf,'PaperPositionMode','auto')
% print('-depsc2', '-loose', '../../figures/FFTHeat1');

%
figure
subplot(1,2,1)
h=waterfall(u(1:10:end,:));
set(h,'LineWidth',2,'FaceAlpha',.5);
colormap(jet/1.5)
view(22,29)
set(gca,'LineWidth',1.5)
set(gca,'XTick',[0 500 1000],'XTickLabels',{})
set(gca,'ZTick',[0 .5 1],'ZTickLabels',{})
set(gca,'YTick',[0 5 10],'YTickLabels',{})

subplot(1,2,2)
imagesc(flipud(u));
set(gca,'LineWidth',1.5)
set(gca,'XTick',[0 500 1000],'XTickLabels',{})
set(gca,'YTick',[0 50 100],'YTickLabels',{})

colormap jet
set(gcf,'Position',[100 100 600 250])
set(gcf,'PaperPositionMode','auto')
% print('-depsc2', '-loose', '../../figures/FFTHeat2');


text:
struct BinCovFactor{T<:Unsigned} <: AbstractFactor{T}
    name::String
    levels::AbstractVector{<:AbstractString}
    bins::AbstractVector{<:Real}
    covariate::AbstractCovariate{<:AbstractFloat}
end

Base.length(var::BinCovFactor{T}) where {T<:Unsigned} = length(var.covariate)

function BinCovFactor(name::String, bins::AbstractVector{T}, covariate::AbstractCovariate{S}) where {T<:Real} where {S<:AbstractFloat}
    bins = issorted(bins) ? bins : sort(bins)
    levelcount = length(bins)
    levels = [@sprintf("[%G,%G%s", bins[i], bins[i + 1], i == levelcount - 1 ? "]" : ")") for i in 1:levelcount - 1]
    if levelcount <= typemax(UInt8)
        BinCovFactor{UInt8}(name, levels, bins, covariate)
    elseif levelcount <= typemax(UInt16)
        BinCovFactor{UInt16}(name, levels, bins, covariate)
    else
        BinCovFactor{UInt32}(name, levels, bins, covariate)
    end
end

function factor(covariate::AbstractCovariate{S}, bins::AbstractVector{T}) where {T<:Real} where {S<:AbstractFloat}
    BinCovFactor(getname(covariate), bins, covariate)
end

function slice(factor::BinCovFactor{T}, fromobs::Integer, toobs::Integer, slicelength::Integer) where {T<:Unsigned}
    bins = factor.bins
    binlenm1 = convert(T, length(bins) - 1)
    z = bins[length(bins)]
    f = (x ->
            begin
                if x == z
                    return binlenm1
                elseif x > z
                    return 0
                else
                    i = searchsortedlast(bins, x)
                    return convert(T, i)
                end
            end
        )
    slicelength = verifyslicelength(fromobs, toobs, slicelength)
    slices = slice(factor.covariate, fromobs, toobs, slicelength)
    mapslice(f, slices, slicelength, T)
end

function Base.map(factor::BinCovFactor{T}, dataframe::AbstractDataFrame) where {T<:Unsigned}
    BinCovFactor(getname(factor), factor.bins, map(factor.covariate, dataframe))
end

function isordinal(factor::BinCovFactor{T}) where {T<:Unsigned}
    true
end

max_stars_repo_name: vikashdwd/JuML

max_issues_repo_name: vikashdwd/JuML

max_forks_repo_name: vikashdwd/JuML


text:
using DimReduc
using Documenter

makedocs(;
    modules=[DimReduc],
    authors="Alex",
    repo="https://github.com/Student2Pro/DimReduc.jl/blob/{commit}{path}#L{line}",
    sitename="DimReduc.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://Student2Pro.github.io/DimReduc.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/Student2Pro/DimReduc.jl",
)


text:
function mgrid(ngrid)

# evaluates the grid in the meridional coordinate y=cos(theta)
# and also returns the corresponding values of latitude. The
# end points at the equator and the pole are excluded.
#
#
# input
# ngrid  - number of grid points
#
# output
# y  - grid values
# lat - grid values in latitude (radian)
#

# Evaluates the grid
dy = 1.0 / (ngrid+1);
y = collect(dy : dy : 1.0-dy);

# corresponding latitude (degree)
lat = asin.(y)*180/pi;

return y, lat

end

max_stars_repo_name: bbuffett/MagRossby

max_issues_repo_name: bbuffett/MagRossby

max_forks_repo_name: bbuffett/MagRossby


text:
function makePlots()

println("Generating plots")

plotsDir = (pwd()[end-3:end] == "docs") ? "build/plots" : "docs/build/plots"
mkdir(plotsDir)
Plots.pyplot()

# PID design functions
P = tf(1,[1,1])^4
gangoffourplot(P,tf(1))
Plots.savefig(plotsDir*"/pidgofplot.svg")

ωp = 0.8
kp,ki,C = loopshapingPI(P,ωp,phasemargin=60, doplot=false)
gangoffourplot(P, [tf(1), C])
Plots.savefig(plotsDir*"/pidgofplot2.svg")
nyquistplot([P, P*C])
Plots.savefig(plotsDir*"/pidnyquistplot.svg")

ωp = 2
kp,ki,C60 = loopshapingPI(P,ωp,rl=1,phasemargin=60, doplot=true)
gangoffourplot(P, [tf(1), C60])
Plots.savefig(plotsDir*"/pidgofplot3.svg")
nyquistplot([P, P*C60])
Plots.savefig(plotsDir*"/pidnyquistplot2.svg")

# Advanced pole placement
ζ = 0.2
ω = 1
B = [1]
A   = [1, 2ζ*ω, ω^2]
P  = tf(B,A)
# Control design
ζ0 = 0.7
ω0 = 2
Am = [1, 2ζ0*ω0, ω0^2]
Ao = conv(2Am, [1/2, 1]) # Observer polynomial
AR = [1,0] # Force the controller to contain an integrator

B⁺  = [1] # The process numerator polynomial can be facored as B = B⁺B⁻ where B⁻ contains the zeros we do not want to cancel (non-minimum phase and poorly damped zeros)
B⁻  = [1]
Bm  = conv(B⁺, B⁻) # In this case, keep the entire numerator polynomial of the process

R,S,T = rstc(B⁺,B⁻,A,Bm,Am,Ao,AR) # Calculate the 2-DOF controller polynomials

Gcl = tf(conv(B,T),zpconv(A,R,B,S)) # Form the closed loop polynomial from reference to output

stepplot([P,Gcl]) # Visualize the open and closed loop responses.
Plots.savefig(plotsDir*"/ppstepplot.svg")
gangoffourplot(P, tf(-S,R)) # Plot the gang of four to check that all tranfer functions are OK
Plots.savefig(plotsDir*"/ppgofplot.svg")

P1 = "exp(-sqrt(s))"
f1 = stabregionPID(P1,logspace(-5,1,1000)); Plots.savefig(plotsDir*"/stab1.svg")
P2 = "100*(s+6).^2./(s.*(s+1).^2.*(s+50).^2)"
f2 = stabregionPID(P2,logspace(-5,2,1000)); Plots.savefig(plotsDir*"/stab2.svg")
P3 = tf(1,[1,1])^4
f3 = stabregionPID(P3,logspace(-5,0,1000)); Plots.savefig(plotsDir*"/stab3.svg")

# PID plots

P = tf([1.],[1., 1])
ζ = 0.5 # Desired damping
ws = logspace(-1,2,8) # A vector of closed-loop bandwidths
kp = 2*ζ*ws-1 # Simple pole placement with PI given the closed-loop bandwidth, the poles are placed in a butterworth pattern
ki = ws.^2
pidplots(P,:nyquist,;kps=kp,kis=ki, ω= logspace(-2,2,500))
Plots.savefig(plotsDir*"/pidplotsnyquist1.svg")
pidplots(P,:gof,;kps=kp,kis=ki, ω= logspace(-2,2,500))
Plots.savefig(plotsDir*"/pidplotgof1.svg")

kp = linspace(-1,1,8) # Now try a different strategy, where we have specified a gain crossover frequency of 0.1 rad/s
ki = sqrt(1-kp.^2)/10
pidplots(P,:nyquist,;kps=kp,kis=ki)
Plots.savefig(plotsDir*"/pidplotsnyquist2.svg")
pidplots(P,:gof,;kps=kp,kis=ki)
Plots.savefig(plotsDir*"/pidplotsgof2.svg")

#AR = [1,0]
#B⁺  = [1]
#B⁻  = [1]
#Bm  = conv(B⁺, B⁻)
#R,S,T = rstc(B⁺,B⁻,A,bm,Am,Ao,AR)
#Gcl = tf(conv(B,T),zpconv(A,R,B,S))
#f1 = stepplot([P,Gcl])
#f2 = gangoffourplot(P, tf(-S,R))
#Plots.savefig(f1, "$plotsDir/rstcstepplot.svg")
#Plots.savefig(f2, "$plotsDir/rstcgofplot.svg")
end


text:
" Implements a Convection operator "
type Convection <: Kernel
    u::Variable

    velocity::Vec{2, Float64}
end

@inline function computeQpResidual(kernel::Convection, qp::Integer, i::Integer)
    u = kernel.u

    return (kernel.velocity ⋅ u.grad[qp]) * u.phi[qp][i]
end

@inline function computeQpJacobian(kernel::Convection, v::Variable, qp::Integer, i::Integer, j::Integer)::Float64
    u = kernel.u

    if u.id == v.id
        return (kernel.velocity ⋅ v.grad_phi[qp][j]) * u.phi[qp][i]
    end

    return 0
end

@inline function coupledVars(kernel::Convection)::Array{Variable}
    return [kernel.u]
end


text:
include("data.jl")

mutable struct supportsresistances
   close::Float64
   date::
   timeframe::String
   asset::String
   regime::String
   pos::String
   supportsresistances(close, date, asset, timeframe) = begin
      close = close
      date = date
      timeframe = timeframe.upper()
      asset = asset
      regime = asset.regime
      pos = length(asset.data[timeframe])
   new(close, date, asset, timeframe)
   end
end

function RSI(Market::market)
   position = length(Market.data)
   if position < 15
      return 0.00
   else
      data = Market.data[position-14:position]
      wins = []
      losses = []
      for i in data
         if i.Close < i.Open
            push!(losses,(i.Open-i.Close))
         elseif i.Close > i.Open
            push!(wins,(i.Close-i.Open))
         end
      end
      winssum = sum(wins)/length(wins)
      lossessum = sum(losses)/length(losses)
      firstrs = winssum/lossessum
      RSI = 100-(100/firstrs)
      return RSI
   end
end

function MACD(ma26::Float64,ma12::Float64,ma9::Float64)
   macdline = ma12 - ma26
   signalline = ma9
   macdhistogram = macdline - signalline
   return macdhistogram
end

function MA9(Market::market)
   position = length(Market.data)
   if position < 10
      return 0.00
   else
      ma9 = sum(Market.data[position-9].Close:Market.data[position].Close)/9
      return ma9
   end
end

function MA12(Market::market)
   position = length(Market.data)
   if position < 13
      return 0.00
   else
      ma12 = sum(Market.data[position-12].Close:Market.data[position].Close)/12
      return ma12
   end
end

function MA26(Market::market)
   position = length(Market.data)
   if position < 27
      return 0.00
   else
      ma26 = sum(Market.data[position-26].Close:Market.data[position].Close)/26
      return ma26
   end
end

function MA52(Market::market)
   position = length(Market.data)
   if position < 53
      return 0.00
   else
      ma52 = sum(Market.data[position-52].Close:Market.data[position].Close)/52
      return ma52
   end
end

max_stars_repo_name: SRelieved/JuliaBacktesting

max_issues_repo_name: MartinGenesse/JuliaBacktesting

max_forks_repo_name: MartinGenesse/JuliaBacktesting


text:
__precompile__()
module SExpressions

using Reexport

include("Lists/Lists.jl")
include("Keywords/Keywords.jl")
include("Parser/Parser.jl")

export SExpression
import .Parser: parse, parseall, parsefile, SExpression
@reexport using .Lists

macro sx_str(x::String)
    QuoteNode(parse(x))
end

export @sx_str

end  # module SExpressions


text:
using ApproxFun, OrdinaryDiffEq, Sundials
using DiffEqDevTools
using LinearAlgebra
using Plots; gr()

function cheb(N)
    N==0 && return (0,1)
    x = cos.(pi*(0:N)/N)
    c = [2; ones(N-1,1); 2].*(-1).^(0:N)
    X = hcat([x for i in 1:N+1]...)
    dX = X-X'
    D  = (c*(1 ./c)')./(dX+I)      # off-diagonal entries
    D  = D .- Diagonal(vec(sum(D,dims=2)))                 # diagonal entries
    D,x
end
N = 128
ChebD2,x = cheb(N)
xx = x
x = x[2:N]
w = .53*x + .47*sin.(-1.5*pi*x) - x # use w = u-x to make BCs homogeneous
u = [1;w+x;-1]

ϵ=0.01
D2=ϵ*(ChebD2^2)[2:N, 2:N]
function allen_cahn(du,u,x,t)
    @. du = (u + x) - (u + x)^3
end

prob = SplitODEProblem(DiffEqArrayOperator(D2), allen_cahn, w, (0.0,5.0), x)
sol  = solve(prob, RadauIIA5(autodiff=false); reltol=1e-14,abstol=1e-14)
test_sol = TestSolution(sol)

tslices=[0.0 1.0 2.0 3.0 5.0]
ys=hcat(([1;x.+sol(t);-1] for t in tslices)...)
labels=["t=$t" for t in tslices]
plot(xx,ys,label=labels)

abstols = 0.1 .^ (5:8)
reltols = 0.1 .^ (1:4)
multipliers =  0.5 .^ (0:3)
setups = [Dict(:alg => IMEXEuler(), :dts => 1e-3 * multipliers),
          Dict(:alg => CNAB2(), :dts => 5e-3 * multipliers),
          Dict(:alg => CNLF2(), :dts => 5e-4 * multipliers),
          Dict(:alg => SBDF2(), :dts => 1e-3 * multipliers)]
labels = ["IMEXEuler" "CNAB2" "CNLF2" "SBDF2"]
@time wp1 = WorkPrecisionSet(prob,abstols,reltols,setups;
                            print_names=true,names=labels,
                            numruns=5,seconds=5,
                            save_everystop=false,appxsol=test_sol,maxiters=Int(1e5));

plot(wp1,label=labels,markershape=:auto,title="IMEX methods, dense linsolve, low order")

abstols = 0.1 .^ (5:8) # all fixed dt methods so these don't matter much
reltols = 0.1 .^ (1:4)
multipliers = 0.5 .^ (0:3)
setups = [Dict(:alg => IMEXEuler(linsolve=LinSolveGMRES()), :dts => 1e-3 * multipliers),
          Dict(:alg => CNAB2(linsolve=LinSolveGMRES()), :dts => 5e-3 * multipliers),
          Dict(:alg => CNLF2(linsolve=LinSolveGMRES()), :dts => 5e-4 * multipliers),
          Dict(:alg => SBDF2(linsolve=LinSolveGMRES()), :dts => 1e-3 * multipliers)]
labels = ["IMEXEuler" "CNAB2" "CNLF2" "SBDF2"]
@time wp1 = WorkPrecisionSet(prob,abstols,reltols,setups;
                            print_names=true, names=labels,
                            numruns=5, error_estimate=:l2,
                            save_everystep=false, appxsol=test_sol, maxiters=Int(1e5));

plot(wp1, label=labels, markershape=:auto, title="IMEX methods, Krylov linsolve, low order")

abstols = 0.1 .^ (5:8) # all fixed dt methods so these don't matter much
reltols = 0.1 .^ (1:4)
multipliers = 0.5 .^ (0:3)
setups = [Dict(:alg => NorsettEuler(), :dts => 1e-3 * multipliers),
          Dict(:alg => NorsettEuler(krylov=true, m=5), :dts => 1e-3 * multipliers),
          Dict(:alg => NorsettEuler(krylov=true, m=20), :dts => 1e-3 * multipliers),
          Dict(:alg => ETDRK2(), :dts => 1e-2 * multipliers),
          Dict(:alg => ETDRK2(krylov=true, m=5), :dts => 1e-2 * multipliers),
          Dict(:alg => ETDRK2(krylov=true, m=20), :dts => 1e-2 * multipliers)]
labels = hcat("NorsettEuler (caching)", "NorsettEuler (m=5)", "NorsettEuler (m=20)",
                            "ETDRK2 (caching)", "ETDRK2 (m=5)", "ETDRK2 (m=20)")
@time wp2 = WorkPrecisionSet(prob,abstols,reltols,setups;
                            print_names=true, names=labels,
                            numruns=5,
                            save_everystep=false, appxsol=test_sol, maxiters=Int(1e5));

plot(wp2, label=labels, markershape=:auto, title="ExpRK methods, low order")

abstols = 0.1 .^ (5:8) # all fixed dt methods so these don't matter much
reltols = 0.1 .^ (1:4)
multipliers = 0.5 .^ (0:3)
setups = [Dict(:alg => CNAB2(), :dts => 5e-3 * multipliers),
          Dict(:alg => CNAB2(linsolve=LinSolveGMRES()), :dts => 5e-3 * multipliers),
          Dict(:alg => ETDRK2(krylov=true, m=5), :dts => 1e-2 * multipliers)]
labels = ["CNAB2 (dense linsolve)" "CNAB2 (Krylov linsolve)" "ETDRK2 (m=5)"]
@time wp3 = WorkPrecisionSet(prob,abstols,reltols,setups;
                            print_names=true, names=labels,
                            numruns=5, error_estimate=:l2,
                            save_everystep=false, appxsol=test_sol, maxiters=Int(1e5));

plot(wp3, label=labels, markershape=:auto, title="Between family, low orders")

abstols = 0.1 .^ (7:13)
reltols = 0.1 .^ (4:10)
setups = [Dict(:alg => KenCarp3()),
          Dict(:alg => KenCarp4()),
          Dict(:alg => KenCarp5()),
          Dict(:alg => ARKODE(Sundials.Implicit(), order=3, linear_solver=:Band, jac_upper=1, jac_lower=1)),
          Dict(:alg => ARKODE(Sundials.Implicit(), order=4, linear_solver=:Band, jac_upper=1, jac_lower=1)),
          Dict(:alg => ARKODE(Sundials.Implicit(), order=5, linear_solver=:Band, jac_upper=1, jac_lower=1))]
labels = hcat("KenCarp3", "KenCarp4", "KenCarp5",
              "ARKODE3", "ARKODE4", "ARKODE5")
@time wp4 = WorkPrecisionSet(prob,abstols,reltols,setups;
                            print_names=true, names=labels,
                            numruns=5, error_estimate=:l2,
                            save_everystep=false, appxsol=test_sol, maxiters=Int(1e5));

plot(wp4, label=labels, markershape=:auto, title="IMEX methods, dense/band linsolve, medium order")

abstols = 0.1 .^ (7:13)
reltols = 0.1 .^ (4:10)
setups = [Dict(:alg => KenCarp3(linsolve=LinSolveGMRES())),
          Dict(:alg => KenCarp4(linsolve=LinSolveGMRES())),
          Dict(:alg => KenCarp5(linsolve=LinSolveGMRES())),
          Dict(:alg => ARKODE(Sundials.Implicit(), order=3, linear_solver=:GMRES)),
          Dict(:alg => ARKODE(Sundials.Implicit(), order=4, linear_solver=:GMRES)),
          Dict(:alg => ARKODE(Sundials.Implicit(), order=5, linear_solver=:GMRES))]
labels = ["KenCarp3" "KenCarp4" "KenCarp5" "ARKODE3" "ARKODE4" "ARKODE5"]
@time wp4 = WorkPrecisionSet(prob,abstols,reltols,setups;
                            print_names=true, names=labels,
                            numruns=5, error_estimate=:l2,
                            save_everystep=false, appxsol=test_sol, maxiters=Int(1e5));

plot(wp4, label=labels, markershape=:auto, title="IMEX methods, Krylov linsolve, medium order")

abstols = 0.1 .^ (7:11) # all fixed dt methods so these don't matter much
reltols = 0.1 .^ (4:8)
multipliers = 0.5 .^ (0:4)
setups = [Dict(:alg => ETDRK3(), :dts => 1e-2 * multipliers),
          Dict(:alg => ETDRK3(krylov=true, m=5), :dts => 1e-2 * multipliers),
          Dict(:alg => ETDRK4(), :dts => 1e-2 * multipliers),
          Dict(:alg => ETDRK4(krylov=true, m=5), :dts => 1e-2 * multipliers),
          Dict(:alg => HochOst4(), :dts => 1e-2 * multipliers),
          Dict(:alg => HochOst4(krylov=true, m=5), :dts => 1e-2 * multipliers)]
labels = hcat("ETDRK3 (caching)", "ETDRK3 (m=5)", "ETDRK4 (caching)",
              "ETDRK4 (m=5)", "HochOst4 (caching)", "HochOst4 (m=5)")
@time wp5 = WorkPrecisionSet(prob,abstols,reltols,setups;
                            print_names=true, names=labels,
                            numruns=5, error_estimate=:l2,
                            save_everystep=false, appxsol=test_sol, maxiters=Int(1e5));

plot(wp5, label=labels, markershape=:auto, title="ExpRK methods, medium order")

abstols = 0.1 .^ (7:11)
reltols = 0.1 .^ (4:8)
multipliers = 0.5 .^ (0:4)
setups = [Dict(:alg => KenCarp5()),
          Dict(:alg => ARKODE(Sundials.Implicit(), order=5, linear_solver=:Dense)),
          Dict(:alg => KenCarp5(linsolve=LinSolveGMRES())),
          Dict(:alg => ARKODE(Sundials.Implicit(), order=5, linear_solver=:GMRES)),
          Dict(:alg => ETDRK3(krylov=true, m=5), :dts => 1e-2 * multipliers),
          Dict(:alg => ETDRK4(krylov=true, m=5), :dts => 1e-2 * multipliers)]
labels = hcat("KenCarp5 (dense linsolve)", "ARKODE (dense linsolve)", "KenCarp5 (Krylov linsolve)",
              "ARKODE (Krylov linsolve)", "ETDRK3 (m=5)", "ETDRK4 (m=5)")
@time wp6 = WorkPrecisionSet(prob,abstols,reltols,setups;
                            print_names=true, names=labels,
                            numruns=5, error_estimate=:l2,
                            save_everystep=false, appxsol=test_sol, maxiters=Int(1e5));

plot(wp6, label=labels, markershape=:auto, title="Between family, medium order")

using DiffEqBenchmarks
DiffEqBenchmarks.bench_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])


text:
using Base.Test
using DataFrames

@testset "All spins down" begin
    const L = 4
    const meas = 10
    const theT = 1

    ferrofileA = tempname()
    ferrofileB = tempname()
    run(`./exe/annealer -j $(rand(UInt32)) -f -l $L -i 1e3 -m $meas $theT $ferrofileA`)
    run(`./exe/annealer -j $(rand(UInt32)) -f -l $L -i 1e3 -m $meas $theT $ferrofileB`)
    const ONES_BY_TRIPLETS = UInt64(0x1249249249249249)
    ONES = [ONES_BY_TRIPLETS for i in 1:L*L*L] :: Vector{UInt64}

    # Overwrite all the spins with "↑"
    fakefileA = tempname()
    fakefileB = tempname()

    # Copy the header
    ioferroA = open(ferrofileA,"r")
    iofakeA  = open(fakefileA,"w")

    ioferroB = open(ferrofileB,"r")
    iofakeB  = open(fakefileB,"w")

    header_A        = read(ioferroA,Int8,8)
    is_big_endian_A = read(ioferroA,Int64)
    rL_A            = read(ioferroA,Int64)
    ntws_A          = read(ioferroA,Int64)
    tws_A           = read(ioferroA,Int64,ntws_A)
    nts_A           = read(ioferroA,Int64)
    ts_A            = read(ioferroA,Int64,nts_A)
    total_meas_A    = read(ioferroA,Int64)
    T_A             = read(ioferroA,Float64)

    header_B        = read(ioferroB,Int8,8)
    is_big_endian_B = read(ioferroB,Int64)
    rL_B            = read(ioferroB,Int64)
    ntws_B          = read(ioferroB,Int64)
    tws_B           = read(ioferroB,Int64,ntws_B)
    nts_B           = read(ioferroB,Int64)
    ts_B            = read(ioferroB,Int64,nts_B)
    total_meas_B    = read(ioferroB,Int64)
    T_B             = read(ioferroB,Float64)

    write(iofakeA, header_A)
    write(iofakeA, is_big_endian_A)
    write(iofakeA, rL_A)
    write(iofakeA, ntws_A)
    write(iofakeA, tws_A)
    write(iofakeA, nts_A)
    write(iofakeA, ts_A)
    write(iofakeA, total_meas_A)
    write(iofakeA, T_A)

    write(iofakeB, header_B)
    write(iofakeB, is_big_endian_B)
    write(iofakeB, rL_B)
    write(iofakeB, ntws_B)
    write(iofakeB, tws_B)
    write(iofakeB, nts_B)
    write(iofakeB, ts_B)
    write(iofakeB, total_meas_B)
    write(iofakeB, T_B)

    for i in 1:3 # J's
        Js_A = read(ioferroA, UInt64, L*L*L)
        Js_B = read(ioferroB, UInt64, L*L*L)
        write(iofakeA, Js_A)
        write(iofakeB, Js_B)
    end

    # Now, for every measure set, copy it with spins down except for
    # the first iteration.
    flagfirst = true
    for i in 1:total_meas_A
        mcA = read(ioferroA,Int64)
        write(iofakeA, mcA)

        mcB = read(ioferroB,Int64)
        write(iofakeB, mcB)

        spinsA = read(ioferroA,UInt64,L*L*L)
        spinsB = read(ioferroB,UInt64,L*L*L)
        if flagfirst
            write(iofakeA, ONES)
            write(iofakeB, ONES)
            flagfirst = false
        else
            write(iofakeA, zeros(UInt64,L*L*L))
            write(iofakeB, zeros(UInt64,L*L*L))
        end
    end

    footerA = read(ioferroA,Int8,3)
    write(iofakeA,footerA)

    footerB = read(ioferroB,Int8,3)
    write(iofakeB,footerB)

    close(ioferroA); close(iofakeA)
    close(ioferroB); close(iofakeB)

    # Measure things
    measfile = tempname()

    run(pipeline(`./exe/measure -cCemq $fakefileA $fakefileB`,stdout=measfile))
    data = readtable(measfile)

    # The first iteration is completelly uncorrelated with everything
    # else except itself.
    # The rest are completelly correlated.
    Cdata = data[data[:observable] .== "Correlation",:]
    for i in 1:size(Cdata,1)
        mc = Cdata[:mc][i]
        t  = Cdata[:parameter][i]
        C  = Cdata[:value][i]
        if mc == 1 && t≠0
            @test C ≈ -1
        else
            @test C ≈ 1
        end
    end

    # The spatial correlation is complete.
    sCdata = data[data[:observable] .== "Spatial_correlation",:]
    @test all(sCdata[:value] .== 1)

    # The energy is always -1.
    Edata = data[data[:observable] .== "Energy",:]
    @test all(Edata[:value] .== -1)

    # The magnetization is always -1 (except for the first iteration).
    Mdata = data[data[:observable] .== "Magnetization",:]
    for i in 1:size(Mdata,1)
        if Mdata[:mc][i] == 1
            @test Mdata[:value][i] == 1
        else
            @test Mdata[:value][i] == -1
        end
    end

    # Overlap is absolute.
    Qdata = data[data[:observable] .== "Overlap",:]
    @test all(Qdata[:value] .== 1)

    # The only ts are the intended ones.
    U_ts  = Cdata[:parameter] |> unique
    @test U_ts == vcat(ts_A)
    @test U_ts == vcat(ts_B)

    # mc = tw + t
    tw = data[:mc]
    tw_plus_t = [tw+t for tw in tws_A, t in vcat(0,ts_A)]
    tw_plus_t = [tw+t for tw in tws_B, t in vcat(0,ts_B)]
    @test tw |> unique |> sort == tw_plus_t |> unique |> sort
end

max_stars_repo_name: RedPointyJackson/tfg

max_issues_repo_name: RedPointyJackson/tfg

max_forks_repo_name: RedPointyJackson/tfg


text:
module TestMOIWrapper

using Test
using MathOptInterface
import Cbc

const MOI = MathOptInterface

function runtests()
    for name in names(@__MODULE__; all = true)
        if startswith("$(name)", "test_")
            @testset "$(name)" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
    return
end

function test_SolverName()
    @test MOI.get(Cbc.Optimizer(), MOI.SolverName()) ==
          "COIN Branch-and-Cut (Cbc)"
end

function test_supports_incremental_interface()
    @test !MOI.supports_incremental_interface(Cbc.Optimizer())
    return
end

function test_runtests()
    model = MOI.Bridges.full_bridge_optimizer(
        MOI.Utilities.CachingOptimizer(
            MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
            Cbc.Optimizer(),
        ),
        Float64,
    )
    MOI.set(model, MOI.Silent(), true)
    MOI.Test.runtests(
        model,
        MOI.Test.Config(
            exclude = Any[
                MOI.ConstraintDual,
                MOI.DualObjectiveValue,
                MOI.ConstraintBasisStatus,
                MOI.VariableBasisStatus,
            ],
        ),
        exclude = [
            # TODO(odow): bug in Cbc.jl
            "test_model_copy_to_UnsupportedAttribute",
            "test_model_ModelFilter_AbstractConstraintAttribute",
            # TODO(odow): bug in MOI
            "test_model_LowerBoundAlreadySet",
            "test_model_UpperBoundAlreadySet",
            # TODO(odow): upstream bug in Cbc
            "_Indicator_",
            "_SOS1_",
            "_SOS2_",
            # Can't prove infeasible.
            "test_conic_NormInfinityCone_INFEASIBLE",
            "test_conic_NormOneCone_INFEASIBLE",
            "test_solve_TerminationStatus_DUAL_INFEASIBLE",
        ],
    )
    return
end

function test_params()
    # Note: we generate a non-trivial problem to ensure that Cbc struggles to
    # find a solution at the root node.
    knapsack_model = MOI.Utilities.Model{Float64}()
    N = 100
    x = MOI.add_variables(knapsack_model, N)
    MOI.add_constraint.(knapsack_model, x, MOI.ZeroOne())
    MOI.add_constraint(
        knapsack_model,
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.([1 + sin(i) for i in 1:N], x),
            0.0,
        ),
        MOI.LessThan(10.0),
    )
    MOI.set(
        knapsack_model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.([cos(i) for i in 1:N], x),
            0.0,
        ),
    )
    model = Cbc.Optimizer(
        maxSol = 1,
        presolve = "off",
        cuts = "off",
        heur = "off",
        logLevel = 0,
    )
    MOI.copy_to(model, knapsack_model)
    MOI.optimize!(model)
    @test MOI.get(model, MOI.TerminationStatus()) == MOI.SOLUTION_LIMIT
    @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
    MOI.empty!(model)
    MOI.copy_to(model, knapsack_model)
    MOI.optimize!(model)
    @test MOI.get(model, MOI.TerminationStatus()) == MOI.SOLUTION_LIMIT
    @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
    return
end

function test_PrimalStatus()
    model = MOI.Utilities.Model{Float64}()
    x = MOI.add_variable(model)
    MOI.add_constraint(model, x, MOI.GreaterThan(1.0))
    MOI.add_constraint(model, x, MOI.LessThan(0.0))
    cbc = Cbc.Optimizer()
    MOI.copy_to(cbc, model)
    MOI.optimize!(cbc)
    MOI.get(cbc, MOI.PrimalStatus()) == MOI.NO_SOLUTION
    return
end

end

TestMOIWrapper.runtests()


text:
module FieldDistributionNonuniformMedium

using ProgressMeter

include("simulation.jl")
include("plot.jl")

end


text:
module C2Julia
include("CNumbers.jl")
# `float` is also exported, which conflicts with `Base.float`. Resolve it with `const float = C2Julia.float`
export bool, char, unsigned_char, short, unsigned_short, int, unsigned_int, long, unsigned_long, long_long, unsigned_long_long, float, double
export @post_inc, @post_dec

include("CFunction.jl")
export @C, @bool, @cfor

include("CPointer.jl")
export @pointer, Pointer, malloc, free

# Export Base symbols
export +, -, *, /, %, &, |, ⊻, <<, >>, ~, !
export <, >, <=, >=, ==, !=
export sizeof

include("stdlib.jl")

end # module


text:
using Random
using DelimitedFiles
using Statistics

# Tested under Julia 1.6.2

function scale(center::Matrix{Float64}, point::Matrix{Float64}, scaling::Float64)
    @assert size(center) == size(point)
    @assert size(center, 1) == 1
    return scaling .* (point .- center) .+ center
end

# the code assumes
# - node id in embedding can be ignored and is a first column
# - communities are in the same order as nodes
# - files are space delimited
# - center of the community is computed using mean
function scale_communities(community_file::AbstractString,
                           embedding_file::AbstractString,
                           out_embedding_file::AbstractString,
                           scaling::Float64)
    communities = readdlm(community_file, Int)
    embedding = readdlm(embedding_file, Float64)[:, 2:end]

    @assert size(communities, 2) == 1
    @assert size(embedding, 1) == size(communities, 1)

    comm_ids = Dict{Int, Vector{Int}}()

    for (i, com) in enumerate(communities)
        push!(get!(comm_ids, com, Int[]), i)
    end

    @assert sum(length.(values(comm_ids))) == size(embedding, 1)
    @assert sort!(reduce(vcat, values(comm_ids))) == axes(embedding, 1)

    new_embedding = fill(NaN, size(embedding))
    for ids in values(comm_ids)
        center = mean(embedding[ids, :], dims=1)
        for i in ids
            @assert all(isnan, new_embedding[i, :])
            new_embedding[i:i, :] = scale(center, embedding[i:i, :], scaling)
            @assert !any(isnan, new_embedding[i, :])
        end
        @assert center ≈ mean(new_embedding[ids, :], dims=1)
        if scaling == 1.0
            @assert embedding[ids, :] ≈ new_embedding[ids, :]
        end
        if scaling == 0.0
            @assert repeat(center, length(ids)) ≈ new_embedding[ids, :]
        end
    end
    @assert !any(isnan, new_embedding)
    open(out_embedding_file, "w") do io
        for (i, row) in enumerate(eachrow(new_embedding))
            println(io, i-1, " ", join(row, " "))
        end
    end
    if scaling == 1
        @assert readdlm(embedding_file, Float64)[:, 2:end] ≈
                readdlm(out_embedding_file, Float64)[:, 2:end]
    end
    return nothing
end

scale_communities(ARGS[1],ARGS[2],ARGS[3], parse(Float64,ARGS[4]))


text:
using LinearAlgebra
using SparseArrays
using Plots

"This code solves the 2D Laplace's equation with nonzero Dirichlet and Neumann
boundary conditions using a finite difference method."

m = 100 # number of points
x = LinRange(-1,1,m+2)
y = LinRange(-1,1,m+2)
xint = x[1:end]
yint = y[2:end-1]
h = x[2]-x[1]

id(i,j,m) = i + (j-1)*(m+2) # return global index from i,j

function FD_matrix_2D(h,m)
    A = spzeros((m+2)*m,(m+2)*m)
    for i = 1:m+2, j = 1:m # loop thru x,y indices
        A[id(i,j,m),id(i,j,m)] = 4/h^2

        # stencil at interior nodes
        if i > 1 # avoids leftmost line of nodes
            A[id(i,j,m),id(i-1,j,m)] = -1/h^2 # x-derivative
        end
        if i < m+2 # avoids rightmost line of nodes
            A[id(i,j,m),id(i+1,j,m)] = -1/h^2
        end
        if j > 1
            A[id(i,j,m),id(i,j-1,m)] = -1/h^2 # y-derivative
        end
        if j < m
            A[id(i,j,m),id(i,j+1,m)] = -1/h^2
        end

        # modify for Neumann BCs on left/right
        # 2*(u_{m+1,j} - u_{m,j})/h^2 at Neumann boundaries in the x coordinate
        if i==1 # on left
            A[id(i,j,m),id(i+1,j,m)] = -2/h^2
        end
        if i==m+2 # on right
            A[id(i,j,m),id(i-1,j,m)] = -2/h^2
        end
    end
    return A
end

A = FD_matrix_2D(h,m)

# f(x,y) = exp(-10*((x-.25)^2 + (y-.75)^2))
f(x,y) = 0.0
b = vec(f.(xint,yint'))

# impose non-zero Dirichlet BCs
u_top(x) = 1.0
u_bottom(x) = -1.0
dudx_left(y) = y^2
dudx_right(y) = 20*y

for i = 1:m+2, j = 1:m
    # modify Neumann BCs first, then add Dirichlet contributions next
    if i==1 # left
        b[id(i,j,m)] += 2*dudx_left(yint[j])/h
    end
    if i==m+2  # right
        b[id(i,j,m)] += 2*dudx_right(yint[j])/h
    end

    if j==1 # bottom
        uij = u_bottom(xint[i])
        b[id(i,j,m)] += uij/h^2
    end
    if j==m # top
        uij = u_top(xint[i])
        b[id(i,j,m)] += uij/h^2
    end
end

# solve Au = f(xi,yj)
u = A\vec(b)
contourf(xint,yint,u) #vec(reshape(u,m,m)'))


text:
# # Getting started
Say you have a `NamedTuple` and you want to update it:
```jldoctest getting_started
julia> x = (greeting="Hello", name="World")
(greeting = "Hello", name = "World")

julia> x.greeting = "Hi"
ERROR: setfield! immutable struct of type NamedTuple cannot be changed
[...]
```
This fails, because named tuples are immutable. Instead you can use Accessors to carry out the update:

```jldoctest getting_started
julia> using Accessors

julia> @set x.greeting = "Hi"
(greeting = "Hi", name = "World")

julia> x # still the same. Accessors did not overwrite x, it just created an updated copy
(greeting = "Hello", name = "World")

julia> x_new = @set x.greeting = "Hi" # typically you will assign a name to the updated copy
(greeting = "Hi", name = "World")
```
Accessors.jl does not only support `NamedTuple`, but arbitrary structs and nested updates.
```jldoctest getting_started
julia> struct HelloWorld
           greeting::String
           name::String
       end

julia> x = HelloWorld("hi", "World")
HelloWorld("hi", "World")

julia> @set x.name = "Accessors" # update a struct
HelloWorld("hi", "Accessors")

julia> x = (a=1, b=(c=3, d=4))
(a = 1, b = (c = 3, d = 4))

julia> @set x.b.c = 10 # nested update
(a = 1, b = (c = 10, d = 4))
```
Accessors.jl does not only support updates of properties, but also index updates.
```jldoctest getting_started
julia> x = (10,20,21)
(10, 20, 21)

julia> @set x[3] = 30
(10, 20, 30)
```
In fact Accessors.jl supports many more notions of update:
```jldoctest getting_started
julia> x = [1,2,3];

julia> x_new = @set eltype(x) = UInt8;

julia> @show x_new;
x_new = UInt8[0x01, 0x02, 0x03]
```
Accessors.jl is very composable, which means different updates can be nested and combined.
```jldoctest getting_started
julia> data = (a = (b = (1,2),), c=3)
(a = (b = (1, 2),), c = 3)

julia> @set data.a.b[end] = 20
(a = (b = (1, 20),), c = 3)

julia> @set splitext("some_file.py")[2] = ".jl"
"some_file.jl"
```


text:
using Test
using MultivariateBases

@testset "Orthogonal" begin
    orthogonal_test(ChebyshevBasis, x -> [
        1,
        x,
        2x^2 - 1,
        4x^3 - 3x,
        8x^4 - 8x^2 + 1
    ], true)
    orthogonal_test(ChebyshevBasisSecondKind, x -> [
        1,
        2x,
        4x^2 - 1,
        8x^3 - 4x,
        16x^4 - 12x^2 + 1
    ], true)
end

@testset "API degree = $degree" for degree in 0:3
    api_test(ChebyshevBasis, degree)
    api_test(ChebyshevBasisSecondKind, degree)
end

max_stars_repo_name: JuliaAlgebra/PolynomialBases

max_issues_repo_name: JuliaAlgebra/PolynomialBases

max_forks_repo_name: JuliaAlgebra/PolynomialBases


text:
"""
Gives a reinterpreted view (of element type T) of the underlying array (of element type S).
If the size of `T` differs from the size of `S`, the array will be compressed/expanded in
the first dimension.
"""
struct ReinterpretArray{T,N,S,A<:AbstractArray{S, N}} <: AbstractArray{T, N}
    parent::A
    global reinterpret
    function reinterpret(::Type{T}, a::A) where {T,N,S,A<:AbstractArray{S, N}}
        function throwbits(::Type{S}, ::Type{T}, ::Type{U}) where {S,T,U}
            @_noinline_meta
            throw(ArgumentError("cannot reinterpret `$(S)` `$(T)`, type `$(U)` is not a bits type"))
        end
        function throwsize0(::Type{S}, ::Type{T})
            @_noinline_meta
            throw(ArgumentError("cannot reinterpret a zero-dimensional `$(S)` array to `$(T)` which is of a different size"))
        end
        function thrownonint(::Type{S}, ::Type{T}, dim)
            @_noinline_meta
            throw(ArgumentError("""
                cannot reinterpret an `$(S)` array to `$(T)` whose first dimension has size `$(dim)`.
                The resulting array would have non-integral first dimension.
            """))
        end
        isbits(T) || throwbits(S, T, T)
        isbits(S) || throwbits(S, T, S)
        (N != 0 || sizeof(T) == sizeof(S)) || throwsize0(S, T)
        if N != 0 && sizeof(S) != sizeof(T)
            dim = size(a)[1]
            rem(dim*sizeof(S),sizeof(T)) == 0 || thrownonint(S, T, dim)
        end
        new{T, N, S, A}(a)
    end
end

parent(a::ReinterpretArray) = a.parent

function size(a::ReinterpretArray{T,N,S} where {N}) where {T,S}
    psize = size(a.parent)
    size1 = div(psize[1]*sizeof(S), sizeof(T))
    tuple(size1, tail(psize)...)
end

unsafe_convert(::Type{Ptr{T}}, a::ReinterpretArray{T,N,S} where N) where {T,S} = Ptr{T}(unsafe_convert(Ptr{S},a.parent))

@inline @propagate_inbounds getindex(a::ReinterpretArray{T,0}) where {T} = reinterpret(T, a.parent[])
@inline @propagate_inbounds getindex(a::ReinterpretArray) = a[1]

@inline @propagate_inbounds function getindex(a::ReinterpretArray{T,N,S}, inds::Vararg{Int, N}) where {T,N,S}
    # Make sure to match the scalar reinterpret if that is applicable
    if sizeof(T) == sizeof(S) && (fieldcount(T) + fieldcount(S)) == 0
        return reinterpret(T, a.parent[inds...])
    else
        ind_start, sidx = divrem((inds[1]-1)*sizeof(T), sizeof(S))
        t = Ref{T}()
        s = Ref{S}()
        GC.@preserve t s begin
            tptr = Ptr{UInt8}(unsafe_convert(Ref{T}, t))
            sptr = Ptr{UInt8}(unsafe_convert(Ref{S}, s))
            i = 1
            nbytes_copied = 0
            # This is a bit complicated to deal with partial elements
            # at both the start and the end. LLVM will fold as appropriate,
            # once it knows the data layout
            while nbytes_copied < sizeof(T)
                s[] = a.parent[ind_start + i, tail(inds)...]
                while nbytes_copied < sizeof(T) && sidx < sizeof(S)
                    unsafe_store!(tptr, unsafe_load(sptr, sidx + 1), nbytes_copied + 1)
                    sidx += 1
                    nbytes_copied += 1
                end
                sidx = 0
                i += 1
            end
        end
        return t[]
    end
end

@inline @propagate_inbounds setindex!(a::ReinterpretArray{T,0,S} where T, v) where {S} = (a.parent[] = reinterpret(S, v))
@inline @propagate_inbounds setindex!(a::ReinterpretArray, v) = (a[1] = v)

@inline @propagate_inbounds function setindex!(a::ReinterpretArray{T,N,S}, v, inds::Vararg{Int, N}) where {T,N,S}
    v = convert(T, v)::T
    # Make sure to match the scalar reinterpret if that is applicable
    if sizeof(T) == sizeof(S) && (fieldcount(T) + fieldcount(S)) == 0
        return setindex!(a.parent, reinterpret(S, v), inds...)
    else
        ind_start, sidx = divrem((inds[1]-1)*sizeof(T), sizeof(S))
        t = Ref{T}(v)
        s = Ref{S}()
        GC.@preserve t s begin
            tptr = Ptr{UInt8}(unsafe_convert(Ref{T}, t))
            sptr = Ptr{UInt8}(unsafe_convert(Ref{S}, s))
            nbytes_copied = 0
            i = 1
            # Deal with any partial elements at the start. We'll have to copy in the
            # element from the original array and overwrite the relevant parts
            if sidx != 0
                s[] = a.parent[ind_start + i, tail(inds)...]
                while nbytes_copied < sizeof(T) && sidx < sizeof(S)
                    unsafe_store!(sptr, unsafe_load(tptr, nbytes_copied + 1), sidx + 1)
                    sidx += 1
                    nbytes_copied += 1
                end
                a.parent[ind_start + i, tail(inds)...] = s[]
                i += 1
                sidx = 0
            end
            # Deal with the main body of elements
            while nbytes_copied < sizeof(T) && (sizeof(T) - nbytes_copied) > sizeof(S)
                while nbytes_copied < sizeof(T) && sidx < sizeof(S)
                    unsafe_store!(sptr, unsafe_load(tptr, nbytes_copied + 1), sidx + 1)
                    sidx += 1
                    nbytes_copied += 1
                end
                a.parent[ind_start + i, tail(inds)...] = s[]
                i += 1
                sidx = 0
            end
            # Deal with trailing partial elements
            if nbytes_copied < sizeof(T)
                s[] = a.parent[ind_start + i, tail(inds)...]
                while nbytes_copied < sizeof(T) && sidx < sizeof(S)
                    unsafe_store!(sptr, unsafe_load(tptr, nbytes_copied + 1), sidx + 1)
                    sidx += 1
                    nbytes_copied += 1
                end
                a.parent[ind_start + i, tail(inds)...] = s[]
            end
        end
    end
    return a
end

max_stars_repo_name: nreigl/julia

max_issues_repo_name: BenChung/Julia-BTI

max_forks_repo_name: BenChung/Julia-BTI


text:
include("../core/eq-core.jl")

function test_compare()
    failing = 0
    if ! expect(Card(Ranks.ACE, Suits.CLUBS), Card(Ranks.KING, Suits.CLUBS), 1, 1, desc="ERROR: higher rank of same suit beats lower rank")
        failing += 1
    end
    if ! expect(Card(Ranks.ACE, Suits.HEARTS), Card(Ranks.KING, Suits.CLUBS), 1, 1, desc="ERROR: higher rank of better suit beats lower rank")
        failing += 1
    end
    if ! expect(Card(Ranks.ACE, Suits.CLUBS), Card(Ranks.KING, Suits.HEARTS), 1, 1, desc="ERROR: higher rank of worse suit beats lower rank")
        failing += 1
    end
    if ! expect(Card(Ranks.ACE, Suits.CLUBS), Card(Ranks.ACE, Suits.CLUBS), 0, 0, desc="ERROR: same rank of same suit is equal")
        failing += 1
    end
    if ! expect(Card(Ranks.ACE, Suits.HEARTS), Card(Ranks.ACE, Suits.CLUBS), 0, 1, desc="ERROR: same rank of better suit is equal (tie_break=false) or beats worse suit")
        failing += 1
    end
    if ! expect(Card(Ranks.ACE, Suits.CLUBS), Card(Ranks.KING, Suits.CLUBS), 1, 1, desc="ERROR: ACE beats KING")
        failing += 1
    end
    if ! expect(Card(Ranks.KING, Suits.CLUBS), Card(Ranks.QUEEN, Suits.CLUBS), 1, 1, desc="ERROR: KING beats QUEEN")
        failing += 1
    end
    if ! expect(Card(Ranks.QUEEN, Suits.CLUBS), Card(Ranks.JACK, Suits.CLUBS), 1, 1, desc="ERROR: QUEEN beats JACK")
        failing += 1
    end
    if ! expect(Card(Ranks.JACK, Suits.CLUBS), Card(Ranks.TEN, Suits.CLUBS), 1, 1, desc="ERROR: JACK beats TEN")
        failing += 1
    end
    if ! expect(Card(Ranks.TEN, Suits.CLUBS), Card(Ranks.NINE, Suits.CLUBS), 1, 1, desc="ERROR: TEN beats NINE")
        failing += 1
    end
    if ! expect(Card(Ranks.NINE, Suits.CLUBS), Card(Ranks.EIGHT, Suits.CLUBS), 1, 1, desc="ERROR: NINE beats EIGHT")
        failing += 1
    end
    if ! expect(Card(Ranks.EIGHT, Suits.CLUBS), Card(Ranks.SEVEN, Suits.CLUBS), 1, 1, desc="ERROR: EIGHT beats SEVEN")
        failing += 1
    end
    if ! expect(Card(Ranks.SEVEN, Suits.CLUBS), Card(Ranks.SIX, Suits.CLUBS), 1, 1, desc="ERROR: SEVEN beats SIX")
        failing += 1
    end
    if ! expect(Card(Ranks.SIX, Suits.CLUBS), Card(Ranks.FIVE, Suits.CLUBS), 1, 1, desc="ERROR: SIX beats FIVE")
        failing += 1
    end
    if ! expect(Card(Ranks.FIVE, Suits.CLUBS), Card(Ranks.FOUR, Suits.CLUBS), 1, 1, desc="ERROR: FIVE beats FOUR")
        failing += 1
    end
    if ! expect(Card(Ranks.FOUR, Suits.CLUBS), Card(Ranks.THREE, Suits.CLUBS), 1, 1, desc="ERROR: FOUR beats THREE")
        failing += 1
    end
    if ! expect(Card(Ranks.THREE, Suits.CLUBS), Card(Ranks.TWO, Suits.CLUBS), 1, 1, desc="THREE beats TWO")
        failing += 1
    end
    return failing
end

function test_compare_hand()
    failing = 0
    if ! expect(Rankings.ROYAL, Rankings.SFLUSH, 1, 1, desc="ERROR: ROYAL beats SFLUSH")
        failing += 1
    end
    if ! expect(Rankings.SFLUSH, Rankings.QUAD, 1, 1, desc="ERROR: SFLUSH beats QUAD")
        failing += 1
    end
    if ! expect(Rankings.QUAD, Rankings.FULL, 1, 1, desc="ERROR: QUAD beats FULL")
        failing += 1
    end
    if ! expect(Rankings.FULL, Rankings.FLUSH, 1, 1, desc="ERROR: FULL beats FLUSH")
        failing += 1
    end
    if ! expect(Rankings.FLUSH, Rankings.STRAIGHT, 1, 1, desc="ERROR: FLUSH beats STRAIGHT")
        failing += 1
    end
    if ! expect(Rankings.STRAIGHT, Rankings.TRIP, 1, 1, desc="ERROR: STRAIGHT beats TRIP")
        failing += 1
    end
    if ! expect(Rankings.TRIP, Rankings.DUBS, 1, 1, desc="ERROR: TRIP beats DUBS")
        failing += 1
    end
    if ! expect(Rankings.DUBS, Rankings.PAIR, 1, 1, desc="ERROR: DUBS beats PAIR")
        failing += 1
    end
    if ! expect(Rankings.PAIR, Rankings.HIGH, 1, 1, desc="ERROR: PAIR beats HIGH")
        failing += 1
    end
    return failing
end

function expect(c1::Card, c2::Card, v1::Int, v2::Int; desc="")
    if ! (compare(c1, c2, tie_break=false) == v1 && compare(c1, c2, tie_break=true) == v2)
        println(desc)
        return false
    end
    if ! (compare(c2, c1, tie_break=false) == v1*-1 && compare(c2, c1, tie_break=true) == v2*-1)
        println(desc)
        return false
    end
    return true
end

function expect(r1::Rankings.Ranking, r2::Rankings.Ranking, v1::Int, v2::Int; desc="")
    if ! (compare(ranked_hand(r1), ranked_hand(r2), tie_break=false) == v1 && compare(ranked_hand(r1), ranked_hand(r2), tie_break=true) == v2)
        println(desc)
        return false
    end
    if ! (compare(ranked_hand(r2), ranked_hand(r1), tie_break=false) == v1*-1 && compare(ranked_hand(r2), ranked_hand(r1), tie_break=true) == v2*-1)
        println(desc)
        return false
    end
    return true
end

function ranked_hand(r::Rankings.Ranking)
    return Hand(0, r, null_card(), null_card(), null_card(), null_card(), null_card())
end

function null_card()
    return Card(Ranks.NULL, Suits.NULL)
end

function run_tests()
    failing = 0
    failing += test_compare()
    failing += test_compare_hand()
    println("Done!")
    if failing == 0
        println("All tests passing!")
    else
        println("$(failing) tests failing")
    end
end

run_tests()

max_stars_repo_name: ttowncompiled/EQPoker

max_issues_repo_name: ttowncompiled/EQPoker

max_forks_repo_name: ttowncompiled/EQPoker


text:
# forward pass of the SDDP algorithm
function noDisruptionBuild(Δt, T, qpopt = false, solveOpt = true)
    # precalculate data
    Rdict = Dict();
    Xdict = Dict();
    for k in fData.brList
        Rdict[k] = fData.g[k]/(fData.g[k]^2 + fData.b[k]^2);
        Xdict[k] = -fData.b[k]/(fData.g[k]^2 + fData.b[k]^2);
    end
    Ω = [ω for ω in keys(pDistr.ωDistrn)];

    # construct the first stage without disruption occurring
    if qpopt
        mp = Model(solver = IpoptSolver(print_level = 0, linear_solver = "ma27"));
    else
        mp = Model(solver = GurobiSolver(GUROBI_ENV,OutputFlag = 0,Threads = 1));
    end

    # set up the variables
    @variable(mp, fData.Pmin[i] <= sp[i in fData.genIDList, t in 1:T] <= fData.Pmax[i]);
    @variable(mp, fData.Qmin[i] <= sq[i in fData.genIDList, t in 1:T] <= fData.Qmax[i]);
    sphatsum = Dict();
    for t in 1:T
        for i in fData.IDList
            sphatsum[i,t] = @expression(mp,0.0);
            if i in keys(fData.LocRev)
                for j in fData.LocRev[i]
                    sphatsum[i,t] += sp[j,t];
                end
            end
        end
    end
    sqhatsum = Dict();
    for t in 1:T
        for i in fData.IDList
            sqhatsum[i,t] = @expression(mp,0.0);
            if i in keys(fData.LocRev)
                for j in fData.LocRev[i]
                    sqhatsum[i,t] += sq[j,t];
                end
            end
        end
    end

    @variable(mp, p[k in fData.brList, t in 1:T]);
    @variable(mp, q[k in fData.brList, t in 1:T]);
    @variable(mp, fData.Vmin[i]^2 <= v[i in fData.IDList, t in 1:T] <= fData.Vmax[i]^2);
    @variable(mp, 0 <= w[i in bData.IDList, t in 0:T] <= bData.cap[i]);
    @variable(mp, y[i in bData.IDList, t in 1:T]);
    @variable(mp, zp[i in bData.IDList, t in 1:T]);
    @variable(mp, zq[i in bData.IDList, t in 1:T]);
    @variable(mp, lpp[i in fData.IDList, t in 1:T] >= 0);
    @variable(mp, lqp[i in fData.IDList, t in 1:T] >= 0);
    @variable(mp, lpm[i in fData.IDList, t in 1:T] >= 0);
    @variable(mp, lqm[i in fData.IDList, t in 1:T] >= 0);
    @variable(mp, 0 <= u[i in bData.IDList] <= bData.uCap[i]);
    @variable(mp, θ[tp in 2:T, ω in Ω] >= 0);

    # set up the constraints
    @constraint(mp, pBalance[i in fData.IDList, t in 1:T], sum(zp[b,t] for b in bData.IDList if bData.Loc[b] == i) + lpp[i,t] - lpm[i,t] +
        sphatsum[i,t] - dData.pd[i][t] == sum(p[k,t] for k in fData.branchDict1[i]));
    @constraint(mp, qBalance[i in fData.IDList, t in 1:T], sum(zq[b,t] for b in bData.IDList if bData.Loc[b] == i) + lqp[i,t] - lqm[i,t] +
        sqhatsum[i,t] - dData.qd[i][t] == sum(q[k,t] for k in fData.branchDict1[i]));
    @constraint(mp, pequal[k in fData.brList, t in 1:T], p[k,t] == -p[(k[2],k[1],k[3]),t]);
    @constraint(mp, qequal[k in fData.brList, t in 1:T], q[k,t] == -q[(k[2],k[1],k[3]),t]);
    if qpopt
        @constraint(mp, lineThermal[k in fData.brList, t in 1:T], p[k,t]^2 + q[k,t]^2 <= fData.rateA[k]^2);
        @constraint(mp, bThermal[i in bData.IDList, t in 1:T], zp[i,t]^2+zq[i,t]^2 <= u[i]^2);
    else
        @constraint(mp, lineThermal[k in fData.brList, t in 1:T], norm([p[k,t],q[k,t]]) <= fData.rateA[k]);
        @constraint(mp, bThermal[i in bData.IDList, t in 1:T], norm([zp[i,t],zq[i,t]]) <= u[i]);
    end
    @constraint(mp, powerflow[k in fData.brList, t in 1:T], v[k[2],t] == v[k[1],t] - 2*(Rdict[k]*p[k,t] + Xdict[k]*q[k,t]));
    @constraint(mp, rampUp[i in fData.genIDList, t in 2:T], sp[i,t] - sp[i,t - 1] <= fData.RU[i]);
    @constraint(mp, rampDown[i in fData.genIDList, t in 2:T], sp[i,t] - sp[i,t - 1] >= fData.RD[i]);
    @constraint(mp, bInv[i in bData.IDList, t in 1:T], w[i,t] == w[i,t-1] - y[i,t]*Δt);
    @constraint(mp, bEfficient[i in bData.IDList, l in 1:length(bData.ηα[i]), t in 1:T], zp[i,t] <= bData.ηα[i][l]*y[i,t] + bData.ηβ[i][l]);
    @constraint(mp, bInvIni[i in bData.IDList], w[i,0] == bData.bInv[i]);

    # set up the cuts, here tp is the disruption time
    for tp in 2:T
        for ω in Ω
            if (tp,ω) in keys(cutDict)
                for l in 1:length(cutDict[tp,ω])
                    @constraint(mp, θ[tp,ω] >= cutDict[tp,ω][l].vhat +
                        sum(cutDict[tp,ω][l].λ[i]*sp[i,tp - 1] for i in fData.genIDList) +
                        sum(cutDict[tp,ω][l].γ[i]*w[i,tp - 1] + cutDict[tp,ω][l].μ[i]*u[i] for i in bData.IDList));
                end
            end
        end
    end

    # set up the objective function
    @variable(mp,fs[i in fData.genIDList, t in 1:T]);
    if qpopt
        @constraint(mp, genCost[i in fData.genIDList, t in 1:T], fs[i,t] ==
            sum(fData.cp[i].params[nco]*sp[i,t]^(fData.cp[i].n - nco) for nco in 1:(fData.cp[i].n - 1)));
    else
        @variable(mp,tAux1[i in fData.genIDList, t in 1:T; fData.cp[i].n == 3] >= 0);
        @variable(mp,tAux2[i in fData.genIDList, t in 1:T; fData.cp[i].n == 3]);
        @variable(mp,tAux3[i in fData.genIDList, t in 1:T; fData.cp[i].n == 3]);

        @constraint(mp,gcAux1[i in fData.genIDList, t in 1:T; fData.cp[i].n == 3], tAux1[i,t] == fs[i,t] +
            (fData.cp[i].params[1] + fData.cp[i].params[2]^2)/(4*fData.cp[i].params[1]));
        @constraint(mp,gcAux2[i in fData.genIDList, t in 1:T; fData.cp[i].n == 3], tAux2[i,t] == fs[i,t] +
            (-fData.cp[i].params[1] + fData.cp[i].params[2]^2)/(4*fData.cp[i].params[1]));
        @constraint(mp,gcAux3[i in fData.genIDList, t in 1:T; fData.cp[i].n == 3], tAux3[i,t] == sqrt(fData.cp[i].params[1])*sp[i,t] +
            fData.cp[i].params[2]/(2*sqrt(fData.cp[i].params[1])));
        @constraint(mp, genCost1[i in fData.genIDList, t in 1:T; fData.cp[i].n == 3],
            norm([tAux2[i,t],tAux3[i,t]]) <= tAux1[i,t]);
        @constraint(mp, genCost2[i in fData.genIDList, t in 1:T; fData.cp[i].n == 2],
            fs[i,t] == fData.cp[i].params[1]*sp[i,t]);
    end

    fDict,lDict,θDict = calDualC1(T, fData, pDistr);
    objExpr = @expression(mp, sum(bData.cost[i]*u[i] for i in bData.IDList) +
        sum(sum(fDict[i,t]*fs[i,t] for i in fData.genIDList) +
        sum(lDict[i,t]*(lpp[i,t] + lpm[i,t] + lqp[i,t] + lqm[i,t]) for i in fData.IDList) for t in 1:T) +
        sum(sum(θDict[t,ω]*θ[t,ω] for ω in Ω) for t in 2:T));

    @objective(mp, Min, objExpr);

    if solveOpt
        # solve the problem
        # optimize!(mp, with_optimizer(Gurobi.Optimizer, GUROBI_ENV, OutputFlag = 0,
        #     QCPDual = 1, NumericFocus = 3, BarQCPConvTol = 1e-9, FeasibilityTol = 1e-9));
        statusMp = solve(mp);
        if statusMp != :Optimal
            mp.solver = GurobiSolver(GUROBI_ENV,OutputFlag = 0,NumericFocus = 3,Threads = 1);
            statusMp = solve(mp);
        end
        #with_optimizer(Ipopt.Optimizer, linear_solver = "ma27", acceptable_tol = 1e-8, print_level = 0, max_iter = 10000));
        mpObj = getobjectivevalue(mp);
        println("First stage, solving status $(statusMp)");
        # obtain the solutions
        solSp = Dict();
        solSq = Dict();
        solw = Dict();
        solu = Dict();
        solLp = Dict();
        solLq = Dict();
        solzp = Dict();
        for i in fData.genIDList
            for t in 1:T
                if abs(getvalue(sp[i,t])) > 1e-5
                    solSp[i,t] = getvalue(sp[i,t]);
                else
                    solSp[i,t] = 0;
                end
                if abs(getvalue(sq[i,t])) > 1e-5
                    solSq[i,t] = getvalue(sq[i,t]);
                else
                    solSq[i,t] = 0;
                end
            end
        end
        for i in bData.IDList
            if abs(getvalue(u[i])) > 1e-5
                solu[i] = getvalue(u[i]);
                for t in 1:T
                    if abs(getvalue(w[i,t])) > 1e-5
                        solw[i,t] = getvalue(w[i,t]);
                    else
                        solw[i,t] = 0;
                    end
                    if abs(getvalue(zp[i,t])) > 1e-6
                        solzp[i,t] = getvalue(zp[i,t]);
                    else
                        solzp[i,t] = 0;
                    end
                end
            else
                solu[i] = 0;
                for t in 1:T
                    solw[i,t] = 0;
                    if abs(getvalue(zp[i,t])) > 1e-6
                        solzp[i,t] = getvalue(zp[i,t]);
                    else
                        solzp[i,t] = 0;
                    end
                end
            end
        end
        for i in fData.IDList
            for t in 1:T
                if (abs(getvalue(lpp[i,t])) > 1e-8)|(abs(getvalue(lpm[i,t])) > 1e-8)
                    solLp[i,t] = getvalue(lpp[i,t]) - getvalue(lpm[i,t]);
                else
                    solLp[i,t] = 0;
                end
                if (abs(getvalue(lqp[i,t])) > 1e-8)|(abs(getvalue(lqm[i,t])) > 1e-8)
                    solLq[i,t] = getvalue(lqp[i,t]) - getvalue(lqm[i,t]);
                else
                    solLq[i,t] = 0;
                end
            end
        end

        sol = solData(solSp,solSq,solw,solu,solLp,solLq,solzp);
        return sol,mpObj;
    else
        return mp;
    end
end

function fBuild(td, ωd, currentSol, τ, Δt, T, qpopt = false, solveOpt = true, hardened = [])
    # precalculate data
    Rdict = Dict();
    Xdict = Dict();
    for k in fData.brList
        Rdict[k] = fData.g[k]/(fData.g[k]^2 + fData.b[k]^2);
        Xdict[k] = -fData.b[k]/(fData.g[k]^2 + fData.b[k]^2);
    end
    Ω = [ω for ω in keys(pDistr.ωDistrn)];
    Bparams = Dict();
    for t in td:T
        # create B parameters
        for k in fData.brList
            # if the line is disrupted and it is within disruption time
            if (((k[1],k[2]) in ωd)|((k[2],k[1]) in ωd))&(t <= td + τ)
                if ((k[1],k[2]) in hardened) | ((k[2],k[1]) in hardened)
                    Bparams[k,t] = 1;
                else
                    Bparams[k,t] = 0;
                end
            else
                Bparams[k,t] = 1;
            end
        end
        for i in fData.genIDList
            if (i in ωd)&(t <= td + τ)
                if i in hardened
                    Bparams[i,t] = 1;
                else
                    Bparams[i,t] = 0;
                end
            else
                Bparams[i,t] = 1;
            end
        end
    end
    for i in fData.genIDList
        Bparams[i,td - 1] = 1;
    end

    if qpopt
        mp = Model(solver = IpoptSolver(print_level = 0, linear_solver = "ma27"));
    else
        mp = Model(solver = GurobiSolver(GUROBI_ENV,OutputFlag = 0,Threads = 1));
    end

    # set up the variables
    @variable(mp, fData.Pmin[i]*Bparams[i,t] <= sp[i in fData.genIDList,t in (td - 1):T] <= fData.Pmax[i]*Bparams[i,t]);
    @variable(mp, fData.Qmin[i]*Bparams[i,t] <= sq[i in fData.genIDList,t in td:T] <= fData.Qmax[i]*Bparams[i,t]);
    sphatsum = Dict();
    for t in td:T
        for i in fData.IDList
            sphatsum[i,t] = @expression(mp,0.0);
            if i in keys(fData.LocRev)
                for j in fData.LocRev[i]
                    sphatsum[i,t] += sp[j,t];
                end
            end
        end
    end
    sqhatsum = Dict();
    for t in td:T
        for i in fData.IDList
            sqhatsum[i,t] = @expression(mp,0.0);
            if i in keys(fData.LocRev)
                for j in fData.LocRev[i]
                    sqhatsum[i,t] += sq[j,t];
                end
            end
        end
    end

    @variable(mp, p[k in fData.brList, t in td:T]);
    @variable(mp, q[k in fData.brList, t in td:T]);
    @variable(mp, fData.Vmin[i]^2 <= v[i in fData.IDList, t in td:T] <= fData.Vmax[i]^2);
    @variable(mp, 0 <= w[i in bData.IDList, t in (td - 1):T] <= bData.cap[i]);
    @variable(mp, y[i in bData.IDList, t in td:T]);
    @variable(mp, zp[i in bData.IDList, t in td:T]);
    @variable(mp, zq[i in bData.IDList, t in td:T]);
    @variable(mp, lpp[i in fData.IDList, t in td:T] >= 0);
    @variable(mp, lqp[i in fData.IDList, t in td:T] >= 0);
    @variable(mp, lpm[i in fData.IDList, t in td:T] >= 0);
    @variable(mp, lqm[i in fData.IDList, t in td:T] >= 0);
    @variable(mp, 0 <= u[i in bData.IDList] <= bData.uCap[i]);
    @variable(mp, θ[tp in (td + τ + 1):T, ω in Ω] >= 0);

    # set up the constraints
    @constraint(mp, pBalance[i in fData.IDList, t in td:T], sum(zp[b,t] for b in bData.IDList if bData.Loc[b] == i) + lpp[i,t] - lpm[i,t] +
        sphatsum[i,t] - dData.pd[i][t] == sum(p[k,t] for k in fData.branchDict1[i]));
    @constraint(mp, qBalance[i in fData.IDList, t in td:T], sum(zq[b,t] for b in bData.IDList if bData.Loc[b] == i) + lqp[i,t] - lqm[i,t] +
        sqhatsum[i,t] - dData.qd[i][t] == sum(q[k,t] for k in fData.branchDict1[i]));
    @constraint(mp, pequal[k in fData.brList, t in td:T], p[k,t] == -p[(k[2],k[1],k[3]),t]);
    @constraint(mp, qequal[k in fData.brList, t in td:T], q[k,t] == -q[(k[2],k[1],k[3]),t]);
    if qpopt
        @constraint(mp, lineThermal1[k in fData.brList, t in td:T;Bparams[k,t] == 1], p[k,t]^2+q[k,t]^2 <= fData.rateA[k]^2);
        @constraint(mp, bThermal[i in bData.IDList, t in td:T], zp[i,t]^2+zq[i,t]^2 <= u[i]^2);
    else
        @constraint(mp, lineThermal1[k in fData.brList, t in td:T;Bparams[k,t] == 1], norm([p[k,t],q[k,t]]) <= fData.rateA[k]);
        @constraint(mp, bThermal[i in bData.IDList, t in td:T], norm([zp[i,t],zq[i,t]]) <= u[i]);
    end
    @constraint(mp, lineThermal2[k in fData.brList, t in td:T;Bparams[k,t] == 0], p[k,t] == 0);
    @constraint(mp, lineThermal3[k in fData.brList, t in td:T;Bparams[k,t] == 0], q[k,t] == 0);
    @constraint(mp, powerflow1[k in fData.brList, t in td:T;Bparams[k,t] == 1], v[k[2],t] == v[k[1],t] - 2*(Rdict[k]*p[k,t] + Xdict[k]*q[k,t]));
    @constraint(mp, rampUp[i in fData.genIDList, t in td:T; Bparams[i,t] == 1], sp[i,t] - sp[i,t - 1] <= fData.RU[i]);
    @constraint(mp, rampDown[i in fData.genIDList, t in td:T; Bparams[i,t] == 1], sp[i,t] - sp[i,t - 1] >= fData.RD[i]);
    @constraint(mp, bInv[i in bData.IDList, t in td:T], w[i,t] == w[i,t-1] - y[i,t]*Δt);
    # @constraint(mp, bThermal1[i in bData.IDList, t in td:T], zp[i,t] <= u[i]);
    # @constraint(mp, bThermal3[i in bData.IDList, t in td:T], -zp[i,t] <= u[i]);
    # #@constraint(mp, bThermal2[i in bData.IDList, t in 1:T], zq[i,t] <= u[i]);
    # @constraint(mp, bThermal2[i in bData.IDList, t in td:T], zq[i,t] == 0);
    @constraint(mp, bEfficient[i in bData.IDList, l in 1:length(bData.ηα[i]), t in td:T], zp[i,t] <= bData.ηα[i][l]*y[i,t] + bData.ηβ[i][l]);
    @constraint(mp, bInvIni[i in bData.IDList], w[i,td - 1] == currentSol.w[i,td - 1]);
    @constraint(mp, spIni[i in fData.genIDList], sp[i,td - 1] == currentSol.sp[i,td - 1]);
    @constraint(mp, uIni[i in bData.IDList], u[i] == currentSol.u[i]);

    # set up the cuts, here tp is the disruption time
    for tp in (td + τ + 1):T
        for ω in Ω
            if (tp,ω) in keys(cutDict)
                for l in 1:length(cutDict[tp,ω])
                    @constraint(mp, θ[tp,ω] >= cutDict[tp,ω][l].vhat +
                        sum(cutDict[tp,ω][l].λ[i]*sp[i,tp - 1] for i in fData.genIDList) +
                        sum(cutDict[tp,ω][l].γ[i]*w[i,tp - 1] + cutDict[tp,ω][l].μ[i]*u[i] for i in bData.IDList));
                end
            end
        end
    end

    # set up the objective function
    @variable(mp,fs[i in fData.genIDList, t in td:T]);
    if qpopt
        @constraint(mp, genCost[i in fData.genIDList, t in td:T], fs[i,t] ==
            sum(fData.cp[i].params[nco]*sp[i,t]^(fData.cp[i].n - nco) for nco in 1:(fData.cp[i].n - 1)));
    else
        @variable(mp,tAux1[i in fData.genIDList, t in td:T; fData.cp[i].n == 3] >= 0);
        @variable(mp,tAux2[i in fData.genIDList, t in td:T; fData.cp[i].n == 3]);
        @variable(mp,tAux3[i in fData.genIDList, t in td:T; fData.cp[i].n == 3]);

        @constraint(mp,gcAux1[i in fData.genIDList, t in td:T; fData.cp[i].n == 3], tAux1[i,t] == fs[i,t] +
            (fData.cp[i].params[1] + fData.cp[i].params[2]^2)/(4*fData.cp[i].params[1]));
        @constraint(mp,gcAux2[i in fData.genIDList, t in td:T; fData.cp[i].n == 3], tAux2[i,t] == fs[i,t] +
            (-fData.cp[i].params[1] + fData.cp[i].params[2]^2)/(4*fData.cp[i].params[1]));
        @constraint(mp,gcAux3[i in fData.genIDList, t in td:T; fData.cp[i].n == 3], tAux3[i,t] == sqrt(fData.cp[i].params[1])*sp[i,t] +
            fData.cp[i].params[2]/(2*sqrt(fData.cp[i].params[1])));
        @constraint(mp, genCost1[i in fData.genIDList, t in td:T; fData.cp[i].n == 3],
            norm([tAux2[i,t],tAux3[i,t]]) <= tAux1[i,t]);
        @constraint(mp, genCost2[i in fData.genIDList, t in td:T; fData.cp[i].n == 2],
            fs[i,t] == fData.cp[i].params[1]*sp[i,t]);
    end

    fDict,lDict,θDict = calDualC(td, τ, T, fData, pDistr);
    objExpr = @expression(mp, sum(sum(fDict[i,t]*fs[i,t] for i in fData.genIDList) +
        sum(lDict[i,t]*(lpp[i,t] + lpm[i,t] + lqp[i,t] + lqm[i,t]) for i in fData.IDList) for t in td:T) +
        sum(sum(θDict[t,ω]*θ[t,ω] for ω in Ω) for t in (td + τ + 1):T));
    @objective(mp, Min, objExpr);

    if solveOpt
        # solve the problem
        # optimize!(mp, with_optimizer(Gurobi.Optimizer, GUROBI_ENV, OutputFlag = 0,
        #     QCPDual = 1, NumericFocus = 3, BarQCPConvTol = 1e-9, FeasibilityTol = 1e-9));
        statusMp = solve(mp);
        if statusMp != :Optimal
            mp.solver = GurobiSolver(GUROBI_ENV,OutputFlag = 0,NumericFocus = 3,Threads = 1);
            statusMp = solve(mp);
        end
        #optimize!(mp, with_optimizer(Ipopt.Optimizer, linear_solver = "ma27", acceptable_tol = 1e-8, print_level = 0, max_iter = 10000));
        mpObj = getobjectivevalue(mp);
        println("Disruption time $(td), scenario $(ωd), solving status $(statusMp)");
        # obtain the solutions
        solSp = Dict();
        solSq = Dict();
        solw = Dict();
        solu = Dict();
        solLp = Dict();
        solLq = Dict();
        solzp = Dict();
        for i in fData.genIDList
            for t in td:T
                if abs(getvalue(sp[i,t])) > 1e-5
                    solSp[i,t] = getvalue(sp[i,t]);
                else
                    solSp[i,t] = 0;
                end
                if abs(getvalue(sq[i,t])) > 1e-5
                    solSq[i,t] = getvalue(sq[i,t]);
                else
                    solSq[i,t] = 0;
                end
            end
        end
        for i in bData.IDList
            if abs(getvalue(u[i])) > 1e-5
                solu[i] = getvalue(u[i]);
                for t in td:T
                    if abs(getvalue(w[i,t])) > 1e-5
                        solw[i,t] = getvalue(w[i,t]);
                    else
                        solw[i,t] = 0;
                    end
                    if abs(getvalue(zp[i,t])) > 1e-6
                        solzp[i,t] = getvalue(zp[i,t]);
                    else
                        solzp[i,t] = 0;
                    end
                end
            else
                solu[i] = 0;
                for t in td:T
                    solw[i,t] = 0;
                    if abs(getvalue(zp[i,t])) > 1e-6
                        solzp[i,t] = getvalue(zp[i,t]);
                    else
                        solzp[i,t] = 0;
                    end
                end
            end
        end
        for i in fData.IDList
            for t in td:T
                if (abs(getvalue(lpp[i,t])) > 1e-8)|(abs(getvalue(lpm[i,t])) > 1e-8)
                    solLp[i,t] = getvalue(lpp[i,t]) - getvalue(lpm[i,t]);
                else
                    solLp[i,t] = 0;
                end
                if (abs(getvalue(lqp[i,t])) > 1e-8)|(abs(getvalue(lqm[i,t]))> 1e-8)
                    solLq[i,t] = getvalue(lqp[i,t]) - getvalue(lqm[i,t]);
                else
                    solLq[i,t] = 0;
                end
            end
        end

        sol = solData(solSp,solSq,solw,solu,solLp,solLq,solzp);
        return sol,mpObj;
    else
        return mp;
    end
end

function constructForwardM(td, ωd, sol, τ, Δt, T, qpopt = false, hardenend = [])
    # construct the math program given the state variables and current stage
    if td == 1
        # if it is the no-disruption problem
        sol,objV = noDisruptionBuild(Δt, T, qpopt);
    else
        # if it is f_{ht}^ω
        sol,objV = fBuild(td, ωd, sol, τ, Δt, T, qpopt, true, hardenend);
    end
    return sol,objV;
end

function buildPath(T, Δt, qpopt = false, pathList = [], hardened = [])
    disT = 1;
    ωd = 0;
    τω = 0;
    costn = 0;
    solHist = [];
    currentLB = 0;
    currentSol = solData(Dict(),Dict(),Dict(),Dict(),Dict(),Dict(),Dict());
    iter = 1;
    while disT <= T
        # solve the current stage problem, state variables are passed
        nowT = disT;
        currentSol,objV = constructForwardM(disT, ωd, currentSol, τω, Δt, T, qpopt, hardened);
        push!(solHist,(currentSol,nowT,ωd,τω));

        # generate disruption
        if pathList == []
            tp,ωd,τω = genScenario(pDistr);
        else
            tp,ωd,τω = pathList[iter];
        end
        iter += 1;
        if nowT == 1
            currentLB = objV;
            disT += tp;
            disT = min(disT, T + 1);
            # calculate the cost of the solution until the next disruption time
            costn += sum(currentSol.u[i]*bData.cost[i] for i in bData.IDList);
            costn = calCostF(costn, currentSol, T, fData, nowT, disT);
        else
            disT += tp + τω;
            disT = min(disT, T + 1);
            # calculate the cost of the solution until the next disruption time
            costn = calCostF(costn, currentSol, T, fData, nowT, disT);
        end
    end
    println("Path Built!");
    return [solHist,currentLB,costn];
end

function exeForward(T, Δt, N, qpopt = false, pathDict = Dict(), hardened = [])
    # execution of forward pass
    # input: N: the number of trial points;
    #       cutDict: set of currently generated cuts (global in every core)
    # output: solList: a list of solution paths
    solDict = Dict();
    costDict = Dict();
    objV = 0;
    currentLB = 0;
    # for n in 1:N
    #     # for each trial path
    #     returnData = buildPath(τ, T, Δt, fData, bData, dData, pDistr, cutDict);
    #     solDict[n] = returnData[1];
    #     costDict[n] = returnData[2];
    # end
    if pathDict == Dict()
        for i in 1:N
            pathDict[i] = [];
        end
    end
    returnData = pmap(i -> buildPath(T, Δt, qpopt, pathDict[i], hardened), 1:N);
    for n in 1:N
        solDict[n] = returnData[n][1];
        costDict[n] = returnData[n][3];
    end
    currentLB = returnData[1][2];
    return solDict, currentLB, costDict;
end

function exeForward_simuOpt(T, Δt, N, iterNo, qpopt = false, hardened = [])
    # execution of forward pass, with samples selected to cover more time
    # input: N: the number of trial points;
    #       cutDict: set of currently generated cuts (global in every core)
    # output: solList: a list of solution paths
    solDict = Dict();
    costDict = Dict();
    objV = 0;
    currentLB = 0;
    genCutsCount = Dict();
    for t in 2:T
        genCutsCount[t] = 0;
    end
    pathDict,pathTimeList = pathSimu_cover(N, T, pDistr, genCutsCount, iterNo);
    for i in 1:N
        # update genCutsCount
        for j in 1:length(pathTimeList[i])
            genCutsCount[pathTimeList[i][j]] += 1;
        end
    end
    returnData = pmap(i -> buildPath(T, Δt, qpopt, pathDict[i], hardened), 1:N);
    for n in 1:N
        solDict[n] = returnData[n][1];
        costDict[n] = returnData[n][3];
    end
    currentLB = returnData[1][2];
    return solDict, currentLB, costDict;
end

function exeForward_last(T, Δt, N, qpopt = false, hardened = [])
    # execution of forward pass, with samples selected to cover more time
    # input: N: the number of trial points;
    #       cutDict: set of currently generated cuts (global in every core)
    # output: solList: a list of solution paths
    solDict = Dict();
    costDict = Dict();
    objV = 0;
    currentLB = 0;
    pathDict,pathTimeList = pathSimu_cover_last(N, T, pDistr);
    returnData = pmap(i -> buildPath(T, Δt, qpopt, pathDict[i], hardened), 1:N);
    for n in 1:N
        solDict[n] = returnData[n][1];
        costDict[n] = returnData[n][3];
    end
    currentLB = returnData[1][2];
    return solDict, currentLB, costDict;
end


text:
using BinDeps
@BinDeps.setup

# download https:////d3js.org/d3.v4.min.js to deps/ if it doesn't exist
local_fn = joinpath(dirname(@__FILE__), "d3.min.js")
if !isfile(local_fn)
	info("Cannot find deps/d3.min.js... downloading latest version.")
	download("https://d3js.org/d3.v4.min.js", local_fn)
end

@BinDeps.install


text:
############
# Analytic second derivative variable transform.

module HessianReparameterization

using ForwardDiff
using DualNumbers
using VariationalModelIndices

VERSION < v"0.4.0-dev" && using Docile

export transform_hessian, get_moment_hess
export get_dx_dy_func, get_d2x_dy2_funcs, get_d2x_dy2
export get_df_dy_func, get_d2f_dy2_func

@doc """
Use forward differentiation to get a jacobian from a variable transformation
function.

Args:
	- y_to_x: A function that takes a vector y to a vector x.
	          Note that it must map generic numeric vectors to numeric vectors.
	- K: The size of the vectors

Returns:
	- A function that gets the Jacobian of the transform evaluated at y.
""" ->
function get_dx_dy_func(y_to_x::Function, K::Int64)

	# Calculate the necessary derivatives:
	dx_dy_func = ForwardDiff.jacobian(y_to_x, mutates=false)

	# The old code transposed the Jecobian, which is misleading -- it was only
	# correct because I accidentally also transposed the Jacobian inverse in
	# the notes.
	# function dx_dy_func{T <: Number}(y::Array{T})
	# 	dx_dy_func_transpose(y)'
	# end

	return dx_dy_func
end

function get_d2x_dy2_funcs(y_to_x::Function, K::Int64)
	function y_to_x{T <: Number}(y::Array{T}, i::Int64)
		y_to_x(y)[i]
	end
	[ ForwardDiff.hessian(y -> y_to_x(y, i), mutates=false) for i=1:K ]
end

function get_df_dy_func(f::Function, K::Int64)
	ForwardDiff.gradient(f, mutates=false)
end

function get_d2f_dy2_func(f::Function, K::Int64)
	ForwardDiff.hessian(f, mutates=false)
end

function get_d2x_dy2(d2x_dy2_funcs::Array{Function}, y::Array{Float64, 1})
	K = length(y)
	result = Float64[ d2x_dy2_funcs[k](y)[i, j] for i=1:K, j=1:K, k=1:K ]
	result
end

function transform_hessian(dx_dy::Array{Float64, 2}, d2x_dy2::Array{Float64, 3},
	                         df_dy::Array{Float64, 1}, d2f_dy2::Array{Float64, 2})
	# Input variables:
	#   dx_dy[i, j] = dx[i] / dy[j]
	#   d2x_dy2[i, j, k] = d2 x[k] / (dy[i] dy[j])
	#   df_dy[i] = df / dy[i]
	#   d2f_dy2[i, j] = d2f / (dy[i] dy[j])
	#
	# Returns:
	#   d2f / dx2

	K = size(dx_dy, 1)
	@assert K == size(dx_dy, 2) == size(d2x_dy2, 1) == size(d2x_dy2, 2) ==
	        size(df_dy, 1) == size(d2f_dy2, 1) == size(d2f_dy2, 2)

	# The same as dx_dy' \ d2f_dy2 * inv(dx_dy).  This is the only term
	# if you are at an optimum since in that case df_dy = 0.
	opt_term = dx_dy' \ (dx_dy' \ d2f_dy2)'

	df_dx = dx_dy' \ df_dy

	# This is dJ'_dx[i, j, k] = dJ[i, j] / dx[k] = d / dx[k] (dx[j] / dy[i])
	dJt_dx = Float64[ (dx_dy' \ slice(d2x_dy2, :, i, j)[:])[k] for i=1:K, j=1:K, k=1:K ]

	j_term = Array(Float64, K, K)
	for j=1:K
		j_term[:, j] = -(dx_dy' \ slice(dJt_dx, :, :, j)[:, :]) * df_dx
	end

	j_term + opt_term
end

@doc """
Given a model parameter hessian, moment transform, and indices, compute the
hessian with respect to moment parameters.
""" ->
function get_moment_hess(model_params::Vector{Float64},
	                       model_grad::Vector{Float64}, model_hess::Matrix{Float64},
	                       indices::ModelIndices, transform::Function)
	# Get the necessary derivatives to do a change of variables in the Hessian.
	# A legend for the variable names to ease legibility:
 	# x = Moment parameters
 	# y = Model parameters
 	# 0 = Untransformed parameters
 	# 1 = Transformed parameters
 	# 2 = Constant parameters

	transformed_param_keys = keys(indices.trans_input);
	untransformed_param_keys = setdiff(keys(indices.model_params), transformed_param_keys);

	# The model indices of the transformed parameters in the correct order.
	transformed_indices = zeros(Int64, indices.num_transformed_indices)
	for trans_param in transformed_param_keys
		transformed_indices[indices.trans_input[trans_param]] =
			indices.model_params[trans_param]
	end

	trans_model_params = model_params[transformed_indices];
 	trans_moment_params =
		VariationalModelIndices.get_transformed_params(model_params, indices, transform)

	dx_dy_func = get_dx_dy_func(transform, indices.num_transformed_indices);
	d2x_dy2_funcs = get_d2x_dy2_funcs(transform, indices.num_transformed_indices);

	dx_dy = dx_dy_func(trans_model_params);

	d2x_dy2 = get_d2x_dy2(d2x_dy2_funcs, trans_model_params);

	# Rename for clarity.
	full_df_dy = model_grad;
	full_d2f_dy2 = model_hess;

	df_dy_1 = full_df_dy[transformed_indices];

	d2f_dy2_11 = full_d2f_dy2[transformed_indices, transformed_indices];
	d2f_dy2_1A = full_d2f_dy2[transformed_indices, :]; # A for "all"

	# Compute the submatrix of the transformed Hessian corresponding to
	# the transformed parameters.
	d2f_dx2_11 = transform_hessian(dx_dy, d2x_dy2, df_dy_1, d2f_dy2_11);

	# Compute the submatrix of the transformed Hessian corresponding to
	# the mixed partials between the transformed and untransformed parameters.
	# Note that the columns are indexed in model order, where the rows
	# are indexed by indices.trans_output.
	d2f_dx2_1A = dx_dy' \ d2f_dy2_1A;
	d2f_dx2_A1 = d2f_dx2_1A';

	# Populate the hessian with respect to the moment parameters.
	moment_hess =
		zeros(Float64, indices.num_moment_indices, indices.num_moment_indices);
	for param0_row in untransformed_param_keys
		param0_row_indices = indices.model_params[param0_row]
		moment0_row_indices = indices.moment_params[param0_row]

		# The untransformed by untransformed block
		for param0_col in untransformed_param_keys
			param0_col_indices = indices.model_params[param0_col]
			moment0_col_indices = indices.moment_params[param0_col]
			moment_hess[moment0_row_indices, moment0_col_indices] =
				full_d2f_dy2[param0_row_indices, param0_col_indices]
		end

		# The untransformed by transformed blocks
		for param1 in keys(indices.trans_output)
			param1_indices = indices.trans_output[param1]
			moment1_indices = indices.moment_params[param1]
			moment_hess[moment0_row_indices, moment1_indices] =
				d2f_dx2_A1[param0_row_indices, param1_indices]
			moment_hess[moment1_indices, moment0_row_indices] =
				d2f_dx2_1A[param1_indices, param0_row_indices]
		end
	end

	# The transformed by transformed block
	for param1_row in keys(indices.trans_output)
		param1_row_indices = indices.trans_output[param1_row]
		moment1_row_indices = indices.moment_params[param1_row]
		for param1_col in keys(indices.trans_output)
			param1_col_indices = indices.trans_output[param1_col]
			moment1_col_indices = indices.moment_params[param1_col]
			moment_hess[moment1_row_indices, moment1_col_indices] =
				d2f_dx2_11[param1_row_indices, param1_col_indices]
		end
	end

	# Finally, get the submatrices corresponding to the "constant" parameters
	# (the priors and the data).
	const_param_hess =
		zeros(Float64, indices.num_const_indices, indices.num_moment_indices);
	const_params = keys(indices.model_const)
	@assert collect(keys(indices.model_const)) == collect(keys(indices.moment_const))

	for param2 in const_params
		param2_indices = indices.model_const[param2]
		moment2_indices = indices.moment_const[param2]
		for param0 in untransformed_param_keys
			param0_indices = indices.model_params[param0]
			moment0_indices = indices.moment_params[param0]
			const_param_hess[moment2_indices, moment0_indices] =
				full_d2f_dy2[param2_indices, param0_indices]
		end
		for param1 in keys(indices.trans_output)
			param1_indices = indices.trans_output[param1]
			moment1_indices = indices.moment_params[param1]
			const_param_hess[moment2_indices, moment1_indices] =
				d2f_dx2_A1[param2_indices, param1_indices]
		end
	end

	# The model indices of the constant parameters in the correct order.
	const_indices = zeros(Int64, indices.num_const_indices)
	for const_param in const_params
		const_indices[indices.moment_const[const_param]] = indices.model_const[const_param]
	end
	const_const_hess = full_d2f_dy2[const_indices, const_indices];

	moment_hess, const_param_hess, const_const_hess,
	dx_dy, d2x_dy2, full_df_dy, full_d2f_dy2, transformed_indices
end

end # module