text:
# run './mpiexecjl -n 4 julia --project scripts-part1/Diffusion_3D.jl' for multi-XPUs
const USE_GPU = false
using ParallelStencil
using ParallelStencil.FiniteDifferences3D
@static if USE_GPU
    @init_parallel_stencil(CUDA, Float64, 3)
else
    @init_parallel_stencil(Threads, Float64, 3)
end
using Plots, Printf, ImplicitGlobalGrid
import MPI

# macros to avoid array allocation
macro qx(ix,iy,iz)  esc(:( -D_dx*(Hτ[$ix+1,$iy+1,$iz+1] - Hτ[$ix  ,$iy+1,$iz+1]) )) end
macro qy(ix,iy,iz)  esc(:( -D_dy*(Hτ[$ix+1,$iy+1,$iz+1] - Hτ[$ix+1,$iy  ,$iz+1]) )) end
macro qz(ix,iy,iz)  esc(:( -D_dz*(Hτ[$ix+1,$iy+1,$iz+1] - Hτ[$ix+1,$iy+1,$iz  ]) )) end

@parallel_indices (ix,iy,iz) function compute_pseudo!(Hτ2, Hτ, H, D_dx, D_dy, D_dz, dt, dτ, _dx, _dy, _dz, size_H1_2, size_H2_2, size_H3_2)
    if (ix<=size_H1_2 && iy<=size_H2_2 && iz<=size_H3_2)
        Hτ2[ix+1,iy+1,iz+1] = Hτ[ix+1,iy+1,iz+1] - dτ*( (Hτ[ix+1,iy+1,iz+1] - H[ix+1,iy+1,iz+1])/dt + (@qx(ix+1,iy,iz)-@qx(ix,iy,iz))*_dx + (@qy(ix,iy+1,iz)-@qy(ix,iy,iz))*_dy + (@qz(ix,iy,iz+1)-@qz(ix,iy,iz))*_dz )
    end
    return
end

@parallel_indices (ix,iy,iz) function compute_residuals!(R_H, Hτ, H, D_dx, D_dy, D_dz, dt, _dx, _dy, _dz, size_H1_2, size_H2_2, size_H3_2)
    if (ix<=size_H1_2 && iy<=size_H2_2 && iz<=size_H3_2)
        R_H[ix,iy,iz] = -(Hτ[ix+1,iy+1,iz+1] - H[ix+1,iy+1,iz+1])/dt - (@qx(ix+1,iy,iz)-@qx(ix,iy,iz))*_dx - (@qy(ix,iy+1,iz)-@qy(ix,iy,iz))*_dy - (@qz(ix,iy,iz+1)-@qz(ix,iy,iz))*_dz
    end
    return
end

@views norm_g(A) = (sum2_l = sum(A.^2); sqrt(MPI.Allreduce(sum2_l, MPI.SUM, MPI.COMM_WORLD)))

@views function diffusion_3D(; do_visu=false)
    do_print = false
    # Physics
    Lx ,Ly ,Lz  = 10.0, 10.0, 10.0  # domain size
    Lxc,Lyc,Lzc = Lx/2, Ly/2, Lz/2  # domain center
    D           = 1.0               # diffusion coefficient
    ttot        = 1.0               # total simulation time
    dt          = 0.2               # physical time step

    # Numerics
    nx, ny, nz = 2^5, 2^5, 2^5      # local grid size
    ncheck     = 100                # check residuals every ncheck pseudo time steps
    iter_max   = 1e6                # max pseudo time iterations
    tol        = 1e-8               # tolerance for converging

    # Derived numerics
    me, dims = init_global_grid(nx, ny, nz)                             # Initialization of MPI and more...
    @static if USE_GPU select_device() end                              # select one GPU per MPI local rank (if >1 GPU per node)
      dx,   dy,   dz = Lx/nx_g(), Ly/ny_g(), Lz/nz_g()                  # local grid steps
     _dx,  _dy,  _dz = 1.0/dx, 1.0/dy, 1.0/dz
    D_dx, D_dy, D_dz = D/dx, D/dy, D/dz
    nx_v, ny_v, nz_v = (nx-2)*dims[1], (ny-2)*dims[2], (nz-2)*dims[3]   # global grid size without halo
                  dτ = min(dx, dy, dz)^2/D/8.0                          # pseudo time step
                  nt = cld(ttot, dt)                                    # total number of physical time steps
    Xind, Yind, Zind = Int(cld(Lxc, dx)), Int(cld(Lyc, dy)), Int(cld(Lzc, dz))

    # Array initialisation
    H   = @zeros(nx, ny, nz)
    H  .= Data.Array([2.0*exp(-((x_g(ix,dx,H)+dx/2 -Lx/2)^2+(y_g(iy,dy,H)+dy/2 -Ly/2)^2+(z_g(iz,dz,H)+dz/2 -Lz/2)^2)/2) for ix=1:size(H,1), iy=1:size(H,2), iz=1:size(H,3)])
    Hτ  = copy(H)
    Hτ2 = copy(H)
    R_H = @zeros(nx-2, ny-2, nz-2)          # residual without halo
    size_H1_2, size_H2_2, size_H3_2 = size(H,1)-2, size(H,2)-2, size(H,3)-2

    # Prepare for 2D visualization
    if do_visu
        xc_v, yc_v = LinRange(dx+dx/2, Lx-dx-dx/2, nx_v), LinRange(dy+dy/2, Ly-dy-dy/2, ny_v)
        H_v   = zeros(nx_v, ny_v, nz_v)     # global H without boundaries (for visualization)
        H_inn = zeros(nx-2, ny-2, nz-2)     # no halo local array
        if (me==0) ENV["GKSwstype"]="nul"; if isdir("plots_out1")==false mkdir("plots_out1") end; loadpath = "./plots_out1/"; anim = Animation(loadpath,String[]); println("Animation directory: $(anim.dir)") end
        if (nx_v*ny_v*sizeof(Data.Number) > 0.8*Sys.free_memory()) error("Not enough memory for visualization.") end
    end

    # Physical time stepping
    for it = 1:nt
        # update H using backward Euler
        H .= Hτ

        t_tic = 0.0; niter = 0
        # pseudo time stepping
        for iter =1:iter_max
            # start timer after warming
            if (iter==11) t_tic = Base.time(); niter = 0 end
            # update Hτ in pseudo time
            @hide_communication (16, 2, 2) begin
                @parallel compute_pseudo!(Hτ2, Hτ, H, D_dx, D_dy, D_dz, dt, dτ, _dx, _dy, _dz, size_H1_2, size_H2_2, size_H3_2)
                Hτ, Hτ2 = Hτ2, Hτ
                update_halo!(Hτ)
            end
            niter += 1
            # compute and check residuals
            if iter % ncheck == 0
                @parallel compute_residuals!(R_H, Hτ, H, D_dx, D_dy, D_dz, dt, _dx, _dy, _dz, size_H1_2, size_H2_2, size_H3_2)
                err = norm_g(Array(R_H))/(nx_v*ny_v*nz_v)
                if (me==0 && do_print) println(">> time = $(round(it*dt, sigdigits=3))s, iterations = $(iter), residual = $(round(err, sigdigits=3))") end
                if err < tol break; end
            end
        end

        # performance evaluation
        t_toc = Base.time() - t_tic
        A_eff = 3/1e9*nx_g()*ny_g()*nz_g()*sizeof(Float64)  # Effective main memory access per iteration [GB]
        t_it  = t_toc/niter                                 # Execution time per iteration [s]
        T_eff = A_eff/t_it                                  # Effective memory throughput [GB/s]
        if (me==0 && do_print) println("-----time = $(round(it*dt, sigdigits=3)) s, T_eff = $(round(T_eff, sigdigits=3)) GB/s-----") end

        # Visualization
        if do_visu && me==0
            H_inn .= Array(Hτ[2:end-1,2:end-1,2:end-1]); gather!(H_inn, H_v)
            opts = (aspect_ratio=1, xlims=(xc_v[1], xc_v[end]), ylims=(yc_v[1], yc_v[end]),xtickfontsize=8,ytickfontsize=8,clims=(0, 1.5), xlabel="Lx", ylabel="Ly",xguidefontsize=8,yguidefontsize=8)
            P1 = heatmap(xc_v, yc_v, Array(H_v[:,:,nx÷2])'              ; opts..., title="0.5Lz (t=$(round(it*dt, sigdigits=3))s)", titlefontsize=10)
            P2 = heatmap(xc_v, yc_v, Array(H_v[:,:,Int(ceil(128*0.3))])'; opts..., title="0.3Lz", titlefontsize=10)
            P3 = heatmap(xc_v, yc_v, Array(H_v[:,:,Int(ceil(128*0.7))])'; opts..., title="0.7Lz", titlefontsize=10)
            display(plot(P2,P1,P3, layout = (1, 3), legend = false, size=(1500,400))); frame(anim)
        end
    end

    # Save animation
    if (do_visu && me==0) gif(anim, "diffusion_3D.gif", fps = 5)  end

    # Return the global array xc_g and H_g
    xc_g = LinRange(dx/2, Lx-dx/2, nx_g())
    H_g  = zeros(nx*dims[1],ny*dims[2],nz*dims[3])
    if (me==0) gather!(Array(Hτ), H_g) end

    finalize_global_grid()
    return H_g, xc_g, H_g[Xind, Yind, Zind]
end

H_g, xc_g, Hc = diffusion_3D(; do_visu=false);


text:
module Org

export OrgDoc, @org_str,                        # types/org.jl, here
    org, term, tableofcontents, html, markdown, # render/*
    parsetree                                   # analysis/parsetree.jl

include("types/org.jl")
include("cache/cache.jl")

include("parse/interpret.jl")
include("parse/markdown.jl")

const org_str_typeflags =
    Dict("d" => OrgDoc,
         "e" => Element,
         "o" => Object,
         "h" => Heading,
         "s" => Section,
         "p" => Paragraph,
         "ts" => Timestamp,
         "m" => TextMarkup,
         "t" => TextPlain)

macro org_str(content::String, type::String="d")
    if !haskey(org_str_typeflags, type)
        throw(ArgumentError(
            string("Unknown Org type flag: $type",
                   "\n\nRecognised types flags:\n",
                   join(map(((short, type)::Pair -> " • $short, $type"),
                            collect(org_str_typeflags)), '\n'), '\n')))
    end
    parse(org_str_typeflags[type], content)
end

include("types/documentation.jl")
include("parse/documentation.jl")

include("render/utils.jl")
include("render/formatting.jl")
include("render/org.jl")
include("render/term.jl")
include("render/html.jl")
include("render/markdown.jl")

include("analysis/parsetree.jl")
include("analysis/diff.jl")

include("precompile.jl")

end

ext: jl


text:
@testset "959.regions-cut-by-slashes.jl" begin
    @test regions_by_slashes([
        "\\/",
        "/\\"
      ]) == 4
    @test regions_by_slashes([
        "/\\",
        "\\/"
      ]) == 5
    @test regions_by_slashes([
        "//",
        "/ "
      ]) == 3
    @test regions_by_slashes([
        " /",
        "  "
      ]) == 1
    @test regions_by_slashes([
        " /",
        "/ "
      ]) == 2
end


text:
using JuliaToJupyter
using Test

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


text:
#!/usr/bin/env julia
include("commandUtils.jl")
activatePkg()
using ArgParse
using JLD
using Distributed
s = ArgParseSettings(description="Calculates field based on boundary saved in file.")
@add_arg_table s begin
    "--simFile", "-s"
        help = "Simulation file to include, must set up boundary named simBoundary, can also edit or set k or other params to avoid the need to
                set in command arguments."
        required = true
    "--boundaryFile", "-b"
        help = "File with calculated kf and field xf"
        required = true
    "--dx"
        help = "Plot dx"
        arg_type = Float64
        default = 0.01
    "--dy"
        help = "Plot dy"
        arg_type = Float64
        default = 0.01
    "--xRange"
        help = "X range (tuple)"
        default = "(-1.5, 1.5)"
    "--yRange"
        help = "Y range (tuple)"
        default = "(-1.5, 1.5)"
    "--procs", "-p"
        help = "Number of processes to use in parallel"
        arg_type = Int
        default = convert(Int, Sys.CPU_THREADS/2)
end

arguments = parse_args(s)
@debug "Arguments" arguments

addprocs(arguments["procs"])

using BoundaryIntegralMethod

boundaryFile = arguments["boundaryFile"]
xRange = eval(Meta.parse(arguments["xRange"]))
yRange = eval(Meta.parse(arguments["yRange"]))
dx = arguments["dx"]
dy = arguments["dy"]

simFile = arguments["simFile"]
simName = split(basename(simFile), ".")[1]
if (!startswith(simFile, "/"))
  simFile = "$(pwd())/$simFile"
end
include(simFile)
b = simBoundary

(k, xf) = jldopen(boundaryFile, "r") do file
  (read(file, "k/full"), read(file, "boundary/full"))
end

n = convert(Int,size(xf,1)/4)

kForFile = string(real(k),"_",imag(k))
outFile = "output/field-$simName-$n-$kForFile.h5"
mkpath("output")
exportTotalField(b, xf, k; outFile=outFile, xRange=xRange, yRange=yRange, dx=dx, dy=dy)


text:
using TOML

using Test
using Dates

const jsnval = Dict{String,Function}(
    "string" =>identity,
    "float"    => (s -> Base.parse(Float64, s)),
    "integer"  => (s -> Base.parse(Int64, s)),
    "datetime" => (s -> Base.parse(DateTime, s, dateformat"yyyy-mm-ddTHH:MM:SSZ")),
    "array"    => (a -> map(jsn2data, a)),
    "bool"     => (b -> b == "true")
)

function jsn2data(jsn)
    if "type" in keys(jsn)
        jsnval[jsn["type"]](jsn["value"])
    elseif jsn isa Vector
        [jsn2data(v) for v in jsn]
    else
        Dict{String,Any}([k => jsn2data(v) for (k, v) in jsn])
    end
end

#########
# Valid #
#########

valid_test_folder = joinpath(@__DIR__, "testfiles", "valid")

function check_valid(f)
    fp = joinpath(valid_test_folder, f)
    jsn = jsn2data(@eval include($fp * ".jl"))
    tml = TOML.parsefile(fp * ".toml")
    return isequal(tml, jsn)
end

@testset "valid" begin

@test check_valid("array-empty")
@test check_valid("array-nospaces")
@test check_valid("array-string-quote-comma-2")
@test check_valid("array-string-quote-comma")
@test check_valid("array-string-with-comma")
@test check_valid("array-table-array-string-backslash")
@test check_valid("arrays-hetergeneous")
@test check_valid("arrays-nested")
@test check_valid("arrays")
@test check_valid("bool")
@test check_valid("comments-at-eof")
@test check_valid("comments-at-eof2")
@test check_valid("comments-everywhere")
@test_broken check_valid("datetime-timezone")
@test_broken check_valid("datetime")
@test check_valid("double-quote-escape")
@test check_valid("empty")
@test check_valid("escaped-escape")
@test check_valid("example")
@test check_valid("exponent-part-float")
@test check_valid("float-exponent")
@test check_valid("float-underscore")
@test check_valid("float")
@test check_valid("implicit-and-explicit-after")
@test check_valid("implicit-and-explicit-before")
@test check_valid("implicit-groups")
@test check_valid("inline-table-array")
@test check_valid("inline-table")
@test check_valid("integer-underscore")
@test check_valid("integer")
@test check_valid("key-equals-nospace")
@test check_valid("key-numeric")
@test check_valid("key-space")
@test check_valid("key-special-chars")
@test check_valid("keys-with-dots")
@test check_valid("long-float")
@test check_valid("long-integer")
@test check_valid("multiline-string")
@test check_valid("nested-inline-table-array")
@test check_valid("newline-crlf")
@test check_valid("newline-lf")
if Sys.iswindows() &&
    # Sometimes git normalizes the line endings
    contains(read(joinpath(valid_test_folder, "raw-multiline-string-win.toml"), String), '\r')
    @test check_valid("raw-multiline-string-win")
else
    @test check_valid("raw-multiline-string")
end
@test check_valid("raw-string")
@test check_valid("right-curly-brace-after-boolean")
@test check_valid("string-empty")
@test check_valid("string-escapes")
@test check_valid("string-nl")
@test check_valid("string-simple")
@test check_valid("string-with-pound")
@test check_valid("table-array-implicit")
@test check_valid("table-array-many")
@test check_valid("table-array-nest")
@test check_valid("table-array-one")
@test check_valid("table-array-table-array")
@test check_valid("table-empty")
@test check_valid("table-no-eol")
@test check_valid("table-sub-empty")
@test check_valid("table-whitespace")
@test check_valid("table-with-literal-string")
@test check_valid("table-with-pound")
@test check_valid("table-with-single-quotes")
@test check_valid("underscored-float")
@test check_valid("underscored-integer")
@test check_valid("unicode-escape")
@test check_valid("unicode-literal")

end

###########
# Invalid #
###########

invalid_test_folder = joinpath(@__DIR__, "testfiles", "invalid")

# TODO: Check error type
function check_invalid(f)
    fp = joinpath(invalid_test_folder, f)
    tml = TOML.tryparsefile(fp * ".toml")
    return tml isa TOML.Internals.ParserError
end

@test check_invalid("datetime-malformed-no-leads")
@test check_invalid("datetime-malformed-no-secs")
@test check_invalid("datetime-malformed-no-t")
@test check_invalid("datetime-malformed-with-milli")
@test check_invalid("duplicate-key-table")
@test check_invalid("duplicate-keys")
@test check_invalid("duplicate-tables")
@test check_invalid("empty-implicit-table")
@test check_invalid("empty-table")
@test check_invalid("float-leading-zero-neg")
@test check_invalid("float-leading-zero-pos")
@test check_invalid("float-leading-zero")
@test check_invalid("float-no-leading-zero")
@test check_invalid("float-no-trailing-digits")
@test check_invalid("float-underscore-after-point")
@test check_invalid("float-underscore-after")
@test check_invalid("float-underscore-before-point")
@test check_invalid("float-underscore-before")
@test check_invalid("inline-table-linebreak")
@test check_invalid("integer-leading-zero-neg")
@test check_invalid("integer-leading-zero-pos")
@test check_invalid("integer-leading-zero")
@test check_invalid("integer-underscore-after")
@test check_invalid("integer-underscore-before")
@test check_invalid("integer-underscore-double")
@test check_invalid("key-after-array")
@test check_invalid("key-after-table")
@test check_invalid("key-empty")
@test check_invalid("key-hash")
@test check_invalid("key-newline")
@test check_invalid("key-no-eol")
@test check_invalid("key-open-bracket")
@test check_invalid("key-single-open-bracket")
@test check_invalid("key-space")
@test check_invalid("key-start-bracket")
@test check_invalid("key-two-equals")
@test check_invalid("llbrace")
@test check_invalid("multi-line-inline-table")
@test check_invalid("multi-line-string-no-close")
@test check_invalid("rrbrace")
@test check_invalid("string-bad-byte-escape")
@test check_invalid("string-bad-codepoint")
@test check_invalid("string-bad-escape")
@test check_invalid("string-bad-slash-escape")
@test check_invalid("string-bad-uni-esc")
@test check_invalid("string-byte-escapes")
@test check_invalid("string-no-close")
@test check_invalid("table-array-implicit")
@test check_invalid("table-array-malformed-bracket")
@test check_invalid("table-array-malformed-empty")
@test check_invalid("table-empty")
@test check_invalid("table-nested-brackets-close")
@test check_invalid("table-nested-brackets-open")
@test check_invalid("table-whitespace")
@test check_invalid("table-with-pound")
@test check_invalid("text-after-array-entries")
@test check_invalid("text-after-integer")
@test check_invalid("text-after-string")
@test check_invalid("text-after-table")
@test check_invalid("text-before-array-separator")
@test check_invalid("text-in-array")


text:
using Test, FreeCell

@test true


text:
# load ImpvolEquilibrium first so that methods are accessible
include("../../equilibrium.jl")
using ImpvolEquilibrium, FileIO, JLD2

## these are needed for data -> parameters mapping
parameters = Dict{Symbol, Any}()

# CES parameters
parameters[:sigma] = 0.999
parameters[:theta] = 2.0
# we also have to change eta so that Gamma function is finite. eta only affects overall level of productivity
parameters[:eta] = 2.0

########## parameters common across scenarios
## these are function of data
# inverse of adjustment cost, 0 if cannot readjust
parameters[:one_over_rho] = 0.0

include("../config.jl")
# change parameters after reading data, but common across scenarios
jldopen("common_parameters.jld2", "w") do file
	file["parameters"] = parameters
end


text:
using BrainDecoder
using HDF5
using Test

@testset "BrainDecoder.jl" begin

    # BData
    voxel_data = [0.0   1.0   2.0   3.0   4.0   5.0   6.0   7.0   8.0   9.0
                  10.0  11.0  12.0  13.0  14.0  15.0  16.0  17.0  18.0  19.0
                  20.0  21.0  22.0  23.0  24.0  25.0  26.0  27.0  28.0  29.0
                  30.0  31.0  32.0  33.0  34.0  35.0  36.0  37.0  38.0  39.0
                  40.0  41.0  42.0  43.0  44.0  45.0  46.0  47.0  48.0  49.0
                  50.0  51.0  52.0  53.0  54.0  55.0  56.0  57.0  58.0  59.0]
    runs = reshape([1, 1, 1, 2, 2, 2], 6, 1)
    stimulus_labels = ["stimulus_02"; "stimulus_01"; "stimulus_03"; "stimulus_01"; "stimulus_03"; "stimulus_02"]

    bdata = BData("./data/test.h5")

    @test bdata.select("VoxelData") == voxel_data
    @test bdata.select("Run") == runs
    @test bdata.get_labels("stimulus_name") == stimulus_labels
    @test bdata.get_label("stimulus_name") == stimulus_labels

    # Features
    # TODO: add better tests
    feat_dir = "/home/nu/data/contents_shared/ImageNetTest/derivatives/features/caffe/bvlc_alexnet"
    layers = ["conv1", "conv2", "conv3", "conv4", "conv5", "fc6", "fc7", "fc8", "norm1", "norm2", "pool1", "pool2", "pool5", "prob", "relu1", "relu2", "relu3", "relu4", "relu5", "relu6", "relu7"]
    labels = ["n01443537_22563", "n01621127_19020", "n01677366_18182", "n01846331_17038", "n01858441_11077", "n01943899_24131", "n01976957_13223", "n02071294_46212", "n02128385_20264", "n02139199_10398", "n02190790_15121", "n02274259_24319", "n02416519_12793", "n02437136_12836", "n02437971_5013", "n02690373_7713", "n02797295_15411", "n02824058_18729", "n02882301_14188", "n02916179_24850", "n02950256_22949", "n02951358_23759", "n03064758_38750", "n03122295_31279", "n03124170_13920", "n03237416_58334", "n03272010_11001", "n03345837_12501", "n03379051_8496", "n03452741_24622", "n03455488_28622", "n03482252_22530", "n03495258_9895", "n03584254_5040", "n03626115_19498", "n03710193_22225", "n03716966_28524", "n03761084_43533", "n03767745_109", "n03941684_21672", "n03954393_10038", "n04210120_9062", "n04252077_10859", "n04254777_16338", "n04297750_25624", "n04387400_16693", "n04507155_21299", "n04533802_19479", "n04554684_53399", "n04572121_3262"]

    features = Features(feat_dir)

    @test features.layers == layers
    @test features.labels == labels
    @test features.get_labels() == labels
    @test size(features.get_features("fc8")) == (50, 1000)
    @test size(features.get_features("conv5")) == (50, 256, 13, 13)
    #@test size(features.get_features("conv1")) == (50, 96, 55, 55)

    # save_array
    data = [1 2 3 4 5
            6 7 8 9 10]
    save_array("test_save_array_dense.mat", "data", data)

    data = [1 0 0 0
            2 2 0 0
            3 3 3 0]
    save_array("test_save_array_sparse_2d.mat", "data", data, sparse=true)

end


text:
function save_merafile(dataobject::DataSetType, filename::String;
                    fmode::String="add",
                    ctype::String="blosc",
                    cn::Int=0,
                    comment::String="",
                    merafile_version::Float64=1.,
                    printfields::Bool=false,
                    verbose::Bool=verbose_mode)

    return save_merafile(dataobject,
                    filename=filename,
                    fmode=fmode,
                    ctype=ctype, cn=cn,
                    comment=comment, merafile_version=merafile_version,
                    printfields=printfields, verbose=verbose)
end

function save_merafile( dataobject::DataSetType;
                    filename::String="",
                    fmode::String="add",
                    ctype::String="blosc",
                    cn::Int=0,
                    comment::String="",
                    merafile_version::Float64=1.,
                    printfields::Bool=false,
                    verbose::Bool=verbose_mode)

                    ###filemode = "w" #r, r+, cw, w
                    ###ctype = "blosc"
                    ###cn = 9 # compression strength: [0,9]

    datasource, use_descriptor, descriptor_names = check_datasource(dataobject)
    filemode, overwdata = check_merafile_mode(fmode, datasource, filename)

    fields = propertynames(dataobject)
    mem = usedmemory(dataobject, false)
    count_fields = 0
    column_names = propertynames(dataobject.data.columns)

    if printfields
        println()
        println("Legend:")
        rr =96; gg=168; bb=48 # green
        print("- type in ", "\e[38;2;$rr;$gg;$bb;249m","green", "\e[0m", " => changed type to store it \n")

        rr =175; gg=120; bb=2 # orange
        print("- type in ", "\e[38;2;$rr;$gg;$bb;249m","orange", "\e[0m", " => marks dictionaries \n")

        rr =200; gg=75; bb=75 # red
        print("- type in ", "\e[38;2;$rr;$gg;$bb;249m","red", "\e[0m", " => skipped undefined field \n")
        println("=========================================")
    end

    # Main loop
    h5open(filename, filemode) do file
        count_fields = 0

        # add/overwrite top-level for given datatype
        file = prep_toplevel(file, datasource, overwdata)

        # create top-level (zero-level) for datatype and add header
        group_header, count_fields = groupheader(file, datasource, comment, mem, merafile_version, ctype, cn, count_fields)

        # add fields and data
        for i in fields
            count_fields += 1
            field_levels = 1

            # split and store juliaDB columns into arrays
            # save column names (Symbol) into a String
            if i == :data
                count_fields = store_database(dataobject, i, field_levels, count_fields, printfields, group_header, ctype, cn)
            else
                if isdefined(dataobject, i) # only store defined fields of top-level
                    # create top-level (zero-level) fields and store their related content
                    group_a, a_fields, a, count_fields  = store_fields(dataobject, i, printfields, field_levels, count_fields, group_header)

                    if length(a_fields) != 0 && typeof(a) != Dict{Any,Any}
                        field_levels = 2
                        count_fields   = store_subfields(dataobject, i, printfields, field_levels, count_fields, group_a, a_fields, a)
                    end

                else # field not defined
                    if printfields
                        r =200; g=75; b=75 # red
                        print("\e[1m","\t \t l", field_levels, ": ", "\e[0m", i, " = ", "\e[38;2;$r;$g;$b;249m","undefined", "\e[0m \n")
                    end
                end # isdefined
            end
        end
        count_fields += 1
        group_header["total_fields"] = count_fields
    end #h5open

    if printfields
        println()
        println()
    end

    if verbose
        println("Directory: ", dataobject.info.path )
        println("merafile_version: ", merafile_version)
        println("Simulation code: ", dataobject.info.simcode)
        println("DataType: ", datasource)
        println("Data variables: ", column_names)
        if use_descriptor
             println("Descriptor: ", descriptor_names)
        end
        println("Total fields: ", count_fields)
        println("Writing mode: ", fmode)
        println("Compression type: ", ctype, " with: ", cn )
        mem = usedmemory(dataobject, false)
        println("Uncompressed memory size: ", round(mem[1], digits=3)," ", mem[2])
        s = filesize(filename)
        svalue, sunit = humanize(Float64(s), 3, "memory")
        println("File size: ", svalue, " ", sunit)
        println("-----------------------------------")
        println()
    end

    return
end

function check_datasource(dataobject::DataSetType)
    if typeof(dataobject) == HydroDataType
        datasource = "hydro"
        use_descriptor = dataobject.info.descriptor.usehydro
        descriptor_names = dataobject.info.descriptor.hydro
    #elseif typeof(dataobject) == GravDataType
    #    datasource = "gravity"
    #    use_descriptor = dataobject.info.descriptor.usegravity
    #    descriptor_names = dataobject.info.descriptor.gravity
    elseif typeof(dataobject) == PartDataType
        datasource = "particles"
        use_descriptor = dataobject.info.descriptor.useparticles
        descriptor_names = dataobject.info.descriptor.particles
    elseif typeof(dataobject) == ClumpDataType
        datasource = "clumps"
        use_descriptor = dataobject.info.descriptor.useclumps
        descriptor_names = dataobject.info.descriptor.clumps
    end

    return datasource, use_descriptor, descriptor_names
end

function check_merafile_mode(fmode::String, datasource::String, filename::String)
    println()
    if !isfile(filename)
        println("Create mera-file:")
    else
        println("Existing mera-file:")
    end

    if fmode == "add"
        filemode = "cw"
        overwdata = false
        println("Adding $datasource data...")
    elseif fmode == "add+"
        filemode = "cw"
        overwdata = true
        #println("Adding/overwriting $datasource data...")
    elseif fmode == "w"
        filemode = "w"
        overwdata = true
        println("Creating/Overwriting mera-file and adding $datasource data...")
    else
        error("""Choose between "w" = create/overwrite file and add data, \n use "add" to add data to an existing file or it will be created. """)
    end

    return filemode, overwdata

end

function prep_toplevel(file::HDF5.HDF5File, datasource::String, overwdata::Bool)
    if exists(file, datasource)
        if overwdata == true
            println("""Overwriting existing data-type "$datasource" with new data...""")
            o_delete(file, datasource)
        else
            println("""Stop adding: found existing data-type "$datasource" """)
            error("""Use fmode="add+" to overwrite existing data-type or fmode="w" to overwrite file""")
        end
    else
        if overwdata == true
            println("Adding $datasource data...")
        end
    end
    println()

    return file
end

function groupheader(file::HDF5.HDF5File, datasource::String, comment::String, mem::Tuple{Float64,String}, merafile_version::Float64, ctype::String, cn::Int, count_fields::Int)
    group_header = g_create(file, datasource)
    group_header["comment"] = comment
    group_header["memory"]  = string.( collect(mem) ) # convert tuple to array
    group_header["merafile_version"] = merafile_version
    group_header["compression_type"] = ctype
    group_header["compression"] = cn
    count_fields += 5
    return group_header, count_fields
end

# convert type  if necessary
# atype = typeof(a)
# a = content
function changetype(atype, a)
    achange = false
    if atype == Bool
        a= string(a); achange=true
    elseif atype == Array{Symbol,1}
        a= string.(a); achange=true
    elseif atype == Array{Bool,1}
        a= string.(a); achange=true
    elseif atype == DateTime
        a = string.(a); achange=true
    end
    return a, achange
end

# split and store juliaDB columns into arrays
# save column names (Symbol) into a String
function store_database(dataobject::DataSetType, i, field_levels::Int, count_fields::Int, printfields::Bool, group_header, ctype::String, cn::Int)

    if printfields println("\e[1m","l", field_levels, ": ","\e[0m", i, " = JuliaDB" ) end

    # create fieldname
    count_fields += 1
    group_data = g_create(group_header, string(i))

    # remember column names of JuliaDB data table
    # convert tuple of symbols to array containing strings
    count_fields += 1
    column_names = propertynames(dataobject.data.columns)
    group_data["names"] = [s for s in string.(column_names)]

    # store each column of the data table in an array
    for n in column_names
        count_fields += 1
        group_data[string(n), ctype, cn] = select(dataobject.data, n)
    end

    return count_fields
end

# store top-level fields
function store_fields(dataobject::DataSetType, i, printfields::Bool, field_levels::Int, count_fields::Int, group_header)

    a = getfield(dataobject, i)
    a_fields = propertynames(a)
    atype = typeof(a)

    if printfields
        if atype != Dict{Any,Any}
            print("\e[1m","l", field_levels, ": ", "\e[0m", i, " = ", atype)
        else
            rr =175; gg=120; bb=2 # orange
            print("\e[1m","l", field_levels, ": ", "\e[0m", i, " = ", "\e[38;2;$rr;$gg;$bb;249m",atype, "\e[0m")
        end
    end

    count_fields += 1
    if length(a_fields) == 0 && atype != Dict{Any,Any}
        a, achange = changetype(atype, a)

        group_header[string(i)] = a
        if printfields
            if achange == false
                print(" => ", typeof(a), "\n")  # stored as type
            else
                rr =96; gg=168; bb=48 # green
                print(" => ", "\e[38;2;$rr;$gg;$bb;249m", typeof(a), "\e[0m", "\n")  # stored as type
            end
        end
    elseif atype == Dict{Any,Any}

        if length(a) == 0
            group_header[string(i)] = 0
            if printfields print(" => ",  "\e[1m", 0, "\e[0m", "\n") end # stored as
        else
            count_fields = store_dict(dataobject, i, printfields, field_levels, count_fields, group_header)
        end

    end
    return group_header, a_fields, a, count_fields
end

function store_subfields(dataobject, i, printfields::Bool, field_levels::Int, count_fields::Int, group_header, a_fields, a)

        count_fields += 1
        group_a = g_create(group_header, string(i))

    if printfields print("\n ") end
    for j in a_fields
        count_fields += 1

        if isdefined(a, j) # only store defined fields
            b = getfield(a, j) # field_levels == 2
            b_fields = propertynames(b)
            btype = typeof(b)
            if printfields
                if btype != Dict{Any,Any}
                    print("\e[1m","\t l", field_levels, ": ", "\e[0m", j, " = ", btype)
                else
                    rr =175; gg=120; bb=2 # orange
                    print("\e[1m","\t l", field_levels, ": ", "\e[0m", j, " = ", "\e[38;2;$rr;$gg;$bb;249m",btype, "\e[0m")
                end
            end

            bchange = false
            if length(b_fields) == 0 && btype != Dict{Any,Any}
                b, bchange = changetype(btype, b)
                group_a[string(j)] = b
                if printfields
                    if bchange == false
                        print(" => ", typeof(b), "\n")  # stored as type

                    else
                        rr =96; gg=168; bb=48 # green
                        print(" => ", "\e[38;2;$rr;$gg;$bb;249m", typeof(b), "\e[0m", "\n")  # stored as type
                    end
                end

            elseif btype == Dict{Any,Any}
                if length(b) == 0
                    group_a[string(j)] = 0
                    if printfields print(" => ",  "\e[1m", 0, "\e[0m", "\n") end # stored as
                else
                    if printfields print("\n") end
                    afield =getfield(dataobject,i)
                    count_fields = store_dict(afield, j, printfields, field_levels, count_fields, group_a)
                end

            elseif btype == DateTime
                b, bchange = changetype(btype, b)

                group_a[string(j)] = b
                if printfields
                    if bchange == false
                        print(" => ", typeof(b), "\n")  # stored as type
                    else
                        rr =96; gg=168; bb=48 # green
                        print(" => ", "\e[38;2;$rr;$gg;$bb;249m", typeof(b), "\e[0m", "\n")  # stored as type
                    end
                end

            else
                field_levels = 3
                count_fields = store_subsubfields(dataobject, j, printfields, field_levels, count_fields, group_a, b_fields, b)
                field_levels = 2

            end

        else # field not defined
            if printfields
                rr =200; gg=75; bb=75 # red
                print("\e[1m","\t \t l", field_levels, ": ", "\e[0m", j, " = ", "\e[38;2;$rr;$gg;$bb;249m","undefined", "\e[0m \n")
            end
        end # isdefined

    end
    return  count_fields
end

function store_subsubfields(dataobject, i, printfields::Bool, field_levels::Int, count_fields::Int, group_header, b_fields, b)

        count_fields += 1
        group_b = g_create(group_header, string(i))

    if printfields print("\n ") end
    for j in b_fields
        count_fields += 1

        if isdefined(b, j) # only store defined fields
            c = getfield(b, j) # field_levels == 3
            c_fields = propertynames(c)
            ctype = typeof(c)
            if printfields
                if ctype != Dict{Any,Any}
                    print("\e[1m","\t \t l", field_levels, ": ", "\e[0m", j, " = ", ctype)
                else
                    rr =175; gg=120; bb=2 #orange
                    print("\e[1m","\t \t l", field_levels, ": ", "\e[0m", j, " = ", "\e[38;2;$rr;$gg;$bb;249m",ctype, "\e[0m")
                end
            end

            cchange = false
            if length(c_fields) == 0 && ctype != Dict{Any,Any}
                c, cchange = changetype(ctype, c)
                group_b[string(j)] = c
                if printfields
                    if cchange == false
                        print(" => ", typeof(c), "\n")  # stored as type
                    else
                        rr =96; gg=168; bb=48 # green
                        print(" => ", "\e[38;2;$rr;$gg;$bb;249m", typeof(c), "\e[0m", "\n")  # stored as type
                        #print(" => ", "\e[1m", typeof(c), "\e[0m", "\n")  # stored as type
                    end
                end
            elseif ctype == Dict{Any,Any}
                if length(c) == 0
                    group_b[string(j)] = 0
                    if printfields print(" => ",  "\e[1m", 0, "\e[0m", "\n") end # stored as
                else
                    if printfields print("\n") end
                # store dict level 3 place holder
                #    afield =getfield(dataobject,i)
                #    store_dict(afield, j, printfields, field_levels, count_fields, group_a)
                end

            else
                if printfields print("\n ") end
                # subfield level 4 placeholder

            end

        else # field not defined
            if printfields
                rr =200; gg=75; bb=75 # red
                print("\e[1m","\t \t l", field_levels, ": ", "\e[0m", j, " = ", "\e[38;2;$rr;$gg;$bb;249m","undefined", "\e[0m \n")
            end
        end # isdefined

    end

    return  count_fields

end

function store_dict(dataobject, i, printfields::Bool, field_levels::Int, count_fields::Int, group_header)

    dictlist = getfield(dataobject, i)
    dictkeys = string.(keys(dictlist))

    if !exists(group_header, string(i))
        count_fields += 1
        group_a = g_create(group_header, string(i))
    else
        group_a = group_header[string(i)]
    end

    if printfields
        if field_levels == 1
            print("\t => keys: ", length(dictkeys), "\n")
        elseif field_levels ==2
            print("\t \t => keys: ",   length(dictkeys), "\e[0m", "\n")
        elseif field_levels ==3
            print("\t \t \t => keys: ",  length(dictkeys), "\n")
        end
    end

    for j in dictkeys
        if length(j) != 0
            if printfields
                println()
                if field_levels == 1
                    print("\t ", j ,"\n")
                elseif field_levels ==2
                    print("\t \t ", j, "\n")
                elseif field_levels ==3
                    print("\t \t \t ", j, "\n")
                end
            end

            count_fields += 1
            group_b = g_create(group_a, string(j)) # create header of PARAMS rubrique

            a = dictlist[j]
            atype = typeof(a)

            if atype != Dict{Any,Any} # store values if no sub-dict
                a, achange = changetype(atype, a)
                count_fields += 1
                group_b[string(j)] = a

            else #store falues if sub-dict
                # store sub-dict e.g. parameters of namelist
                subdict = dictlist[j]
                subtype = typeof(subdict)
                subdictkeys = string.(keys(subdict))

                for k in subdictkeys
                    count_fields += 1
                    group_c = g_create(group_b, string(k)) # create header of PARAMS rubrique

                    subdict[k], achange = changetype(subtype, subdict[k])
                    group_c[string(k)] = subdict[k] # write variable with corresponding parameters
                    if printfields
                        if field_levels == 1
                            print("\t ", k , " = SubDict{Any,Any} => ", typeof(subdict[k]), "\n")
                        elseif field_levels ==2
                            print("\t \t ", k, " = SubDict{Any,Any} => ", typeof(subdict[k]), "\n")
                        elseif field_levels ==3
                            print("\t \t \t", " = SubDict{Any,Any} => ", typeof(subdict[k]), "\n")
                        end
                    end
                end

            end

        end # skip if key is zero length
    end

    return count_fields
end


text:
# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.

import Logging: min_enabled_level, shouldlog, handle_message
using Base.CoreLogging, Logging, Printf

struct TimeLogger <: AbstractLogger
    initial_time::Float64
    file::Union{Nothing, IOStream}
    screen_log_level
    io_log_level
end

function TimeLogger(;
                    initial_time::Float64,
                    file::Union{Nothing, IOStream} = nothing,
                    screen_log_level = CoreLogging.Info,
                    io_log_level = CoreLogging.Info,
                   ) :: TimeLogger
    return TimeLogger(initial_time, file, screen_log_level, io_log_level)
end

min_enabled_level(logger::TimeLogger) = logger.io_log_level
shouldlog(logger::TimeLogger, level, _module, group, id) = true

function handle_message(logger::TimeLogger,
                        level,
                        message,
                        _module,
                        group,
                        id,
                        filepath,
                        line;
                        kwargs...)
    elapsed_time = time() - logger.initial_time
    time_string = @sprintf("[%12.3f] ", elapsed_time)
    if level >= logger.screen_log_level
        print(time_string)
        println(message)
    end
    if logger.file !== nothing && level >= logger.io_log_level
        write(logger.file, time_string)
        write(logger.file, message)
        write(logger.file, "\n")
        flush(logger.file)
    end
end

export TimeLogger


text:
module jlgr

import GR

@static if VERSION >= v"0.7.0-DEV.3476"
    using Serialization
    using Sockets
end

const None = Union{}

@static if VERSION < v"0.7.0-DEV.3137"
    const Nothing = Void
end

@static if VERSION < v"0.7.0-DEV.3155"
    const popfirst! = shift!
end

@static if VERSION >= v"0.7.0-DEV.3272"
    function search(s::AbstractString, c::Char)
        result = findfirst(isequal(c), s)
        result != nothing ? result : 0
    end
end

@static if VERSION < v"0.7.0-DEV.4534"
    reverse(a::AbstractArray; dims=nothing) =
        dims===nothing ? Base.reverse(a) : Base.flipdim(a, dims)
end

macro _tuple(t)
    :( Tuple{$t} )
end

const PlotArg = Union{AbstractString, AbstractVector, AbstractMatrix, Function}

const gr3 = GR.gr3

const plot_kind = [:line, :scatter, :stem, :hist, :contour, :contourf, :hexbin, :heatmap, :wireframe, :surface, :plot3, :scatter3, :imshow, :isosurface, :polar, :trisurf, :tricont]

const arg_fmt = [:xys, :xyac, :xyzc]

const kw_args = [:accelerate, :alpha, :backgroundcolor, :color, :colormap, :figsize, :isovalue, :labels, :location, :nbins, :rotation, :size, :tilt, :title, :xflip, :xlabel, :xlim, :xlog, :yflip, :ylabel, :ylim, :ylog, :zflip, :zlabel, :zlim, :zlog]

const colors = [
    [0xffffff, 0x000000, 0xff0000, 0x00ff00, 0x0000ff, 0x00ffff, 0xffff00, 0xff00ff] [0x282c34, 0xd7dae0, 0xcb4e42, 0x99c27c, 0x85a9fc, 0x5ab6c1, 0xd09a6a, 0xc57bdb] [0xfdf6e3, 0x657b83, 0xdc322f, 0x859900, 0x268bd2, 0x2aa198, 0xb58900, 0xd33682] [0x002b36, 0x839496, 0xdc322f, 0x859900, 0x268bd2, 0x2aa198, 0xb58900, 0xd33682]
    ]

const distinct_cmap = [ 0, 1, 984, 987, 989, 983, 994, 988 ]

@static if VERSION > v"0.7-"
  function linspace(start, stop, length)
    range(start, stop=stop, length=length)
  end
  repmat(A::AbstractArray, m::Int, n::Int) = repeat(A::AbstractArray, m::Int, n::Int)
end

mutable struct PlotObject
  obj
  args
  kvs
end

function Figure(width=600, height=450)
    obj = Dict()
    args = @_tuple(Any)
    kvs = Dict()
    kvs[:size] = (width, height)
    kvs[:ax] = false
    kvs[:subplot] = [0, 1, 0, 1]
    kvs[:clear] = true
    kvs[:update] = true
    PlotObject(obj, args, kvs)
end

function gcf()
    plt.kvs
end

plt = Figure()
ctx = Dict()
scheme = 0
background = 0xffffff

isrowvec(x::AbstractArray) = ndims(x) == 2 && size(x, 1) == 1 && size(x, 2) > 1

isvector(x::AbstractVector) = true
isvector(x::AbstractMatrix) = size(x, 1) == 1

function set_viewport(kind, subplot)
    mwidth, mheight, width, height = GR.inqdspsize()
    if haskey(plt.kvs, :figsize)
        w = 0.0254 *  width * plt.kvs[:figsize][1] / mwidth
        h = 0.0254 * height * plt.kvs[:figsize][2] / mheight
    else
        dpi = width / mwidth * 0.0254
        if dpi > 200
            w, h = [x * dpi / 100 for x in plt.kvs[:size]]
        else
            w, h = plt.kvs[:size]
        end
    end
    viewport = zeros(4)
    vp = copy(float(subplot))
    if w > h
        ratio = float(h) / w
        msize = mwidth * w / width
        GR.setwsviewport(0, msize, 0, msize * ratio)
        GR.setwswindow(0, 1, 0, ratio)
        vp[3] *= ratio
        vp[4] *= ratio
    else
        ratio = float(w) / h
        msize = mheight * h / height
        GR.setwsviewport(0, msize * ratio, 0, msize)
        GR.setwswindow(0, ratio, 0, 1)
        vp[1] *= ratio
        vp[2] *= ratio
    end
    if kind in (:wireframe, :surface, :plot3, :scatter3, :trisurf)
        extent = min(vp[2] - vp[1], vp[4] - vp[3])
        vp1 = 0.5 * (vp[1] + vp[2] - extent)
        vp2 = 0.5 * (vp[1] + vp[2] + extent)
        vp3 = 0.5 * (vp[3] + vp[4] - extent)
        vp4 = 0.5 * (vp[3] + vp[4] + extent)
    else
        vp1, vp2, vp3, vp4 = vp
    end
    viewport[1] = vp1 + 0.125 * (vp2 - vp1)
    viewport[2] = vp1 + 0.925 * (vp2 - vp1)
    viewport[3] = vp3 + 0.125 * (vp4 - vp3)
    viewport[4] = vp3 + 0.925 * (vp4 - vp3)

    if kind in (:contour, :contourf, :hexbin, :heatmap, :surface, :trisurf)
        viewport[2] -= 0.1
    end
    GR.setviewport(viewport[1], viewport[2], viewport[3], viewport[4])

    plt.kvs[:viewport] = viewport
    plt.kvs[:vp] = vp
    plt.kvs[:ratio] = ratio

    if haskey(plt.kvs, :backgroundcolor)
        GR.savestate()
        GR.selntran(0)
        GR.setfillintstyle(GR.INTSTYLE_SOLID)
        GR.setfillcolorind(plt.kvs[:backgroundcolor])
        if w > h
          GR.fillrect(subplot[1], subplot[2],
                      ratio * subplot[3], ratio * subplot[4])
        else
          GR.fillrect(ratio * subplot[1], ratio * subplot[2],
                      subplot[3], subplot[4])
        end
        GR.selntran(1)
        GR.restorestate()
    end

    if kind == :polar
        xmin, xmax, ymin, ymax = viewport
        xcenter = 0.5 * (xmin + xmax)
        ycenter = 0.5 * (ymin + ymax)
        r = 0.5 * min(xmax - xmin, ymax - ymin)
        GR.setviewport(xcenter - r, xcenter + r, ycenter - r, ycenter + r)
    end
end

function fix_minmax(a, b)
    if a == b
        a -= a != 0 ? 0.1 * a : 0.1
        b += b != 0 ? 0.1 * b : 0.1
    end
    a, b
end

function given(a)
    a != Nothing && a != "Nothing"
end

function minmax()
    xmin = ymin = zmin = typemax(Float64)
    xmax = ymax = zmax = typemin(Float64)
    for (x, y, z, c, spec) in plt.args
        if given(x)
            xmin = min(minimum(x), xmin)
            xmax = max(maximum(x), xmax)
        else
            xmin, xmax = 0, 1
        end
        if given(y)
            ymin = min(minimum(y), ymin)
            ymax = max(maximum(y), ymax)
        else
            ymin, ymax = 0, 1
        end
        if given(z)
            zmin = min(minimum(z), zmin)
            zmax = max(maximum(z), zmax)
        end
    end
    xmin, xmax = fix_minmax(xmin, xmax)
    ymin, ymax = fix_minmax(ymin, ymax)
    zmin, zmax = fix_minmax(zmin, zmax)
    if haskey(plt.kvs, :xlim)
        plt.kvs[:xrange] = plt.kvs[:xlim]
    else
        plt.kvs[:xrange] = xmin, xmax
    end
    if haskey(plt.kvs, :ylim)
        plt.kvs[:yrange] = plt.kvs[:ylim]
    else
        plt.kvs[:yrange] = ymin, ymax
    end
    if haskey(plt.kvs, :zlim)
        plt.kvs[:zrange] = plt.kvs[:zlim]
    else
        plt.kvs[:zrange] = zmin, zmax
    end
end

function set_window(kind)
    scale = 0
    if kind != :polar
        get(plt.kvs, :xlog, false) && (scale |= GR.OPTION_X_LOG)
        get(plt.kvs, :ylog, false) && (scale |= GR.OPTION_Y_LOG)
        get(plt.kvs, :zlog, false) && (scale |= GR.OPTION_Z_LOG)
        get(plt.kvs, :xflip, false) && (scale |= GR.OPTION_FLIP_X)
        get(plt.kvs, :yflip, false) && (scale |= GR.OPTION_FLIP_Y)
        get(plt.kvs, :zflip, false) && (scale |= GR.OPTION_FLIP_Z)
    end

    minmax()

    if kind in (:wireframe, :surface, :plot3, :scatter3, :polar, :trisurf)
        major_count = 2
    else
        major_count = 5
    end

    xmin, xmax = plt.kvs[:xrange]
    if scale & GR.OPTION_X_LOG == 0
        if !haskey(plt.kvs, :xlim)
            xmin, xmax = GR.adjustlimits(xmin, xmax)
        end
        majorx = major_count
        xtick = GR.tick(xmin, xmax) / majorx
    else
        xtick = majorx = 1
    end
    if scale & GR.OPTION_FLIP_X == 0
        xorg = (xmin, xmax)
    else
        xorg = (xmax, xmin)
    end
    plt.kvs[:xaxis] = xtick, xorg, majorx

    ymin, ymax = plt.kvs[:yrange]
    if kind in (:stem, :hist) && !haskey(plt.kvs, :ylim)
        ymin = 0
    end
    if scale & GR.OPTION_Y_LOG == 0
        if !haskey(plt.kvs, :ylim)
            ymin, ymax = GR.adjustlimits(ymin, ymax)
        end
        majory = major_count
        ytick = GR.tick(ymin, ymax) / majory
    else
        ytick = majory = 1
    end
    if scale & GR.OPTION_FLIP_Y == 0
        yorg = (ymin, ymax)
    else
        yorg = (ymax, ymin)
    end
    plt.kvs[:yaxis] = ytick, yorg, majory

    if kind in (:wireframe, :surface, :plot3, :scatter3, :trisurf)
        zmin, zmax = plt.kvs[:zrange]
        if scale & GR.OPTION_Z_LOG == 0
            if !haskey(plt.kvs, :zlim)
                zmin, zmax = GR.adjustlimits(zmin, zmax)
            end
            majorz = major_count
            ztick = GR.tick(zmin, zmax) / majorz
        else
            ztick = majorz = 1
        end
        if scale & GR.OPTION_FLIP_Z == 0
            zorg = (zmin, zmax)
        else
            zorg = (zmax, zmin)
        end
        plt.kvs[:zaxis] = ztick, zorg, majorz
    end

    plt.kvs[:window] = xmin, xmax, ymin, ymax
    if kind != :polar
        GR.setwindow(xmin, xmax, ymin, ymax)
    else
        GR.setwindow(-1, 1, -1, 1)
    end
    if kind in (:wireframe, :surface, :plot3, :scatter3, :trisurf)
        rotation = get(plt.kvs, :rotation, 40)
        tilt = get(plt.kvs, :tilt, 70)
        GR.setspace(zmin, zmax, rotation, tilt)
    end

    plt.kvs[:scale] = scale
    GR.setscale(scale)
end

function draw_axes(kind, pass=1)
    viewport = plt.kvs[:viewport]
    vp = plt.kvs[:vp]
    ratio = plt.kvs[:ratio]
    xtick, xorg, majorx = plt.kvs[:xaxis]
    ytick, yorg, majory = plt.kvs[:yaxis]

    GR.setlinecolorind(1)
    diag = sqrt((viewport[2] - viewport[1])^2 + (viewport[4] - viewport[3])^2)
    GR.setlinewidth(1)
    charheight = max(0.018 * diag, 0.012)
    GR.setcharheight(charheight)
    ticksize = 0.0075 * diag
    if kind in (:wireframe, :surface, :plot3, :scatter3, :trisurf)
        ztick, zorg, majorz = plt.kvs[:zaxis]
        if pass == 1
            GR.grid3d(xtick, 0, ztick, xorg[1], yorg[2], zorg[1], 2, 0, 2)
            GR.grid3d(0, ytick, 0, xorg[1], yorg[2], zorg[1], 0, 2, 0)
        else
            GR.axes3d(xtick, 0, ztick, xorg[1], yorg[1], zorg[1], majorx, 0, majorz, -ticksize)
            GR.axes3d(0, ytick, 0, xorg[2], yorg[1], zorg[1], 0, majory, 0, ticksize)
        end
    else
        if kind == :heatmap
            ticksize = -ticksize
        else
            GR.grid(xtick, ytick, 0, 0, majorx, majory)
        end
        GR.axes(xtick, ytick, xorg[1], yorg[1], majorx, majory, ticksize)
        GR.axes(xtick, ytick, xorg[2], yorg[2], -majorx, -majory, -ticksize)
    end

    if haskey(plt.kvs, :title)
        GR.savestate()
        GR.settextalign(GR.TEXT_HALIGN_CENTER, GR.TEXT_VALIGN_TOP)
        text(0.5 * (viewport[1] + viewport[2]), vp[4], plt.kvs[:title])
        GR.restorestate()
    end
    if kind in (:wireframe, :surface, :plot3, :scatter3, :trisurf)
        xlabel = get(plt.kvs, :xlabel, "")
        ylabel = get(plt.kvs, :ylabel, "")
        zlabel = get(plt.kvs, :zlabel, "")
        GR.titles3d(xlabel, ylabel, zlabel)
    else
        if haskey(plt.kvs, :xlabel)
            GR.savestate()
            GR.settextalign(GR.TEXT_HALIGN_CENTER, GR.TEXT_VALIGN_BOTTOM)
            text(0.5 * (viewport[1] + viewport[2]), vp[3] + 0.5 * charheight, plt.kvs[:xlabel])
            GR.restorestate()
        end
        if haskey(plt.kvs, :ylabel)
            GR.savestate()
            GR.settextalign(GR.TEXT_HALIGN_CENTER, GR.TEXT_VALIGN_TOP)
            GR.setcharup(-1, 0)
            text(vp[1] + 0.5 * charheight, 0.5 * (viewport[3] + viewport[4]), plt.kvs[:ylabel])
            GR.restorestate()
        end
    end
end

function draw_polar_axes()
    viewport = plt.kvs[:viewport]
    diag = sqrt((viewport[2] - viewport[1])^2 + (viewport[4] - viewport[3])^2)
    charheight = max(0.018 * diag, 0.012)

    window = plt.kvs[:window]
    rmin, rmax = window[3], window[4]

    GR.savestate()
    GR.setcharheight(charheight)
    GR.setlinetype(GR.LINETYPE_SOLID)

    tick = 0.5 * GR.tick(rmin, rmax)
    n = round(Int, (rmax - rmin) / tick + 0.5)
    for i in 0:n
        r = float(i) / n
        if i % 2 == 0
            GR.setlinecolorind(88)
            if i > 0
                GR.drawarc(-r, r, -r, r, 0, 359)
            end
            GR.settextalign(GR.TEXT_HALIGN_LEFT, GR.TEXT_VALIGN_HALF)
            x, y = GR.wctondc(0.05, r)
            GR.text(x, y, string(signif(rmin + i * tick, 12)))
        else
            GR.setlinecolorind(90)
            GR.drawarc(-r, r, -r, r, 0, 359)
        end
    end
    for alpha in 0:45:315
        a = alpha + 90
        sinf = sin(a * π / 180)
        cosf = cos(a * π / 180)
        GR.polyline([sinf, 0], [cosf, 0])
        GR.settextalign(GR.TEXT_HALIGN_CENTER, GR.TEXT_VALIGN_HALF)
        x, y = GR.wctondc(1.1 * sinf, 1.1 * cosf)
        GR.textext(x, y, string(alpha, "^o"))
    end
    GR.restorestate()
end

function inqtext(x, y, s)
    if length(s) >= 2 && s[1] == '$' && s[end] == '$'
        GR.inqmathtex(x, y, s[2:end-1])
    elseif search(s, '\\') != 0 || search(s, '_') != 0 || search(s, '^') != 0
        GR.inqtextext(x, y, s)
    else
        GR.inqtext(x, y, s)
    end
end

function text(x, y, s)
    if length(s) >= 2 && s[1] == '$' && s[end] == '$'
        GR.mathtex(x, y, s[2:end-1])
    elseif search(s, '\\') != 0 || search(s, '_') != 0 || search(s, '^') != 0
        GR.textext(x, y, s)
    else
        GR.text(x, y, s)
    end
end

function draw_legend()
    viewport = plt.kvs[:viewport]
    location = get(plt.kvs, :location, 1)
    num_labels = length(plt.kvs[:labels])
    GR.savestate()
    GR.selntran(0)
    GR.setscale(0)
    w = 0
    h = 0
    for label in plt.kvs[:labels]
        tbx, tby = inqtext(0, 0, label)
        w  = max(w, tbx[3])
        h += max(tby[3] - tby[1], 0.03)
    end
    if location in (8, 9, 10)
        px = 0.5 * (viewport[1] + viewport[2] - w)
    elseif location in (2, 3, 6)
        px = viewport[1] + 0.11
    else
        px = viewport[2] - 0.05 - w
    end
    if location in (5, 6, 7, 10)
        py = 0.5 * (viewport[3] + viewport[4] + h) - 0.03
    elseif location in (3, 4, 8)
        py = viewport[3] + h
    else
        py = viewport[4] - 0.06
    end
    GR.setfillintstyle(GR.INTSTYLE_SOLID)
    GR.setfillcolorind(0)
    GR.fillrect(px - 0.08, px + w + 0.02, py + 0.03, py - h)
    GR.setlinetype(GR.LINETYPE_SOLID)
    GR.setlinecolorind(1)
    GR.setlinewidth(1)
    GR.drawrect(px - 0.08, px + w + 0.02, py + 0.03, py - h)
    i = 1
    GR.uselinespec(" ")
    for (x, y, z, c, spec) in plt.args
        if i <= num_labels
            label = plt.kvs[:labels][i]
            tbx, tby = inqtext(0, 0, label)
            dy = max((tby[3] - tby[1]) - 0.03, 0)
            py -= 0.5 * dy
        end
        GR.savestate()
        mask = GR.uselinespec(spec)
        mask in (0, 1, 3, 4, 5) && GR.polyline([px - 0.07, px - 0.01], [py, py])
        mask & 0x02 != 0 && GR.polymarker([px - 0.06, px - 0.02], [py, py])
        GR.restorestate()
        GR.settextalign(GR.TEXT_HALIGN_LEFT, GR.TEXT_VALIGN_HALF)
        if i <= num_labels
            text(px, py, label)
            py -= 0.5 * dy
            i += 1
        end
        py -= 0.03
    end
    GR.selntran(1)
    GR.restorestate()
end

function colorbar(off=0, colors=256)
    GR.savestate()
    viewport = plt.kvs[:viewport]
    zmin, zmax = plt.kvs[:zrange]
    if get(plt.kvs, :zflip, false)
        options = (GR.inqscale() | GR.OPTION_FLIP_Y) & ~GR.OPTION_FLIP_X
        GR.setscale(options)
    elseif get(plt.kvs, :yflip, false)
        options = GR.inqscale() & ~GR.OPTION_FLIP_Y & ~GR.OPTION_FLIP_X
        GR.setscale(options)
    else
        options = GR.inqscale() & ~GR.OPTION_FLIP_X
        GR.setscale(options)
    end
    GR.setwindow(0, 1, zmin, zmax)
    GR.setviewport(viewport[2] + 0.02 + off, viewport[2] + 0.05 + off,
                   viewport[3], viewport[4])
    l = zeros(Int32, 1, colors)
    l[1,:] = Int[round(Int, _i) for _i in linspace(1000, 1255, colors)]
    GR.cellarray(0, 1, zmax, zmin, 1, colors, l)
    GR.setlinecolorind(1)
    diag = sqrt((viewport[2] - viewport[1])^2 + (viewport[4] - viewport[3])^2)
    charheight = max(0.016 * diag, 0.012)
    GR.setcharheight(charheight)
    if plt.kvs[:scale] & GR.OPTION_Z_LOG == 0
        ztick = 0.5 * GR.tick(zmin, zmax)
        GR.axes(0, ztick, 1, zmin, 0, 1, 0.005)
    else
        GR.setscale(GR.OPTION_Y_LOG)
        GR.axes(0, 2, 1, zmin, 0, 1, 0.005)
    end
    GR.restorestate()
end

function colormap()
    rgb = zeros(256, 3)
    for colorind in 1:256
        color = GR.inqcolor(999 + colorind)
        rgb[colorind, 1] = float( color        & 0xff) / 255.0
        rgb[colorind, 2] = float((color >> 8)  & 0xff) / 255.0
        rgb[colorind, 3] = float((color >> 16) & 0xff) / 255.0
    end
    rgb
end

function to_rgba(value, cmap)
    if !isnan(value)
        r, g, b = cmap[round(Int, value * 255 + 1), :]
        a = 1.0
    else
        r, g, b, a = zeros(4)
    end
    round(UInt32, a * 255) << 24 + round(UInt32, b * 255) << 16 +
    round(UInt32, g * 255) << 8  + round(UInt32, r * 255)
end

function create_context(kind, dict)
    plt.kvs[:kind] = kind
    plt.obj = copy(plt.kvs)
    for (k, v) in dict
        if ! (k in kw_args)
            println("Invalid keyword: $k")
        end
    end
    merge!(plt.kvs, dict)
end

function restore_context()
    global ctx
    ctx = copy(plt.kvs)
    plt.kvs = copy(plt.obj)
end

function figure(; kv...)
    global plt
    plt = Figure()
    merge!(plt.kvs, Dict(kv))
    plt
end

function hold(flag)
    global ctx
    if plt.args != @_tuple(Any)
        plt.kvs[:ax] = flag
        plt.kvs[:clear] = !flag
        for k in (:window, :scale, :xaxis, :yaxis, :zaxis)
            if haskey(ctx, k)
                plt.kvs[k] = ctx[k]
            end
        end
    else
        println("Invalid hold state")
    end
    flag
end

function usecolorscheme(index)
    global scheme
    if 1 <= index <= 4
        scheme = index
    else
        println("Invalid color sheme")
    end
end

function subplot(nr, nc, p)
    xmin, xmax, ymin, ymax = 1, 0, 1, 0
    for i in collect(p)
        r = nr - div(i-1, nc)
        c = (i-1) % nc + 1
        xmin = min(xmin, (c-1)/nc)
        xmax = max(xmax, c/nc)
        ymin = min(ymin, (r-1)/nr)
        ymax = max(ymax, r/nr)
    end
    plt.kvs[:subplot] = [xmin, xmax, ymin, ymax]
    plt.kvs[:clear] = collect(p)[1] == 1
    plt.kvs[:update] = collect(p)[end] == nr * nc
end

function plot_img(I)
    viewport = plt.kvs[:vp][:]
    if haskey(plt.kvs, :title)
        viewport[4] -= 0.05
    end
    vp = plt.kvs[:vp]

    if isa(I, AbstractString)
        width, height, data = GR.readimage(I)
    else
        width, height = size(I)
        data = (float(I) .- minimum(I)) ./ (maximum(I) .- minimum(I))
        data = Int32[round(Int32, 1000 + _i * 255) for _i in data]
    end

    if width  * (viewport[4] - viewport[3]) <
        height * (viewport[2] - viewport[1])
        w = float(width) / height * (viewport[4] - viewport[3])
        xmin = max(0.5 * (viewport[1] + viewport[2] - w), viewport[1])
        xmax = min(0.5 * (viewport[1] + viewport[2] + w), viewport[2])
        ymin = viewport[3]
        ymax = viewport[4]
    else
        h = float(height) / width * (viewport[2] - viewport[1])
        xmin = viewport[1]
        xmax = viewport[2]
        ymin = max(0.5 * (viewport[4] + viewport[3] - h), viewport[3])
        ymax = min(0.5 * (viewport[4] + viewport[3] + h), viewport[4])
    end

    GR.selntran(0)
    GR.setscale(0)
    if get(plt.kvs, :xflip, false)
        tmp = xmax; xmax = xmin; xmin = tmp;
    end
    if get(plt.kvs, :yflip, false)
        tmp = ymax; ymax = ymin; ymin = tmp;
    end
    if isa(I, AbstractString)
        GR.drawimage(xmin, xmax, ymin, ymax, width, height, data)
    else
        GR.cellarray(xmin, xmax, ymin, ymax, width, height, data)
    end

    if haskey(plt.kvs, :title)
        GR.savestate()
        GR.settextalign(GR.TEXT_HALIGN_CENTER, GR.TEXT_VALIGN_TOP)
        text(0.5 * (viewport[1] + viewport[2]), vp[4], plt.kvs[:title])
        GR.restorestate()
    end
    GR.selntran(1)
end

function plot_iso(V)
    viewport = plt.kvs[:viewport]

    if viewport[4] - viewport[3] < viewport[2] - viewport[1]
        width = viewport[4] - viewport[3]
        centerx = 0.5 * (viewport[1] + viewport[2])
        xmin = max(centerx - 0.5 * width, viewport[1])
        xmax = min(centerx + 0.5 * width, viewport[2])
        ymin = viewport[3]
        ymax = viewport[4]
    else
        height = viewport[2] - viewport[1]
        centery = 0.5 * (viewport[3] + viewport[4])
        xmin = viewport[1]
        xmax = viewport[2]
        ymin = max(centery - 0.5 * height, viewport[3])
        ymax = min(centery + 0.5 * height, viewport[4])
    end

    GR.selntran(0)
    values = round.(UInt16, (V .- minimum(V)) ./ (maximum(V) .- minimum(V)) .* (2^16-1))
    nx, ny, nz = size(V)
    isovalue = get(plt.kvs, :isovalue, 0.5)
    rotation = get(plt.kvs, :rotation, 40) * π / 180.0
    tilt = get(plt.kvs, :tilt, 70) * π / 180.0
    r = 2.5
    gr3.clear()
    mesh = gr3.createisosurfacemesh(values, (2/(nx-1), 2/(ny-1), 2/(nz-1)),
                                    (-1., -1., -1.),
                                    round(Int64, isovalue * (2^16-1)))
    if haskey(plt.kvs, :color)
        color = plt.kvs[:color]
    else
        color = (0.0, 0.5, 0.8)
    end
    gr3.setbackgroundcolor(1, 1, 1, 0)
    gr3.drawmesh(mesh, 1, (0, 0, 0), (0, 0, 1), (0, 1, 0), color, (1, 1, 1))
    gr3.cameralookat(r*sin(tilt)*sin(rotation), r*cos(tilt), r*sin(tilt)*cos(rotation), 0, 0, 0, 0, 1, 0)
    gr3.drawimage(xmin, xmax, ymin, ymax, 500, 500, gr3.DRAWABLE_GKS)
    gr3.deletemesh(mesh)
    GR.selntran(1)
end

function plot_polar(θ, ρ)
    window = plt.kvs[:window]
    rmin, rmax = window[3], window[4]
    ρ = (ρ - rmin) / (rmax - rmin)
    n = length(ρ)
    x, y = zeros(n), zeros(n)
    for i in 1:n
        x[i] = ρ[i] * cos(θ[i])
        y[i] = ρ[i] * sin(θ[i])
    end
    GR.polyline(x, y)
end

function RGB(color)
    rgb = zeros(3)
    rgb[1] = float((color >> 16) & 0xff) / 255.0
    rgb[2] = float((color >> 8)  & 0xff) / 255.0
    rgb[3] = float( color        & 0xff) / 255.0
    rgb
end

to_double(a) = Float64[float(el) for el in a]
to_int(a) = Int32[el for el in a]

function send_meta()
    handle = GR.openmeta()
    if handle != C_NULL
        GR.sendmeta(handle, "o(")
        for (k, v) in plt.kvs
            GR.sendmetaref(handle, string(k), 's', string(v))
        end
        for (x, y, z, c, spec) in plt.args
            given(x) && GR.sendmetaref(handle, "x", 'D', to_double(x))
            given(y) && GR.sendmetaref(handle, "y", 'D', to_double(y))
            given(z) && GR.sendmetaref(handle, "z", 'D', to_double(z))
            given(c) && GR.sendmetaref(handle, "c", 'I', to_int(c))
            given(spec) && GR.sendmetaref(handle, "spec", 's', spec)
        end
        GR.sendmeta(handle, ")")
        GR.closemeta(handle)
    end
end

function send_serialized(target)
    handle = connect(target, 8001)
    io = IOBuffer()
    serialize(io, Dict("kvs" => plt.kvs, "args" => plt.args))
    write(handle, io.data)
    close(handle)
end

function plot_data(flag=true)
    global scheme, background

    if plt.args == @_tuple(Any)
        return
    end

    target = GR.displayname()
    if flag && target != None
        if target != "meta"
            send_serialized()
        else
            send_meta()
        end
        return
    end

    kind = get(plt.kvs, :kind, :line)

    plt.kvs[:clear] && GR.clearws()

    if scheme != 0
        for colorind in 1:8
            color = colors[colorind, scheme]
            if colorind == 1
                background = color
            end
            r, g, b = RGB(color)
            GR.setcolorrep(colorind - 1, r, g, b)
            if scheme != 1
                GR.setcolorrep(distinct_cmap[colorind], r, g, b)
            end
        end
        r, g, b = RGB(colors[1, scheme])
        rdiff, gdiff, bdiff = RGB(colors[2, scheme]) - [r, g, b]
        for colorind in 1:12
            f = (colorind - 1) / 11.0
            GR.setcolorrep(92 - colorind, r + f*rdiff, g + f*gdiff, b + f*bdiff)
        end
    end

    set_viewport(kind, plt.kvs[:subplot])
    if !plt.kvs[:ax]
        set_window(kind)
        if kind == :polar
            draw_polar_axes()
        elseif kind != :imshow && kind != :isosurface
            draw_axes(kind)
        end
    end

    if haskey(plt.kvs, :colormap)
        GR.setcolormap(plt.kvs[:colormap])
    else
        GR.setcolormap(GR.COLORMAP_VIRIDIS)
    end

    GR.uselinespec(" ")
    for (x, y, z, c, spec) in plt.args
        GR.savestate()
        if haskey(plt.kvs, :alpha)
            GR.settransparency(plt.kvs[:alpha])
        end
        if kind == :line
            mask = GR.uselinespec(spec)
            mask in (0, 1, 3, 4, 5) && GR.polyline(x, y)
            mask & 0x02 != 0 && GR.polymarker(x, y)
        elseif kind == :scatter
            GR.setmarkertype(GR.MARKERTYPE_SOLID_CIRCLE)
            if given(z) || given(c)
                if given(c)
                    c = (c - minimum(c)) / (maximum(c) - minimum(c))
                    cind = Int[round(Int, 1000 + _i * 255) for _i in c]
                end
                for i in 1:length(x)
                    given(z) && GR.setmarkersize(z[i] / 100.0)
                    given(c) && GR.setmarkercolorind(cind[i])
                    GR.polymarker([x[i]], [y[i]])
                end
            else
                GR.polymarker(x, y)
            end
        elseif kind == :stem
            GR.setlinecolorind(1)
            GR.polyline([plt.kvs[:window][1]; plt.kvs[:window][2]], [0; 0])
            GR.setmarkertype(GR.MARKERTYPE_SOLID_CIRCLE)
            GR.uselinespec(spec)
            for i = 1:length(y)
                GR.polyline([x[i]; x[i]], [0; y[i]])
                GR.polymarker([x[i]], [y[i]])
            end
        elseif kind == :hist
            ymin = plt.kvs[:window][3]
            for i = 1:length(y)
                GR.setfillcolorind(989)
                GR.setfillintstyle(GR.INTSTYLE_SOLID)
                GR.fillrect(x[i], x[i+1], ymin, y[i])
                GR.setfillcolorind(1)
                GR.setfillintstyle(GR.INTSTYLE_HOLLOW)
                GR.fillrect(x[i], x[i+1], ymin, y[i])
            end
        elseif kind == :contour
            zmin, zmax = plt.kvs[:zrange]
            GR.setspace(zmin, zmax, 0, 90)
            h = linspace(zmin, zmax, 20)
            if length(x) == length(y) == length(z)
                x, y, z = GR.gridit(x, y, z, 200, 200)
            end
            GR.contour(x, y, h, z, 1000)
            colorbar(0, 20)
        elseif kind == :contourf
            zmin, zmax = plt.kvs[:zrange]
            GR.setspace(zmin, zmax, 0, 90)
            scale = plt.kvs[:scale]
            GR.setscale(scale)
            if length(x) == length(y) == length(z)
                x, y, z = GR.gridit(x, y, z, 200, 200)
                z = reshape(z, 200, 200)
            end
            GR.surface(x, y, z, GR.OPTION_CELL_ARRAY)
            colorbar()
        elseif kind == :hexbin
            nbins = get(plt.kvs, :nbins, 40)
            cntmax = GR.hexbin(x, y, nbins)
            if cntmax > 0
                plt.kvs[:zrange] = 0, cntmax
                colorbar()
            end
        elseif kind == :heatmap
            w, h = size(z)
            cmap = colormap()
            data = (float(z) .- minimum(z)) ./ (maximum(z) .- minimum(z))
            if get(plt.kvs, :xflip, false)
                data = reverse(data, dims=1)
            end
            if get(plt.kvs, :yflip, false)
                data = reverse(data, dims=2)
            end
            rgba = [to_rgba(value, cmap) for value = data]
            GR.drawimage(0.5, w + 0.5, h + 0.5, 0.5, w, h, rgba)
            colorbar()
        elseif kind == :wireframe
            if length(x) == length(y) == length(z)
                x, y, z = GR.gridit(x, y, z, 50, 50)
            end
            GR.setfillcolorind(0)
            GR.surface(x, y, z, GR.OPTION_FILLED_MESH)
            draw_axes(kind, 2)
        elseif kind == :surface
            if length(x) == length(y) == length(z)
                x, y, z = GR.gridit(x, y, z, 200, 200)
            end
            if get(plt.kvs, :accelerate, true)
                gr3.clear()
                GR.gr3.surface(x, y, z, GR.OPTION_COLORED_MESH)
            else
                GR.surface(x, y, z, GR.OPTION_COLORED_MESH)
            end
            draw_axes(kind, 2)
            colorbar(0.05)
        elseif kind == :plot3
            GR.polyline3d(x, y, z)
            draw_axes(kind, 2)
        elseif kind == :scatter3
            GR.polymarker3d(x, y, z)
            draw_axes(kind, 2)
        elseif kind == :imshow
            plot_img(z)
        elseif kind == :isosurface
            plot_iso(z)
        elseif kind == :polar
            GR.uselinespec(spec)
            plot_polar(x, y)
        elseif kind == :trisurf
            GR.trisurface(x, y, z)
            draw_axes(kind, 2)
            colorbar(0.05)
        elseif kind == :tricont
            zmin, zmax = plt.kvs[:zrange]
            levels = linspace(zmin, zmax, 20)
            GR.tricontour(x, y, z, levels)
        end
        GR.restorestate()
    end

    if kind in (:line, :scatter, :stem) && haskey(plt.kvs, :labels)
        draw_legend()
    end

    if plt.kvs[:update]
        GR.updatews()
        if GR.isinline()
            restore_context()
            return GR.show()
        end
    end

    flag && restore_context()
    return
end

function plot_args(args; fmt=:xys)
    args = Any[args...]
    parsed_args = Any[]

    while length(args) > 0
        local x, y, z, c
        a = popfirst!(args)
        if isa(a, AbstractVecOrMat) || isa(a, Function)
            elt = eltype(a)
            if elt <: Complex
                x = real(a)
                y = imag(a)
                z = Nothing
                c = Nothing
            elseif elt <: Real || isa(a, Function)
                if fmt == :xys
                    if length(args) >= 1 &&
                       (isa(args[1], AbstractVecOrMat) && eltype(args[1]) <: Real || isa(args[1], Function))
                        x = a
                        y = popfirst!(args)
                        z = Nothing
                        c = Nothing
                    else
                        y = a
                        n = isrowvec(y) ? size(y, 2) : size(y, 1)
                        if haskey(plt.kvs, :xlim)
                            xmin, xmax = plt.kvs[:xlim]
                            x = linspace(xmin, xmax, n)
                        else
                            x = linspace(1, n, n)
                        end
                        z = Nothing
                        c = Nothing
                    end
                elseif fmt == :xyac || fmt == :xyzc
                    if length(args) >= 3 &&
                        isa(args[1], AbstractVecOrMat) && eltype(args[1]) <: Real &&
                       (isa(args[2], AbstractVecOrMat) && eltype(args[2]) <: Real || isa(args[2], Function)) &&
                       (isa(args[3], AbstractVecOrMat) && eltype(args[3]) <: Real || isa(args[3], Function))
                        x = a
                        y = popfirst!(args)
                        z = popfirst!(args)
                        c = popfirst!(args)
                    elseif length(args) >= 2 &&
                        isa(args[1], AbstractVecOrMat) && eltype(args[1]) <: Real &&
                       (isa(args[2], AbstractVecOrMat) && eltype(args[2]) <: Real || isa(args[2], Function))
                        x = a
                        y = popfirst!(args)
                        z = popfirst!(args)
                        c = Nothing
                    elseif fmt == :xyac && length(args) >= 1 &&
                       (isa(args[1], AbstractVecOrMat) && eltype(args[1]) <: Real || isa(args[1], Function))
                        x = a
                        y = popfirst!(args)
                        z = Nothing
                        c = Nothing
                    elseif fmt == :xyzc && length(args) == 0
                        z = a
                        nx, ny = size(z)
                        if haskey(plt.kvs, :xlim)
                            xmin, xmax = plt.kvs[:xlim]
                            x = linspace(xmin, xmax, nx)
                        else
                            x = linspace(1, nx, nx)
                        end
                        if haskey(plt.kvs, :ylim)
                            ymin, ymax = plt.kvs[:ylim]
                            y = linspace(ymin, ymax, ny)
                        else
                            y = linspace(1, ny, ny)
                        end
                        c = Nothing
                    end
                end
            else
                error("expected Real or Complex")
            end
        else
            error("expected array or function")
        end
        if isa(y, Function)
            f = y
            y = Float64[f(a) for a in x]
        end
        if isa(z, Function)
            f = z
            z = Float64[f(a,b) for b in y, a in x]
        end
        spec = ""
        if fmt == :xys && length(args) > 0 && isa(args[1], AbstractString)
            spec = popfirst!(args)
        end
        push!(parsed_args, (x, y, z, c, spec))
    end

    pltargs = Any[]

    for arg in parsed_args
        x, y, z, c, spec = arg

        isa(x, UnitRange) && (x = collect(x))
        isa(y, UnitRange) && (y = collect(y))
        isa(z, UnitRange) && (z = collect(z))
        isa(c, UnitRange) && (c = collect(c))

        isvector(x) && (x = vec(x))

        if typeof(y) == Function
            y = [y(a) for a in x]
        else
            isvector(y) && (y = vec(y))
        end
        if given(z)
            if fmt == :xyzc && typeof(z) == Function
                z = [z(a,b) for a in x, b in y]
            else
                isvector(z) && (z = vec(z))
            end
        end
        if given(c)
            isvector(c) && (c = vec(c))
        end

        local xyzc
        if !given(z)
            if isa(x, AbstractVector) && isa(y, AbstractVector)
                xyzc = [ (x, y, z, c) ]
            elseif isa(x, AbstractVector)
                xyzc = length(x) == size(y, 1) ?
                       [ (x, view(y,:,j), z, c) for j = 1:size(y, 2) ] :
                       [ (x, view(y,i,:), z, c) for i = 1:size(y, 1) ]
            elseif isa(y, AbstractVector)
                xyzc = size(x, 1) == length(y) ?
                       [ (view(x,:,j), y, z, c) for j = 1:size(x, 2) ] :
                       [ (view(x,i,:), y, z, c) for i = 1:size(x, 1) ]
            else
                @assert size(x) == size(y)
                xyzc = [ (view(x,:,j), view(y,:,j), z, c) for j = 1:size(y, 2) ]
            end
        elseif isa(x, AbstractVector) && isa(y, AbstractVector) &&
               (isa(z, AbstractVector) || typeof(z) == Array{Float64,2} ||
                typeof(z) == Array{Int32,2} || typeof(z) == Array{Any,2})
            xyzc = [ (x, y, z, c) ]
        else
            xyzc = [ (vec(float(x)), vec(float(y)), vec(float(z)), c) ]
        end
        for (x, y, z, c) in xyzc
            push!(pltargs, (x, y, z, c, spec))
        end
    end

    pltargs
end

function plot(args::PlotArg...; kv...)
    create_context(:line, Dict(kv))

    if plt.kvs[:ax]
        plt.args = append!(plt.args, plot_args(args))
    else
        plt.args = plot_args(args)
    end

    plot_data()
end

function oplot(args::PlotArg...; kv...)
    create_context(:line, Dict(kv))

    plt.args = append!(plt.args, plot_args(args))

    plot_data()
end

function scatter(args...; kv...)
    create_context(:scatter, Dict(kv))

    plt.args = plot_args(args, fmt=:xyac)

    plot_data()
end

function stem(args...; kv...)
    create_context(:stem, Dict(kv))

    plt.args = plot_args(args)

    plot_data()
end

function hist(x, nbins::Integer=0)
    if nbins <= 1
        nbins = round(Int, 3.3 * log10(length(x))) + 1
    end

    xmin, xmax = extrema(x)
    edges = linspace(xmin, xmax, nbins + 1)
    counts = zeros(nbins)
    buckets = Int[max(2, min(searchsortedfirst(edges, xᵢ), length(edges)))-1 for xᵢ in x]
    for b in buckets
        counts[b] += 1
    end
    collect(edges), counts
end

function histogram(x; kv...)
    create_context(:hist, Dict(kv))

    nbins = get(plt.kvs, :nbins, 0)
    x, y = hist(x, nbins)
    plt.args = [(x, y, Nothing, Nothing, "")]

    plot_data()
end

function contour(args...; kv...)
    create_context(:contour, Dict(kv))

    plt.args = plot_args(args, fmt=:xyzc)

    plot_data()
end

function contourf(args...; kv...)
    create_context(:contourf, Dict(kv))

    plt.args = plot_args(args, fmt=:xyzc)

    plot_data()
end

function hexbin(args...; kv...)
    create_context(:hexbin, Dict(kv))

    plt.args = plot_args(args)

    plot_data()
end

function heatmap(D; kv...)
    create_context(:heatmap, Dict(kv))

    if ndims(D) == 2
        z = D'
        width, height = size(z)
        if !haskey(plt.kvs, :xlim) plt.kvs[:xlim] = (0.5, width + 0.5) end
        if !haskey(plt.kvs, :ylim) plt.kvs[:ylim] = (0.5, height + 0.5) end

        plt.args = [(1:width, 1:height, z, Nothing, "")]

        plot_data()
    else
        error("expected 2-D array")
    end
end

function wireframe(args...; kv...)
    create_context(:wireframe, Dict(kv))

    plt.args = plot_args(args, fmt=:xyzc)

    plot_data()
end

function surface(args...; kv...)
    create_context(:surface, Dict(kv))

    plt.args = plot_args(args, fmt=:xyzc)

    plot_data()
end

function plot3(args...; kv...)
    create_context(:plot3, Dict(kv))

    plt.args = plot_args(args, fmt=:xyzc)

    plot_data()
end

function scatter3(args...; kv...)
    create_context(:scatter3, Dict(kv))

    plt.args = plot_args(args, fmt=:xyzc)

    plot_data()
end

function title(s)
    plt.kvs[:title] = s
end

function xlabel(s)
    plt.kvs[:xlabel] = s
end

function ylabel(s)
    plt.kvs[:ylabel] = s
end

function legend(args::AbstractString...; kv...)
    plt.kvs[:labels] = args
end

function xlim(a)
    plt.kvs[:xlim] = a
end

function ylim(a)
    plt.kvs[:ylim] = a
end

function savefig(filename)
    GR.beginprint(filename)
    plot_data(false)
    GR.endprint()
end

function meshgrid(vx, vy)
    m, n = length(vy), length(vx)
    vx = reshape(vx, 1, n)
    vy = reshape(vy, m, 1)
    (repmat(vx, m, 1), repmat(vy, 1, n))
end

function meshgrid(vx, vy, vz)
    m, n, o = length(vy), length(vx), length(vz)
    vx = reshape(vx, 1, n, 1)
    vy = reshape(vy, m, 1, 1)
    vz = reshape(vz, 1, 1, o)
    om = ones(Int, m)
    on = ones(Int, n)
    oo = ones(Int, o)
    (vx[om, :, oo], vy[:, on, oo], vz[om, on, :])
end

function peaks(n=49)
    x = linspace(-2.5, 2.5, n)
    y = x
    x, y = meshgrid(x, y)
    3 * (1 .- x).^2 .* exp.(-(x.^2) .- (y.+1).^2) .- 10*(x/5 .- x.^3 .- y.^5) .* exp.(-x.^2 .- y.^2) .- 1/3 * exp.(-(x.+1).^2 .- y.^2)
end

function imshow(I; kv...)
    create_context(:imshow, Dict(kv))

    plt.args = [(Nothing, Nothing, I, Nothing, "")]

    plot_data()
end

function isosurface(V; kv...)
    create_context(:isosurface, Dict(kv))

    plt.args = [(Nothing, Nothing, V, Nothing, "")]

    plot_data()
end

function cart2sph(x, y, z)
    azimuth = atan2.(y, x)
    elevation = atan2.(z, sqrt.(x.^2 + y.^2))
    r = sqrt.(x.^2 + y.^2 + z.^2)
    azimuth, elevation, r
end

function sph2cart(azimuth, elevation, r)
    x = r .* cos.(elevation) .* cos.(azimuth)
    y = r .* cos.(elevation) .* sin.(azimuth)
    z = r .* sin.(elevation)
    x, y, z
end

function polar(args...; kv...)
    create_context(:polar, Dict(kv))

    plt.args = plot_args(args)

    plot_data()
end

function trisurf(args...; kv...)
    create_context(:trisurf, Dict(kv))

    plt.args = plot_args(args, fmt=:xyzc)

    plot_data()
end

function tricont(args...; kv...)
    create_context(:tricont, Dict(kv))

    plt.args = plot_args(args, fmt=:xyzc)

    plot_data()
end

function mainloop()
    server = listen(8001)
    try
        while true
            sock = accept(server)
            while isopen(sock)
                io = IOBuffer()
                write(io, read(sock))
                seekstart(io)

                obj = deserialize(io)
                merge!(plt.kvs, obj["kvs"])
                plt.args = obj["args"]

                plot_data(false)
            end
        end
    catch
        true
    end
end

end # module


text:
abstract DataFiles
function bread(df :: DataFiles)
  error("No function for abstract type")
end

export bread


text:
# Export StrTab API
export StrTab,
       handle, strtab_lookup

"""
    StrTab

This type encapsulates a string table within an object file, enabling queries
against the string table for symbol names, section names, etc... The list of
available API operations is given below, with methods that subclasses must
implement marked in emphasis:

### Creation
  - *StrTab()*

### Accessors
  - *handle()*
  - *strtab_lookup()*
"""
abstract type StrTab{H <: ObjectHandle} end

"""
    handle(s::StrTab)

Return the `ObjectHandle` this `StrTab` belongs to.
"""
@mustimplement handle(strtab::StrTab)

"""
    strtab_lookup(s::StrTab, index)

Reads a string from the given `StrTab` at `index`.
"""
@mustimplement strtab_lookup(s::StrTab, index)


text:
using LinearAlgebra
using LinearMaps
using ReverseDiff
using ForwardDiff
using IterativeSolvers

# Gets all the necessary derivatives
function get_closures(f)

  function ∇x_f(x, y)
    return ReverseDiff.gradient((x, y) -> f(x, y), (x, y))[1]
  end

  function ∇y_f(x, y)
    return ReverseDiff.gradient((x, y) -> f(x, y), (x, y))[2]
  end

  function ∇y∇y_f(v, x, y)
    function ∇y_f(y, x)
      return ReverseDiff.gradient((y, x) -> f(x, y), (y, x))[1]
    end
    return ForwardDiff.derivative(h -> ∇y_f(y + h * v, x ), zero(eltype(v)))
  end

  function ∇y∇x_f(v, x, y)
    function ∇y_f(x, y)
      return ReverseDiff.gradient((x, y) -> f(x, y), (x, y))[2]
    end
    return ForwardDiff.derivative(h -> ∇y_f(x + h * v, y ), zero(eltype(v)))
  end

  function ∇x∇x_f(v, x, y)
    function ∇x_f(x, y)
      return ReverseDiff.gradient((x, y) -> f(x, y), (x, y))[1]
    end
    return ForwardDiff.derivative(h -> ∇x_f(x + h * v, y ), zero(eltype(v)))
  end

  function ∇x∇y_f(v, x, y)
    function ∇x_f(y, x)
      return ReverseDiff.gradient((y, x) -> f(x, y), (y, x))[2]
    end
    return ForwardDiff.derivative(h -> ∇x_f(y + h * v, x ), zero(eltype(v)))
  end

  return ∇x_f, ∇y_f, ∇x∇x_f, ∇x∇y_f, ∇y∇y_f, ∇y∇x_f
end

# Gets all the necessary derivatives just using forward mode
# autodiff. Slower for large problems, only used for debugging
function get_closures_forward(f)

  function ∇x_f(x, y)
    return ForwardDiff.gradient(x -> f(x, y), x)
  end

  function ∇y_f(x, y)
    return ForwardDiff.gradient(y -> f(x, y), y)
  end

  function ∇y∇y_f(v, x, y)
    function ∇y_f(y)
      return ForwardDiff.gradient(y -> f(x, y), y)
    end

    return ForwardDiff.derivative(h -> ∇y_f(y + h * v), zero(eltype(v)))
  end

  function ∇y∇x_f(v, x, y)
    function ∇y_f(x)
      return ForwardDiff.gradient(y -> f(x, y), y)
    end

    return ForwardDiff.derivative(h -> ∇y_f(x + h * v), zero(eltype(v)))
  end

  function ∇x∇x_f(v, x, y)
    function ∇x_f(x)
      return ForwardDiff.gradient(x -> f(x, y), x)
    end

    return ForwardDiff.derivative(h -> ∇x_f(x + h * v), zero(eltype(v)))
  end

  function ∇x∇y_f(v, x, y)
    function ∇x_f(y)
      return ForwardDiff.gradient(x -> f(x, y), x)
    end

    return ForwardDiff.derivative(h -> ∇x_f(y + h * v), zero(eltype(v)))
  end

  return ∇x_f, ∇y_f, ∇x∇x_f, ∇x∇y_f, ∇y∇y_f, ∇y∇x_f
end

function GDA!(x, y, ∇x_f, ∇y_g, η)
  ∇x_f = ∇x_f(x,y)
  ∇y_g = ∇y_g(x,y)
  x .-= η .* ∇x_f
  y .-= η .* ∇y_g
  return 2
end

function GDA!(x, y, ∇x_f, ∇y_g, ηx, ηy)
  ∇x_f = ∇x_f(x,y)
  ∇y_g = ∇y_g(x,y)
  x .-= ηx * ∇x_f
  y .-= ηy * ∇y_g
  return 2
end

function OGDA!(x, y, Δx, Δy, ∇x_f, ∇y_g, η )
  ∇x_f = ∇x_f(x,y)
  ∇y_g = ∇y_g(x,y)
  x .-= 2 * η * ∇x_f - η .* Δx
  y .-= 2 * η * ∇y_g - η .* Δy
  Δx .= ∇x_f
  Δy .= ∇y_g
  return 2
end

function OGDA!(x, y, Δx, Δy, ∇x_f, ∇y_g, ηx, ηy )
  ∇x_f = ∇x_f(x,y)
  ∇y_g = ∇y_g(x,y)
  x .-= 2 * ηx * ∇x_f - ηx * Δx
  y .-= 2 * ηy * ∇y_g - ηy * Δy
  Δx .= ∇x_f
  Δy .= ∇y_g
  return 2
end

#Only coincides with SGA in the case where f = -g
function SGA!(x, y, ∇x_f, ∇y_g, ∇x∇y_f, ∇y∇x_g, η, γ)
  ∇x_f = ∇x_f(x,y)
  ∇y_g = ∇y_g(x,y)
  Δx = ∇x_f - γ * ∇x∇y_f(∇y_g, x, y)
  Δy = ∇y_g - γ * ∇y∇x_g(∇x_f, x, y)
  x .-= η * Δx
  y .-= η * Δy
  return 4
end

function SGA!(x, y, ∇x_f, ∇y_g, ∇x∇y_f, ∇y∇x_g, ηx, ηy, γ)
  ∇x_f = ∇x_f(x,y)
  ∇y_g = ∇y_g(x,y)
  Δx = ∇x_f - γ * ∇x∇y_f(∇y_g, x, y)
  Δy = ∇y_g - γ * ∇y∇x_g(∇x_f, x, y)

  x .-= ηx * Δx
  y .-= ηy * Δy
  return 4
end

function conOpt!(x, y, ∇x_f, ∇y_g, ∇x∇x_f, ∇x∇y_f, ∇y∇y_g, ∇y∇x_g, η, γ)
  ∇x_f = ∇x_f(x,y)
  ∇y_g = ∇y_g(x,y)
  Δx = ∇x_f - γ * ∇x∇y_f(∇y_g, x, y) + γ * ∇x∇x_f(∇x_f, x, y)
  Δy = ∇y_g - γ * ∇y∇x_g(∇x_f, x, y) + γ * ∇y∇y_g(∇y_g, x, y)

  x .-= η * Δx
  y .-= η * Δy
  return 6
end

function conOpt!(x, y, ∇x_f, ∇y_g, ∇x∇x_f, ∇x∇y_f, ∇y∇y_g, ∇y∇x_g, ηx, ηy, γ)
  ∇x_f = ∇x_f(x,y)
  ∇y_g = ∇y_g(x,y)
  Δx = ∇x_f - γ * ∇x∇y_f(∇y_g, x, y) + γ * ∇x∇x_f(∇x_f, x, y)
  Δy = ∇y_g - γ * ∇y∇x_g(∇x_f, x, y) + γ * ∇y∇y_g(∇y_g, x, y)

  x .-= ηx * Δx
  y .-= ηy * Δy
  return 6
end

function LCGDA!(x, y, ∇x_f, ∇x∇y_f, ∇y_g, ∇y∇x_g, η)
  ∇x_f = ∇x_f(x, y)
  ∇y_g = ∇y_g(x, y)

  Δx = (∇x_f - η * ∇x∇y_f(∇y_g, x, y))
  Δy = (∇y_g - η * ∇y∇x_g(∇x_f, x, y))

  x .-= η * Δx
  y .-= η * Δy
  return 4
end

function CG_CGDA!(x, y, Δx, Δy, ∇x_f, ∇x∇y_f, ∇y_g, ∇y∇x_g, η, ε)
  ∇x_f = ∇x_f(x, y)
  ∇y_g = ∇y_g(x, y)

  #Defining the matvec for the solution of the system of player x
  function MVec_x(inp)
    yTemp = ∇y∇x_g(inp, x, y)
    return inp - η^2 * ∇x∇y_f(yTemp, x, y)
  end
  M_x = LinearMap{eltype(x)}(MVec_x, size(x,1), size(x,1))
  rhs_x = - (∇x_f - η * ∇x∇y_f(∇y_g, x, y))

  # Adjusting the tolerance of the algorithm, only necessary due to
  # limitations of package, hence not tracked.
  tolad_x = norm(rhs_x) / norm(MVec_x(Δx) - rhs_x)

  Δx .= rhs_x
  Δx, xhist = cg!(Δx, M_x, rhs_x; tol = tolad_x * ε, log = true)

  #Defining the matvec for the solution of the system of player y
  function MVec_y(inp)
    yTemp = ∇x∇y_f(inp, x, y)
    return inp - η^2 * ∇y∇x_g(yTemp, x, y)
  end
  M_y = LinearMap{eltype(x)}(MVec_y, size(y,1), size(y,1))
  rhs_y = - (∇y_g - η * ∇y∇x_g(∇x_f, x, y))

  # Adjusting the tolerance of the algorithm, only necessary due to
  # limitations of package, hence not tracked.
  tolad_y = norm(rhs_y) / norm(MVec_y(Δy) - rhs_y)

  Δy .= rhs_y
  Δy, yhist = cg!(Δy, M_y, rhs_y; tol = tolad_y * ε, log = true)

  x .+= η * Δx
  y .+= η * Δy
  return 2 * (xhist.iters  + yhist.iters) + 4
end

function CG_CGDA!(x, y, Δx, Δy, ∇x_f, ∇x∇y_f, ∇y_g, ∇y∇x_g, ηx, ηy, ε)
  ∇x_f = ∇x_f(x, y)
  ∇y_g = ∇y_g(x, y)

  #Defining the matvec for the solution of the system of player x
  function MVec_x(inp)
    yTemp = ηy * ∇y∇x_g(inp, x, y)
    return ηx \ inp - ∇x∇y_f(yTemp, x, y)
  end
  M_x = LinearMap{eltype(x)}(MVec_x, size(x,1), size(x,1))
  rhs_x = - (∇x_f - ∇x∇y_f(ηy * ∇y_g, x, y))

  # Adjusting the tolerance of the algorithm, only necessary due to
  # limitations of package, hence not tracked in computational complexity
  tolad_x = norm(rhs_x) / norm(MVec_x(Δx) - rhs_x)

  Δx .= rhs_x
  Δx, xhist = cg!(Δx, M_x, rhs_x; tol = tolad_x * ε, log = true)

  #Defining the matvec for the solution of the system of player y
  function MVec_y(inp)
    yTemp = ηx * ∇x∇y_f(inp, x, y)
    return ηy \ inp - ∇y∇x_g(yTemp, x, y)
  end
  M_y = LinearMap{eltype(x)}(MVec_y, size(y,1), size(y,1))
  rhs_y = - (∇y_g - ∇y∇x_g(ηx * ∇x_f, x, y))

  # Adjusting the tolerance of the algorithm, only necessary due to
  # limitations of package, hence not tracked.
  tolad_y = norm(rhs_y) / norm(MVec_y(Δy) - rhs_y)

  Δy .= rhs_y
  Δy, yhist = cg!(Δy, M_y, rhs_y; tol = tolad_y * ε, log = true)

  x .+= Δx
  y .+= Δy
  return 2 * (xhist.iters  + yhist.iters) + 4
end

#CGDA only computing the equilibrium for the first player
function CG_CGDA_x!(x, y, Δx, Δy, ∇x_f, ∇x∇y_f, ∇y_g, ∇y∇x_g, η, ε)
  ∇x_f = ∇x_f(x, y)
  ∇y_g = ∇y_g(x, y)

  #Defining the matvec for the solution of the system of player x
  function MVec_x(inp)
    yTemp = ∇y∇x_g(inp, x, y)
    return inp - η^2 * ∇x∇y_f(yTemp, x, y)
  end
  M_x = LinearMap{eltype(x)}(MVec_x, size(x,1), size(x,1))
  rhs_x = - (∇x_f - η * ∇x∇y_f(∇y_g, x, y))

  # Adjusting the tolerance of the algorithm, only necessary due to
  # limitations of package, hence not tracked.
  tolad_x = norm(rhs_x) / norm(MVec_x(Δx) - rhs_x)
  Δx, xhist = cg!(Δx, M_x, rhs_x; tol = tolad_x * ε, log = true)

  Δy .= - (∇y_g + η * ∇y∇x_g(Δx, x, y))

  x .+= η * Δx
  y .+= η * Δy
  return 2 * (xhist.iters ) + 4
end

function CG_CGDA_x!(x, y, Δx, Δy, ∇x_f, ∇x∇y_f, ∇y_g, ∇y∇x_g, ηx, ηy, ε)
  ∇x_f = ∇x_f(x, y)
  ∇y_g = ∇y_g(x, y)

  #Defining the matvec for the solution of the system of player x
  function MVec_x(inp)
    yTemp = ηy * ∇y∇x_g(inp, x, y)
    return ηx \ inp - ∇x∇y_f(yTemp, x, y)
  end
  M_x = LinearMap{eltype(x)}(MVec_x, size(x,1), size(x,1))
  rhs_x = - (∇x_f - ∇x∇y_f(ηy * ∇y_g, x, y))

  # Adjusting the tolerance of the algorithm, only necessary due to
  # limitations of package, hence not tracked.
  tolad_x = norm(rhs_x) / norm(MVec_x(Δx) - rhs_x)
  Δx, xhist = cg!(Δx, M_x, rhs_x; tol = tolad_x * ε, log = true)

  Δy .= - ηy * (∇y_g + ∇y∇x_g(Δx, x, y))

  x .+= Δx
  y .+= Δy
  return 2 * (xhist.iters ) + 4
end

#CGDA only computing the equilibrium for the second player
function CG_CGDA_y!(x, y, Δx, Δy, ∇x_f, ∇x∇y_f, ∇y_g, ∇y∇x_g, η, ε)
  ∇x_f = ∇x_f(x, y)
  ∇y_g = ∇y_g(x, y)

  #Defining the matvec for the solution of the system of player y
  function MVec_y(inp)
    yTemp = ∇x∇y_f(inp, x, y)
    return inp - η^2 * ∇y∇x_g(yTemp, x, y)
  end
  M_y = LinearMap{eltype(x)}(MVec_y, size(y,1), size(y,1))
  rhs_y = - (∇y_g - η * ∇y∇x_g(∇x_f, x, y))

  # Adjusting the tolerance of the algorithm, only necessary due to
  # limitations of package, hence not tracked.
  tolad_y = norm(rhs_y) / norm(MVec_y(Δy) - rhs_y)
  Δy, yhist = cg!(Δy, M_y, rhs_y; tol = tolad_y * ε, log = true)

  Δx .= - (∇x_f + η * ∇x∇y_f(Δy, x, y))

  x .+= η * Δx
  y .+= η * Δy
  return 2 * (yhist.iters) + 4
end

function CG_CGDA_y!(x, y, Δx, Δy, ∇x_f, ∇x∇y_f, ∇y_g, ∇y∇x_g, ηx, ηy, ε)
  ∇x_f = ∇x_f(x, y)
  ∇y_g = ∇y_g(x, y)

  #Defining the matvec for the solution of the system of player y
  function MVec_y(inp)
    yTemp = ηx * ∇x∇y_f(inp, x, y)
    return ηy \ inp - ∇y∇x_g(yTemp, x, y)
  end
  M_y = LinearMap{eltype(x)}(MVec_y, size(y,1), size(y,1))
  rhs_y = - (∇y_g - ∇y∇x_g(ηx * ∇x_f, x, y))

  # Adjusting the tolerance of the algorithm, only necessary due to
  # limitations of package, hence not tracked.
  tolad_y = norm(rhs_y) / norm(MVec_y(Δy) - rhs_y)
  Δy, yhist = cg!(Δy, M_y, rhs_y; tol = tolad_y * ε, log = true)

  Δx .= - ηx * (∇x_f + ∇x∇y_f(Δy, x, y))

  x .+= Δx
  y .+= Δy

  return 2 * (yhist.iters) + 4
end

max_stars_repo_name: f-t-s/CGD


text:
"""
	module BoundaryConditions
"""
module BoundaryConditions

export noBC
export DirichletBC
export NeumannBC
export FloquetBC
export MatchedBC
export AbstractBC
export AbstractRealBC
export AbstractComplexBC
export LocalBC
export NonLocalBC
export BCLocality

abstract type AbstractBC{SIDE} end
abstract type AbstractRealBC{SIDE} <: AbstractBC{SIDE} end
abstract type AbstractComplexBC{SIDE} <: AbstractBC{SIDE} end

struct noBC{SIDE} <: AbstractRealBC{SIDE} end
struct DirichletBC{SIDE} <: AbstractRealBC{SIDE} end
struct NeumannBC{SIDE} <: AbstractRealBC{SIDE} end
struct FloquetBC{SIDE} <: AbstractComplexBC{SIDE} end
struct MatchedBC{SIDE} <: AbstractComplexBC{SIDE}
	in::Array{Int,1}
	out::Array{Int,1}
end
MatchedBC{SIDE}(;in::Array{Int,1}=Int[],out::Array{Int,1}=Int[]) where SIDE = MatchedBC{SIDE}(in,out)

struct LocalBC{SIDE} end
struct NonLocalBC{SIDE} end

BCLocality(::Union{noBC{SIDE},DirichletBC{SIDE},NeumannBC{SIDE}}) where SIDE = LocalBC{SIDE}
BCLocality(::AbstractBC{SIDE}) where SIDE = NonLocalBC{SIDE}

for TBC ∈ (noBC,DirichletBC,NeumannBC,FloquetBC)
	@eval Base.show(io::IO,::$TBC{SIDE}) where SIDE = print(io,$TBC," (SIDE=",SIDE,")")
end
function Base.show(io::IO,mbc::MatchedBC{SIDE}) where SIDE
	print(io,"MatchedBC (SIDE=",SIDE,";")
	print(io," in:", mbc.in)
	print(io," out:", mbc.out)
	print(io,")")
end

Base.conj(bc::AbstractRealBC) = bc
Base.conj(bc::MatchedBC{SIDE}) where SIDE = MatchedBC{SIDE}(in=bc.out,out=bc.in)

"""
	struct noBC{SIDE}
------
	noBC{SIDE}(...) -> bc
"""
noBC

"""
	struct DirichletBC{SIDE}
------
	DirichletBC{SIDE}(depth) -> bc
"""
DirichletBC

"""
	struct NeumannBC{SIDE}
------
	NeumannBC{SIDE}(depth) -> bc
"""
NeumannBC

"""
	struct FloquetBC{SIDE}
------
	FloquetBC{SIDE}(depth) -> bc
"""
FloquetBC

"""
	struct MatchedBC{SIDE}
------
	MatchedBC{SIDE}(depth) -> bc
"""
MatchedBC

end # module