text:
import Pkg
Pkg.rm(["Grids", "GridGames"])
Pkg.develop(["Grids", "GridGames"])

using CUDA

# Make sure no errors or
# warnings with CUDA
CUDA.versioninfo()

using Base.Threads: @spawn
using MuZero
using GridGames
cd(@__DIR__)

# You need at least 5
# threads.
Threads.nthreads()

# User defined state and game implementations
include("game.jl")

# Play a few sample games to make sure the
# game dynamics are correct. Need to copy and
# paste "play" function into REPL, if using VSCode.
#
# play(Tac)

# Initialize MuZero configuration
include("muzero_config.jl")
config = getconfig(Tac)
config.num_unroll_steps = 5
config.num_actors = 1

# Initialize the policy tests
include("inference_tests.jl")
rprv_tests = makerrvptests(Tac)
rprv_tests.period = 1000
pv_tests = makepolicyvaluetests(Tac)
pv_tests.period = 1000
#valuefn_debug = getvaluefn(Tac)

# Configure Learner
config.num_actors
config.num_unroll_steps
learner = Learner(
    config;
    uuid="b234fc40-a07b-4418-9fbc-adb0278dd4e1",
    dir="~/models/MuZeroModels.jl/tictactoe/",
    tests=[pv_tests, rprv_tests], #
    schedule=tacschedule!,#cycleschedule(min_lr=1e-6, max_lr=9e-3),#
    #valuefn_debug=valuefn_debug
)

xyz = randrrvptest(learner)
xyz.state |> display
xyz.state_next |> display

getpolicyvaluetests(learner)[300] |> svg

l = Learner(Tac)
Learner(Tac; uuid="20B4993D-D1B8-41B7-ABB2-9BC8C70EA5F9")
Learner(Tac; dir="~/models/MuZeroModels.jl/tictactoe/")
Learner(
    Tac;
    uuid="8bcecce1-c870-4e34-8ae9-df6433d033de",
    dir="./models/MuZeroModels.jl/tictactoe/"
)

Learner(config)
Learner(config; uuid="20B4993D-D1B8-41B7-ABB2-9BC8C70EA5F9")
Learner(config; dir="~/models/MuZeroModels.jl/tictactoe/")
Learner(
    config;
    uuid="8bcecce1-c870-4e34-8ae9-df6433d033de",
    dir="./models/MuZeroModels.jl/tictactoe/"
)

Learner(Tac, config)
Learner(Tac, config; uuid="20B4993D-D1B8-41B7-ABB2-9BC8C70EA5F9")
Learner(Tac, config; dir="~/models/MuZeroModels.jl/tictactoe/")
Learner(
    Tac,
    config;
    uuid="8bcecce1-c870-4e34-8ae9-df6433d033de",
    dir="./models/MuZeroModels.jl/tictactoe/"
)

getpvtestresult(learner, 300) |> svg
acc, fails, value_mse = runpvtests(learner)
reward_mse, value_mse, acc, fails = runrrvptests(learner)
getpolicy(learner, getstate(Tac([2,8,6]))) |> display
getvalue(learner, getstate(Tac([2,8,6])), 1) |> display

# Let's make sure our policy test cases look reasonable
randpolicyvaluetest(learner) |> svg
randpolicytest(learner) |> svg
randpolicytest(learner) |> svg
randpolicytest(learner) |> svg

showrrvptestresults(learner)

# The predictions from this untrained model should be random
# and mostly incorrect.
showpolicyvaluetestresults(learner)
randpolicyvaluetestresult(learner) |> svg

# Self-play 10 sample games to ensure proper game implementation.
selfplay!(learner; limit=10)

# Show a random self-played game
g = randgame(learner)
g |> display
g |> svg
randtrainingsample(learner, g) |> display

# Get a random batch of training samples
randtrainingbatch(learner) |> display
randtrainingbatch(learner)[1] |> svg

# If game play seems to be working properly, then start indefinite self play on multiple threads.
r1 = @spawn play!(learner)

# Let's make sure threading is working propery and the games are being recorded
length(learner.replay_buffers)
getavggamestep(learner)
getgamestepquality(learner)
get_steps(learner)
r2 = @spawn fit!(learner, steps=1100)

learner.is_playing = true
res3 = @spawn selfplay!(learner)
learner.is_playing = false

updatecallback!(
    learner,
    LearnerCallback(
        "SaveGamePlayStats",
        savegameplaystats!,
        period=1000,
        event=:afterstep
    )
)

learner.is_training = false
updateschedule!(learner, tacschedule!)
updateschedule!(learner, newcycleschedule(learner, min_lr=1e-5, max_lr=3e-5, cycle_steps=10_000))

# Train three steps to make sure there are no errors. This 1st step
# may take a several minutes if your machine has no usable
# Nvidia GPU.
@show @time fit!(learner, steps=3)

# Let's ensure the initial loss was stored, which means the model is training
showlosses(learner)

# View training progress
showresults(learner)

# Train network for a few more steps
fit!(learner, steps=20)

# View training progress. We should see steady decrease
# in the loss and a steady increase in the accuracy.
showresults(learner)

# Let's train indefinitely on a separate thread and periodically check the
# progress.
res = @spawn fit!(learner; steps=1000000)

fit!(learner; steps=3)

# function losses_svg()
#     function Base.show(io::IO, m::MIME"image/svg+xml", g::GridSVG)
#         draw(g)
#         print(io, read(output_file(g), String))
#     end
# end

# # Visualize any game from batch using the encoded state key
# df = get_df(batch)
# state = TicTacToeState(df."State Key"[1])
# state |> show

# # Or, use the following convenience function to get the state
# # @todo Stopped here!
# get_state(batch, 1) |> show

# Learning rate finder
# @todo implement

# Open up "losses.png" and "accuracy.png" in separate windows
# to view training loss and testing accuracy graphs in real-time.

# Alternatively, print the "losses" and/or "accuracy" arrays
# in the REPL to view progress.
losses
accuracy

# If losses is empty and/or not increasing, see status.txt
# to check if there was an error during training.
`cat ./status.txt` |> run

# Also, you can check "fetch(result)" to see if there was an exception.
# If no exception prints (and the call to fetch hangs), then the training is still
# starting, and you can interrupt the call to fetch
fetch(result)

# Export Flux model
MuZero.export(learner)

@show @time selfplay!(learner; limit=3)
@show @time selfplay!(learner; limit=3)

using BenchmarkTools
@btime selfplay!(learner; limit=3)

@btime selfplay!(learner; limit=10)

@spawn selfplay!(learner)
@spawn selfplay!(learner)
@spawn selfplay!(learner)

length(learner.replay_buffer)

@show @time fit!(learner; steps=3)
@spawn fit!(learner)

learner.stats.losses |> length
learner |> gettraindf |> display
learner |> gettraindf |> mean1000 |> display
learner |> gettraindfall |> mean1000 |> display
stoptraining!(learner)
stopplaying!(learner)
clearreplaybuffer!(learner)
length(learner.replay_buffer)

learner |> show_pv_tests
learner |> show_reward_tests

#@learnertest 4854 |> svg
#@learnersample 33

selfplay!(learner; limit=10)

learner |> show_game

# still a fence post problem here
s = learner |> randtrainingsample
learner |> randtrainingsample |> svg

learner |> randtrainingbatch

using Base.Threads: @spawn
for i = 1:num_actors(learner)
    @spawn selfplay!(learner)
end

length(learner.replay_buffer)

@show get_steps(learner), learner.stats.losses, learner.stats.accuracies, learner.times

@time fit!(learner, steps=1)

learner |> gettraindfall |> display

@show learner.stats.losses, learner.stats.accuracies, learner.times

get_steps(learner), get_accuracy(learner)

# STOPPED HERE

tor = "/home/ubuntu/models/MuZeroModels.jl/tictactoe/eb14fc40-a07b-4418-9fbc-adb0278dd4e1/replay_buffers/e4ffabbe-1638-4580-8cc6-b57c8e7c6dbd/buffer.bin"
tor2 = "/home/ubuntu/models/MuZeroModels.jl/tictactoe/eb14fc40-a07b-4418-9fbc-adb0278dd4e1/replay_buffers/d5281c62-455a-4bd4-a145-0ff86597d058/buffer.bin"
reanalyzebuffer(learner, tor2)

# DEBUGGING PV FAILS
getpolicyvaluetestresults(learner)[1935:1935]
fit!(learner, steps=100) ; learner.stats.policy_test_fails

pvdf = getdf(getpolicyvaluetestresults(learner))
pvdf[learner.stats.policy_test_fails, :]

getpolicyvaluetestresults(learner)[1935:1935]
pt1935 = getpolicyvaluetests(learner)[1935]
getpolicy(learner, pt1935.state)
getvalue(learner, pt1935.state, 5)
getvalue(learner, pt1935.state, 2)
getvalue(learner, pt1935.state, 3)
getvalue(learner, pt1935.state, 7)

pt1154 = getpolicyvaluetests(learner)[1154]
pt1154 |> svg
getpolicy(learner, pt2438.state)

pvdf[pvdf."P/F" .== "Fail", :]

using Plots
plot(getlosses(learner))

# LR finder

writable = learner.replay_buffers.d[learner.replay_buffers.writable_uuid]
g1 = writable[end]
g1.net_step
g1
g1 |> svg
g1s5 = get_state(g1, 5)
g1s5 |> typeof

g2 = writable[end-1]

using Profile
@time play_game(learner)

@profile play_game(learner)

open("/tmp/prof.txt", "w") do s
    Profile.print(IOContext(s))
end

learner.play.representation_network[1].weight
learner.play.representation_network[2]
learner.play.representation_network[2].μ
learner.play.representation_network[2].μ
learner.play.representation_network[2].active
learner.train.representation_network[2].active
learner.test.representation_network[2].active


text:
using Test
using MIPVerify

@testset "flatten.jl" begin

    @testset "Flatten" begin
        @testset "default permutation" begin
            p = Flatten(5)
            @test p.n_dim == 5
            @test all(p.perm .== [5, 4, 3, 2, 1])
        end
        @testset "specified permutation, legal" begin
            p = Flatten([3, 1, 2, 4])
            @test p.n_dim == 4
            @test all(p.perm .== [3, 1, 2, 4])
        end
        @testset "specified permutation, illegal" begin
            @test_throws DomainError Flatten([3, 1, 4])
        end
        @testset "Base.show" begin
            p = Flatten([3, 1, 2, 4])
            io = IOBuffer()
            Base.show(io, p)
            @test String(take!(io)) ==
                  "Flatten(): flattens 4 dimensional input, with dimensions permuted according to the order [3, 1, 2, 4]"
        end
    end

    @testset "flatten" begin
        xs = reshape(collect(1:16), (2, 2, 2, 2))
        p1 = Flatten([3, 1, 2, 4])
        @test p1(xs) == [1, 5, 2, 6, 3, 7, 4, 8, 9, 13, 10, 14, 11, 15, 12, 16]
        p2 = Flatten([1, 3, 4, 2])
        @test p2(xs) == [1, 2, 5, 6, 9, 10, 13, 14, 3, 4, 7, 8, 11, 12, 15, 16]
        p3 = Flatten([1, 2, 3, 4])
        @test p3(xs) == [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16]
    end

end

ext: jl


text:
module julia_pkgdev_example

greet() = print("Hello World!")
include("add_fun.jl")
export my_f

end # module

# add include("isbfviub.jl") type scripts
# if you want these accessable when you  use the module, export them


text:
################################################# FILE DESCRIPTION #########################################################

# This typical subgraphing process is as follows:
# 1. If an edge subset is provided, use the index table to map these edges onto their indices in the edge dataframe.
#
# 2. If a vertex subset is provided, first compute a list of all edges involving these vertices, and the indices of these
#    edges in the edge dataframe. Next compute a subset of the index table involving these vertices, and reorder it.

# 3. If both a vertex subset and an edge subset are provided, compute the edge subgraph first, and then the vertex subset
#    on the edge subgraph.
#
# 4. Use vector indexing to get a sub-dataframe containing rows for the input vertices.
#
# 5. Use vector indexing to get a sub-dataframe containing rows for the input edges, using the indices calculate earlier.
#
# 6. If property subsets are provided, extract the specified columns from the dataframes.

################################################# IMPORT/EXPORT ############################################################
export subgraph

################################################# IMPLEMENTATION ###########################################################

""" Compute the maximal subgraph """
subgraph(g::Graph) = copy(g)

"""
Compute a graph containing only the input vertices and the edges
between them.
"""
function subgraph(g::Graph, vs::VertexList)
   sv,erows = subgraph(indxs(g), vs)
   Graph(length(vs), nnz(sv), sv, subgraph(vdata(g), vs), subgraph(edata(g), erows), subgraph(lmap(g), vs))
end

"""
Compute a graph containing only the input properties in the vertex dataframe.
"""
function subgraph(g::Graph, ::Colon, vprops::Vector{Symbol})
   Graph(nv(g), ne(g), copy(indxs(g)), vdata(g)[vprops], copy(edata(g)), copy(lmap(g)))
end

"""
Compute a graph containing only the input vertices and the edges
between them. Preserve only the input properties in the vertex dataframe.
"""
function subgraph(g::Graph, vs::VertexList, vprops::Vector{Symbol})
   sv,erows = subgraph(indxs(g), vs)
   Graph(length(vs), nnz(sv), sv, subgraph(vdata(g), vs, vprops), subgraph(edata(g), erows), subgraph(lmap(g), vs))
end

""" Compute a graph containing only the input edges. """
function subgraph(g::Graph, es::EdgeList)
   sv,erows = subgraph(indxs(g), es)
   Graph(nv(g), nnz(sv), sv, copy(vdata(g)), subgraph(edata(g), erows), copy(lmap(g)))
end

"""
Compute a graph containing only the input properties in the edge dataframe.
"""
function subgraph(g::Graph, ::Colon, ::Colon, eprops::Vector{Symbol})
   Graph(nv(g), ne(g), copy(indxs(g)), copy(vdata(g)), edata(g)[eprops], copy(lmap(g)))
end

"""
Compute a graph containing only the input edges. Preserve only the input
edge properties in the edge dataframe
"""
function subgraph(g::Graph, es::EdgeList, eprops::Vector{Symbol})
   sv,erows = subgraph(indxs(g), es)
   Graph(nv(g), nnz(sv), sv, copy(vdata(g)), subgraph(edata(g), erows, eprops), copy(lmap(g)))
end

""" Compute a graph containing only the input vertices, and the input edges """
function subgraph(g::Graph, vs::VertexList, es::EdgeList)
   sv,erows = subgraph(indxs(g), vs, es)
   Graph(length(vs), nnz(sv), sv, subgraph(vdata(g), vs), subgraph(edata(g), erows), subgraph(lmap(g), vs))
end

""" Compute a subgraph on all subset properties """
function subgraph(
   g::Graph,
   vs::VertexList,
   es::EdgeList,
   vprops::Vector{Symbol},
   eprops::Vector{Symbol}
   )
   sv,erows = subgraph(indxs(g), vs, es)
   Graph(
      length(vs),
      nnz(sv),
      sv,
      subgraph(vdata(g), vs, vprops),
      subgraph(edata(g), erows, eprops),
      subgraph(lmap(g), vs)
   )
end


text:
# Types

"""
Gender
"""
@enum Gender male=1 female=2

function make_gender(g :: T) where T <: Integer
    if g == 1
        return male
    elseif g == 2
        return female
    else
        error("Invalid g: $g");
    end
end

"""
Country
"""
abstract type Country end

struct SourceCountry <: Country
    name :: String
    prefix :: String
end

const Congo = SourceCountry("Congo", "rdc");
const Senegal = SourceCountry("Senegal", "sn");
const Ghana = SourceCountry("Ghana", "gh");
# Congo=1 Ghana=2 Senegal=3 Netherlands=4 UK=5

# function is_source_country(ctry :: Country)
#     return Int(ctry) <= 3
# end

# function is_host_country(ctry :: Country)
#     return !is_source_country(ctry)
# end

"""
Person
"""
mutable struct Person
    id :: String
    sex :: Gender
    birthYear :: Int
    birthCountry :: Country
end

# --------------  Spell History

abstract type Spell end
abstract type SpellHistory end

function n_spells(jh :: SpellHistory)
    return length(jh.spellV)
end

function get_spell(mh :: SpellHistory, idx :: T) where T <: Integer
    return mh.spellV[idx]
end

function start_year(mh :: SpellHistory, idx :: Int)
    return mh.spellV[idx].startYear
end

function end_year(mh :: SpellHistory, idx :: Int)
    return mh.spellV[idx].endYear
end

function spell_years(jh :: SpellHistory)
    ns = n_spells(jh);
    syV = Vector{Int}(undef, ns);
    eyV = Vector{Int}(undef, ns);
    if ns > 0
        for i1 = 1 : ns
            syV[i1] = jh.spellV[i1].startYear;
            eyV[i1] = jh.spellV[i1].endYear;
            if i1 > 1
                if syV[i1] < eyV[i1-1]
                    @warn "Start year before previous end year"
                    println("start years: $syV");
                    println("end years: $eyV");
                    error("invalid")
                end
            end
        end
        @assert all(eyV .>= syV)  "End years must be > start years"
    end

    return syV, eyV
end

"""
	$(SIGNATURES)

Return Migration index that matches a year range.
0 if before first spell. (n+1) if after last spell.
If a spell does not contain the year range: return negative.
"""
function spell_from_year_range(mh :: SpellHistory, year1 :: T, year2 :: T) where T <: Integer
    syV, eyV = spell_years(mh);
    nMig = length(syV);
    if nMig < 1
        return 0
    elseif year1 >= eyV[nMig]
        return nMig + 1
    elseif year2 <= syV[1]
        return 0
    else
        # Interior
        idxV = findall((year1 .>= syV)  .&  (year2 .<= eyV));
        if isempty(idxV)
            return -1
        elseif length(idxV) > 1
            if year2 != year1
                @warn "Spells of more than 1 year should have unique indices"
                println("Years: $year1, $year2");
                println("Indices: $idxV");
                error("Aborted")
            end
            return -1
        else
            return idxV[1]
        end
    end
end

# -------------  Job history

"""
Job
"""
mutable struct Job <: Spell
    startYear :: Int
    endYear :: Int
    ctry :: Int
    income :: Float64
    # currency +++++
end

"""
Job history
"""
mutable struct JobHistory <: SpellHistory
    spellV :: Vector{Job}
end

# --------------------  Migration histories

"""
	$(SIGNATURES)

Migration event.
Country code is the Int code from the Migration file.
"""
mutable struct Migration <: Spell
    startYear :: Int
    endYear :: Int
    ctry :: Int
end

"""
	$(SIGNATURES)

Migration history
"""
mutable struct MigrationHistory <: SpellHistory
    spellV :: Vector{Migration}
end

# function n_migrations(mh :: MigrationHistory)
#     return length(mh.migV)
# end

"""
	$(SIGNATURES)

Person histories: Fixed info and spells with jobs
"""
mutable struct PersonHistory <: SpellHistory
    p :: Person
    spellV :: Vector{Job}
end

# --------


text:
module DiffEqUncertainty

using DiffEqBase, Statistics, Distributions, Quadrature

include("probints.jl")
include("koopman.jl")

export ProbIntsUncertainty,AdaptiveProbIntsUncertainty
export koopman_expectation, montecarlo_expectation

end


text:
# using DifferentialEquations,
using BandedMatrices, DiffEqOperators, FastGaussQuadrature, FileIO, ForwardDiff, Interpolations, JLD2, LinearAlgebra
using ModelConstructors, NLsolve, NLPModelsIpopt, OrderedCollections, OrdinaryDiffEq, Printf, Random, Roots
using SparseArrays, StatsBase, UnPack, VectorizedRoutines.Matlab

using EconPDEs: StateGrid, implicit_timestep, finiteschemesolve
using NLPModels: ADNLSModel, LLSModel , FeasibilityFormNLS

import Base: eltype, getindex
import DiffEqBase: initialize!, solve
import Base.show

function Base.show(io::IO,  ::MIME"text/plain", sg::StateGrid)
    @printf io "%s-element %s\n" string(length(sg)) string(typeof(sg))
end

# src/ directory
include("abstract_NLCT_model.jl")
include("util.jl")

# auxiliary/
include("auxiliary/initialize_grid.jl")
include("auxiliary/investment.jl")

# solve/
include("solve/solve.jl")
include("solve/pseudo_transient_relaxation.jl")
include("solve/differentiate.jl")
include("solve/init_derivatives.jl")
include("solve/util.jl")

# SLM/
include("SLM/slm.jl")
include("SLM/eval.jl")
include("SLM/property.jl")
include("SLM/solve_slm_system.jl")
include("SLM/util.jl")

# models/li2020
include("models/li2020/li2020.jl")
include("models/li2020/augment_variables.jl")
include("models/li2020/eqcond.jl")
include("models/li2020/initialize.jl")
include("models/li2020/subspecs.jl")

# models/brusan
include("models/brusan/brusan.jl")
include("models/brusan/augment_variables.jl")
include("models/brusan/eqcond.jl")
include("models/brusan/initialize.jl")
include("models/brusan/subspecs.jl")


text:
gfn() = Grid(
    6,
    7,
    (x,y)->string(rand(1:9)),
    line_color="red",
    title=Label("Title", color="#666", size=32),
    title_right=Label("Title Right", color="green", size=32),
    header=Label("Header", color="orange", size=32),
    header_right=Label("Header Right")
)
gfn()

g2fn() = Grid(
    6,
    7,
    (x,y)->Cell(
        string(rand(1:9)),
        tl="a",
        tr="b",
        bl="c",
        br="c"
    ),
    line_color="green",
    title=Label("Title2", color="#666", size=32),
    title_right=Label("Title Right2", color="green", size=32),
    header=Label("Header2", color="orange", size=32),
    header_right=Label("Header Right2")
)
g2fn()

# Chess pieces
# https://commons.wikimedia.org/wiki/Category:SVG_chess_pieces
# By Cburnett - Own work, CC BY-SA 3.0, https://commons.wikimedia.org/w/index.php?curid=1499803

rows = [
    [17:32],
    [14:32],
    [13:32],
    [12:32],
    [12:32],
    [12:32],
    [12:32],
    [13:32],
    [14:32],
    [15:30],
    [15:27],
    [15:26],
    [15:24],
    [15:23],
    [14:23],
    [13:23],
    [12:23],
    [11:23],
    [10:23],
    [9:23],
    [8:23],
    [7:23],
    [6:23],
    [5:23],
    [4:23],
    [3:23],
    [2:23],
    [1:23],
    [1:23],
    [1:23],
]

function getgridmat(rows)
    m = falses(30, 32)
    for (row, ranges) in enumerate(rows)
        for range in ranges
            m[row, range] .= true
        end
    end
    m
end

mat = getgridmat(rows)

# 30x32
function gettrackfn(mat)
    function racetrack(row,col)
        if !mat[row,col]
            return nothing
        end
        bg_color = nothing
        if row == size(mat, 1)
            bg_color = "#E79280"
        elseif col == size(mat, 2)
            bg_color = "#92BD8E"
        end
        Cell(
            nothing,
            line_thickness=1,
            bg_color=bg_color
        )
    end
end

g4fn() = Grid(
    size(mat, 1),
    size(mat, 2),
    gettrackfn(mat),
    line_thickness=0,
    title_right=Label("Finish Line", color="#333", size=48),
    header=Label("Starting Line", color="#333", size=48),
    margin_bottom=60,
    margin_left=20,
    margin_right=20,
)
g4fn()

t1fn() = Table(
    3,
    1,
    (x,y)->gfn(),
    margin_top=0,
    margin_bottom=0,
    margin_left=50,
    margin_right=50,
)
t1 = t1fn()

t2fn() = Table(
    1,
    3,
    (x,y)->y == 2 ? t1fn() : TableCell(gfn(), valign=:top),
)
t2 = t2fn()


text:
module InverseSolve

using KrylovMethods
using jInv.Utils
using jInv.Mesh
using jInv.ForwardShare
using jInv.LinearSolvers

import jInv.ForwardShare.ForwardProbType

export getName, AbstractMisfit

abstract type  AbstractMisfit end

export AbstractModel,AbstractModelDerivative,AbstractModelTransform

abstract type AbstractModel end
abstract type AbstractModelDerivative end
abstract type AbstractModelTransform end

using Distributed
using SparseArrays
using Printf
using LinearAlgebra

include("HessianPreconditioners.jl")

include("misfitParam.jl")

"""
mutable struct jInv.InverseSolve.InverseParam

Type storing parameters for Inversion.

Fields:
    Minv::AbstractMesh
    modelfun              - model function (evaluated by main worker), see models.jl
    regularizer           - regularizer(s), see regularizer.jl
    alpha                 - regularization parameter(s)
    mref                  - reference model(s)
    boundsLow::Vector     - lower bounds for model
    boundsHigh::Vector    - upper bounds for model
    maxStep::Real         - maximum step in optimization
    pcgMaxIter::Int       - maximum number of PCG iterations
    pcgTol::Real          - tolerance for PCG
    minUpdate::Real       - stopping criteria
    maxIter::Int          - maximum number of iterations
    HesPrec               - preconditioner for the Hessian.
    metaData::Dict        - Optional dictonary for storeing additional information
Constructor:
    getInverseParam

Example:
    Minv = getRegularMesh(domain,n)
    modelfun = expMod
    regularizer(m,mref,Minv) = wdiffusionReg(m,mref,Minv)
    alpha   = 1e-3
    mref    = zeros(Minv.nc)
    pInv = getInverseParam(Minv,modelfun,regularizer,alpha,mref)
"""
mutable struct  InverseParam
    MInv::AbstractMesh
    modelfun::Function
    regularizer::Union{Function,Array{Function}}
    alpha::Union{Float64,Array{Float64}}
    mref::Array
    boundsLow::Vector
    boundsHigh::Vector
    maxStep::Real
    pcgMaxIter::Int
    pcgTol::Real
    minUpdate::Real
    maxIter::Int
    HesPrec::HessianPreconditioner
    metaData::Dict
end

function Base.display(pInv::InverseParam)
    println("---jInv.InverseSolve.InverseParam---")
    println("inverse mesh type:    $(typeof(pInv.MInv))")
    println("number of cells:      $(pInv.MInv.nc)")
    println("model function:       $(pInv.modelfun)")
    println("maxStep:              $(pInv.maxStep)")
    println("pcgMaxIter:           $(pInv.pcgMaxIter)")
    println("pcgTol:               $(pInv.pcgTol)")
    println("minUpdate:            $(pInv.minUpdate)")
    println("maxIter:              $(pInv.maxIter)")
end

"""
function jInv.InverseSolve.getInverseParam(...)

Constructs an InverseParam

Required Input:

    Minv::AbstractMesh    - mesh of model
    modFun::Function      - model
    regularizer::Function - regularizer, see regularizer.jl
    alpha::Real           - regularization parameter
    mref                  - reference model
    boundsLow::Vector     - lower bounds for model
    boundsHigh::Vector    - upper bounds for model

Optional Inputs:

    maxStep::Real=1.0     - maximum step in optimization
    pcgMaxIter::Int=10    - maximum number of PCG iterations
    pcgTol::Real          - tolerance for PCG
    minUpdate::Real=1e-4  - stopping criteria
    maxIter::Int=10       - maximum number of iterations
"""
function getInverseParam(MInv::AbstractMesh,
                         modFun::Function,
                         regularizer::Union{Function,Array{Function}},
                         alpha::Union{Float64,Array{Float64}},
                         mref::Array,
                         boundsLow::Vector,
                         boundsHigh::Vector;
                         maxStep::Real=1.0,
                         pcgMaxIter::Int=10,
                         pcgTol::Real=1e-1,
                         minUpdate::Real=1e-4,
                         maxIter::Int=10,
                         HesPrec::HessianPreconditioner=getSSORRegularizationPreconditioner(1.0,1e-15,10),
                         metaData::Dict=Dict())

    return InverseParam(MInv, modFun, regularizer, alpha, mref, boundsLow, boundsHigh, maxStep,
                        pcgMaxIter, pcgTol, minUpdate, maxIter, HesPrec, metaData)
end

export InverseParam, getInverseParam

include("globalToLocal.jl")
include("models.jl")
include("misfit.jl")
include("regularizers.jl")
include("GNsolve.jl")
include("GNhis.jl")
include("projGN.jl")
include("projSD.jl")
include("projPCG.jl")
include("barrierGNCG.jl")
include("computeMisfit.jl")
include("computeGradMisfit.jl")
include("HessMatVec.jl")
include("iteratedTikhonov.jl")
include("getHessian.jl")
end


text:
#!/usr/bin/julia
# memscan.jl PID HEAD TAIL

pid = parse(Int64, ARGS[1])
head = parse(Int64, ARGS[2])
tail = parse(Int64, ARGS[3])
needle = Int64(233)

# open memory mapping
# cat /proc/pic/maps for mem info
mem = open("/proc/$(pid)/mem", "r")

for offset in head:0x8:tail-0x8
	seek(mem, offset)
	num = read(mem, Int64)
	if num == needle
		println(hex(offset), " ", num)
	else
		println(hex(offset), " ", num)
	end
end

max_stars_repo_name: cdluminate/MyNotes

max_issues_repo_name: cdluminate/MyNotes

max_forks_repo_name: cdluminate/MyNotes


text:
using PyPlot, Glob, JLD2

defaultcolors = plt.rcParams["axes.prop_cycle"].by_key()["color"]
removespine(side) = gca().spines[side].set_visible(false)
removespines(sides...) = [removespine(side) for side in sides]

include("ConvergenceTests/ConvergenceTests.jl")

filenames = glob("data/taylor_green*.jld2")

errors = ConvergenceTests.compute_errors(
            (x, y, z, t) -> ConvergenceTests.DoublyPeriodicTaylorGreen.u(x, y, t),
            filenames...)

sizes = ConvergenceTests.extract_sizes(filenames...)

Nx = map(sz -> sz[1], sizes)
L₁ = map(err -> err.L₁, errors)
L∞ = map(err -> err.L∞, errors)

close("all")
fig, ax = subplots()

ax.loglog(Nx, L₁, linestyle="None", marker="o", label="error, \$L_1\$-norm")
ax.loglog(Nx, L∞, linestyle="None", marker="^", label="error, \$L_\\infty\$-norm")
ax.loglog(Nx, L₁[end] * (Nx[end] ./ Nx).^2, "k-", linewidth=1, alpha=0.6,  label=L"\sim N_x^{-2}")

legend()

title("Convergence for freely-decaying Taylor-Green vortex")
removespines("top", "right")
ylabel("Norms of the absolute error, \$ | u_\\mathrm{simulation} - u_\\mathrm{analytical} | \$")
xlabel(L"N_x")
xticks(sort(Nx), ["\$ 2^{$(round(Int, log2(n)))} \$" for n in sort(Nx)])

savefig("figs/taylor_green_convergence.png", dpi=480)


text:
# Script for creating Figure 13 in scale manuscript

using CSV
using DataFrames
using IcanProj
using PyPlot

# Settings

path_results = joinpath(dirname(pathof(IcanProj)), "..", "data")

path_figure = joinpath(dirname(pathof(IcanProj)), "..", "plots", "correlation_plot")

# Load results

file_swe = joinpath(path_results, "correlation_error_swe.csv")

df_swe = CSV.read(file_swe, delim = ";") |> DataFrame

file_latmo = joinpath(path_results, "correlation_error_latmo.csv")

df_latmo = CSV.read(file_latmo, delim = ";") |> DataFrame

file_hatmo = joinpath(path_results, "correlation_error_hatmo.csv")

df_hatmo = CSV.read(file_hatmo, delim = ";") |> DataFrame

file_rnet = joinpath(path_results, "correlation_error_rnet.csv")

df_rnet = CSV.read(file_rnet, delim = ";") |> DataFrame

# Dataframe for bias and rmse

df_rmse = DataFrame(swe = df_swe.r_swe_elev_for_rmse,
                    latmo = df_latmo.r_latmo_elev_for_rmse,
                    hatmo = df_hatmo.r_hatmo_elev_for_rmse,
                    rnet = df_rnet.r_rnet_elev_for_rmse)

df_bias = DataFrame(swe = df_swe.r_swe_elev_for_bias,
                    latmo = df_latmo.r_latmo_elev_for_bias,
                    hatmo = df_hatmo.r_hatmo_elev_for_bias,
                    rnet = df_rnet.r_rnet_elev_for_bias)

# Plot results

fig = figure(figsize = (5, 2.7))

boxplot(convert(Array{Float64}, df_rmse).^2, 0, "")
ylim([0, 1])
yticks(0:0.2:1)
ylabel("Squared correlation\ncoefficient (-)")
xticks(collect(1:4), uppercase.(String.(names(df_rmse))))

savefig(joinpath(path_figure, "correlation_plot.pdf"), dpi = 600)

# Print summary statistics

describe(df_rmse)


text:
import Base: exp
export exp!

# Integrating factor

intfact(x, a) = exp(-2a)besseli(x,2a)

intfact(x, y,a) = exp(-4a)besseli(x,2a)besseli(y,2a)

intfact(x, y, z, a) = exp(-6a)besseli(x,2a)besseli(y,2a)besseli(z,2a)

"""
    plan_intfact(a::Real,dims::Tuple,[fftw_flags=FFTW.ESTIMATE])

Constructor to set up an operator for evaluating the integrating factor with
real-valued parameter `a`. This can then be applied with the `*` or `\` operation on
data of the appropriate size.

Note that `a` can be positive, negative, or zero. However, if `a` is negative,
then only the `\` operation is actually correct; the `*` operation merely
returns the identity to avoid excessive (and noisy) calculation. Similarly, if
`a` is positive, the `\` operation returns the identity. Thus, these operations are
not inverses of one another. If `a` is zero, both operations return the identity.

The `dims` argument can be replaced with data of type `ScalarGridData` to
specify the size of the domain.

# Example

```jldoctest
julia> w = Nodes(Dual,(6,6));

julia> w[4,4] = 1.0;

julia> E = plan_intfact(1.0,(6,6))
Integrating factor with parameter 1.0 on a (nx = 6, ny = 6) grid

julia> E*w
Nodes{Dual,6,6,Float64} data
Printing in grid orientation (lower left is (1,1))
6×6 Array{Float64,2}:
 0.00268447   0.00869352  0.0200715   0.028765    0.0200715   0.00869352
 0.00619787   0.0200715   0.0463409   0.0664124   0.0463409   0.0200715
 0.00888233   0.028765    0.0664124   0.0951774   0.0664124   0.028765
 0.00619787   0.0200715   0.0463409   0.0664124   0.0463409   0.0200715
 0.00268447   0.00869352  0.0200715   0.028765    0.0200715   0.00869352
 0.000828935  0.00268447  0.00619787  0.00888233  0.00619787  0.00268447
```
"""
function plan_intfact end

"""
    plan_intfact!(a::Real,dims::Tuple,[fftw_flags=FFTW.ESTIMATE][,nthreads=length(Sys.cpu_info())])

Same as [`plan_intfact`](@ref), but the resulting operator performs an in-place
operation on data. The number of threads `threads` defaults to the number of
logical CPU cores on the system.
"""
function plan_intfact! end

abstract type IFExpSign end
abstract type PosExp <: IFExpSign end
abstract type NegExp <: IFExpSign end
abstract type ZeroExp <: IFExpSign end
NonNegExp = Union{PosExp,ZeroExp}
NonPosExp = Union{NegExp,ZeroExp}

struct IntFact{NX, NY, PA<:IFExpSign, inplace}
    a::Float64
    conv::Union{CircularConvolution{NX, NY},Nothing}
end

for (lf,inplace) in ((:plan_intfact,false),
                     (:plan_intfact!,true))

    @eval function $lf(a::Real,dims::Tuple{Int,Int};fftw_flags = FFTW.ESTIMATE, nthreads = length(Sys.cpu_info()))
        NX, NY = dims

        if a == 0
          return IntFact{NX, NY, ZeroExp, $inplace}(0.0,nothing)
        elseif a < 0
          a_internal = abs(convert(Float64,a))
          signType = NegExp
        else
          a_internal = convert(Float64,a)
          signType = PosExp
        end

        #qtab = [intfact(x, y, a) for x in 0:NX-1, y in 0:NY-1]
        Nmax = 0
        while abs(intfact(Nmax,0,a_internal)) > eps(Float64)
          Nmax += 1
        end
        qtab = [max(x,y) <= Nmax ? intfact(x, y, a_internal) : 0.0 for x in 0:NX-1, y in 0:NY-1]
        #IntFact{NX, NY, a, $inplace}(Nullable(CircularConvolution(qtab, fftw_flags)))
        IntFact{NX, NY, signType, $inplace}(a_internal,CircularConvolution(qtab, fftw_flags,nthreads=nthreads))
      end

      # Base the size on the dual grid associated with any grid data, since this
      # is what the efficient grid size in PhysicalGrid has been established with
      @eval $lf(a::Real,::GridData{NX,NY}; fftw_flags = FFTW.ESTIMATE, nthreads = length(Sys.cpu_info())) where {T<:CellType,NX,NY} =
          $lf(a,(NX,NY), fftw_flags = fftw_flags, nthreads = nthreads)

end

function Base.show(io::IO, E::IntFact{NX, NY, signType, inplace}) where {NX, NY, signType, inplace}
    nodedims = "(nx = $NX, ny = $NY)"
    atxt = signType == NegExp ? "$(-E.a)" : "$(E.a)"
    isinplace = inplace ? "In-place integrating factor" : "Integrating factor"
    print(io, "$isinplace with parameter $atxt on a $nodedims grid")
end

"""
    exp(L::Laplacian,a[,Nodes(Dual)])

Create the integrating factor exp(L*a). The default size of the operator is
the one appropriate for dual nodes; another size can be specified by supplying
grid data in the optional third argument. Note that, if `L` contains a factor,
it scales the exponent with this factor.
"""
exp(L::Laplacian{NX,NY},a,prototype=Nodes(Dual,(NX,NY))) where {NX,NY} = plan_intfact(L.factor*a,prototype)

"""
    exp!(L::Laplacian,a[,Nodes(Dual)])

Create the in-place version of the integrating factor exp(L*a).
"""
exp!(L::Laplacian{NX,NY},a,prototype=Nodes(Dual,(NX,NY))) where {NX,NY} = plan_intfact!(L.factor*a,prototype)

for (datatype) in (:Nodes, :XEdges, :YEdges)
  @eval function mul!(out::$datatype{T,NX, NY},
                     E::IntFact{MX, MY, PosExp, inplace},
                     s::$datatype{T, NX, NY}) where {T <: CellType, NX, NY, MX, MY, inplace}

      mul!(out.data, E.conv, s.data)
      out
  end

  @eval function mul!(out::$datatype{T,NX, NY},
                     E::IntFact{MX, MY, NegExp, inplace},
                     s::$datatype{T, NX, NY}) where {T <: CellType, NX, NY, MX, MY, inplace}

      out .= deepcopy(s)
  end

  @eval function ldiv!(out::$datatype{T,NX, NY},
                     E::IntFact{MX, MY, PosExp, inplace},
                     s::$datatype{T, NX, NY}) where {T <: CellType, NX, NY, MX, MY, inplace}

      out .= deepcopy(s)
  end

  @eval function ldiv!(out::$datatype{T,NX, NY},
                     E::IntFact{MX, MY, NegExp, inplace},
                     s::$datatype{T, NX, NY}) where {T <: CellType, NX, NY, MX, MY, inplace}

      mul!(out.data, E.conv, s.data)
      out
  end

  @eval function mul!(out::$datatype{T,NX, NY},
                     E::IntFact{MX, MY, ZeroExp, inplace},
                     s::$datatype{T, NX, NY}) where {T <: CellType, NX, NY, MX, MY, inplace}
      out .= deepcopy(s)
  end

  @eval function ldiv!(out::$datatype{T,NX, NY},
                     E::IntFact{MX, MY, ZeroExp, inplace},
                     s::$datatype{T, NX, NY}) where {T <: CellType, NX, NY, MX, MY, inplace}
      out .= deepcopy(s)
  end

end

for (op) in (:mul!,:ldiv!)
  @eval function $op(out::Edges{C,NX,NY},E::IntFact,s::Edges{C,NX,NY}) where {C,NX,NY}
    $op(out.u,E,s.u)
    $op(out.v,E,s.v)
    out
  end

  @eval function $op(out::EdgeGradient{C,D,NX,NY},E::IntFact,s::EdgeGradient{C,D,NX,NY}) where {C,D,NX,NY}
    $op(out.dudx,E,s.dudx)
    $op(out.dvdx,E,s.dvdx)
    $op(out.dudy,E,s.dudy)
    $op(out.dvdy,E,s.dvdy)
    out
  end

end

#=
function mul!(out::Nodes{T,NX, NY},
                   E::IntFact{MX, MY, a, inplace},
                   s::Nodes{T, NX, NY}) where {T <: CellType, NX, NY, MX, MY, a, inplace}

    mul!(out.data, E.conv, s.data)
    out
end

function mul!(out::Nodes{T,NX, NY},
                   E::IntFact{MX, MY, 0.0, inplace},
                   s::Nodes{T, NX, NY}) where {T <: CellType, NX, NY, MX, MY, inplace}
    out .= deepcopy(s)
end
=#

*(E::IntFact{MX,MY,signType,false},s::G) where {MX,MY,signType,G<:GridData} =
  mul!(G(), E, s)

*(E::IntFact{MX,MY,signType,true},s::GridData) where {MX,MY,signType} =
    mul!(s, E, deepcopy(s))

\(E::IntFact{MX,MY,signType,false},s::G) where {MX,MY,signType,G<:GridData} =
  ldiv!(G(), E, s)

\(E::IntFact{MX,MY,signType,true},s::GridData) where {MX,MY,signType} =
    ldiv!(s, E, deepcopy(s))

max_stars_repo_name: JuliaIBPM/CartesianGrids

max_issues_repo_name: JuliaIBPM/CartesianGrids

max_forks_repo_name: JuliaIBPM/CartesianGrids


text:
using CapacityExpansion
using Clp
exact_scenarios=Dict{String,OptResult}()
for (state, years) in [["GER_1", [2016]],["CA_1", [2016]]]
       # laod data
       ts_input_data = load_timeseries_data_provided(state; T=24, years=years) #CEP
       cep_data = load_cep_data_provided(state)
       ## CLUSTERING ##
       ts_clust_data = run_clust(ts_input_data;method="hierarchical",representation="centroid",n_init=1,n_clust=3)
       ts_full_data = run_clust(ts_input_data;method="hierarchical",representation="centroid",n_init=1,n_clust=30)
       ## OPTIMIZATION ##
       optimizer=Clp.Optimizer
       exact_scenarios["$state-$years-co2"]=run_opt(ts_clust_data.best_results,cep_data,optimizer;descriptor="co2",co2_limit=1000)
       exact_scenarios["$state-$years-slack"]=run_opt(ts_clust_data.best_results,cep_data,optimizer;descriptor="slack",lost_el_load_cost=1e6, lost_CO2_emission_cost=700)
       exact_scenarios["$state-$years-ex"]=run_opt(ts_clust_data.best_results,cep_data,optimizer;descriptor="ex",existing_infrastructure=true)
       exact_scenarios["$state-$years-simple"]=run_opt(ts_clust_data.best_results,cep_data,optimizer;descriptor="simple storage",storage="simple")
       exact_scenarios["$state-$years-seasonal"]=run_opt(ts_clust_data.best_results,cep_data,optimizer;descriptor="seasonal storage",storage="seasonal")
       design_result=run_opt(ts_clust_data.best_results,cep_data,optimizer;descriptor="des&op")
       exact_scenarios["$state-$years-des&op"]=run_opt(ts_full_data.best_results,cep_data,design_result.opt_config,get_cep_design_variables(design_result),optimizer;lost_el_load_cost=1e6,lost_CO2_emission_cost=700)
end
for (state, years) in [["GER_18", [2016]],["CA_14", [2016]]]
       # laod data
       ts_input_data = load_timeseries_data_provided(state; T=24, years=years) #CEP
       cep_data = load_cep_data_provided(state)
       ## CLUSTERING ##
       ts_clust_data = run_clust(ts_input_data;method="hierarchical",representation="centroid",n_init=1,n_clust=3)
       ## OPTIMIZATION ##
       optimizer=Clp.Optimizer
       exact_scenarios["$state-$years-trans"]=run_opt(ts_clust_data.best_results,cep_data,optimizer;descriptor="trans",transmission=true)
end

@save normpath(joinpath(dirname(@__FILE__),"cep_exact_data.jld2")) exact_scenarios


text:
function matchregret(I::Int64, r::Vector{Vector{Float64}}, gs::GameSet)
    ni_stage, na_stage = gs.ni_stage, gs.na_stage
    na = numactions(I, ni_stage, na_stage)
    σ = Vector{Float64}(undef, na)
    for a in 1:na
        denom = sum(max(r[I][b], 0.0) for b in 1:na)
        if denom > 0.0
            σ[a] = max(r[I][a], 0.0) / denom
        else
            σ[a] = 1.0 / na
        end
    end
    return σ
end

function updateutility!(player::Int64, I::Int64, h::SVector, σ::Vector{Float64},
        u::Float64, πi::Float64, πo::Float64,
        r::Vector{Vector{Float64}}, s::Vector{Vector{Float64}},
        g::AirDefenseGame, gs::GameSet, depth::Int64)
    A, ni_stage, na_stage = gs.A, gs.ni_stage, gs.na_stage
    na = numactions(I, ni_stage, na_stage)
    m = Vector{Float64}(undef, na)
    nextactions = actions(depth, A)
    # @show na[I]
    for a in 1:na
        ha = setindex(h, nextactions[a], depth)
        if getplayer(depth) == player
            πip = σ[a] * πi
            # up = updatetree!(ha, player, πip, πo, r, s, g, gs)
            up = updatetree!(ha, player, πip, πo, r, s, g, gs, depth + 1)
            m[a] = up
            u = u + σ[a] * up
        else
            πop = σ[a] * πo
            # up = updatetree!(ha, player, πi, πop, r, s, g, gs)
            up = updatetree!(ha, player, πi, πop, r, s, g, gs, depth + 1)
            u = u + up
        end
    end
    return m, u
end

function updateregret!(player::Int64, I::Int64, h::SVector, σ::Vector{Float64},
        u::Float64, m::Vector{Float64}, πi::Float64,
        r::Vector{Vector{Float64}}, s::Vector{Vector{Float64}}, gs::GameSet, depth ::Int64)
    ni_stage, na_stage = gs.ni_stage, gs.na_stage
    if getplayer(depth) == player
        for a in 1:numactions(I, ni_stage, na_stage)
            r[I][a] = r[I][a] + m[a] - u
            s[I][a] = s[I][a] + πi * σ[a]
        end
    end
end

function updatetree!(h::SVector, player::Int64, πi::Float64, πo::Float64,
    r::Vector{Vector{Float64}}, s::Vector{Vector{Float64}}, g::AirDefenseGame,
    gs::GameSet, depth::Int64)
    # πo is short for \vec{π}_{-i}
    A = gs.A
    (terminal(depth)) && (return leafutility(h, player, πi, πo, g, gs))
    if getplayer(depth) == 3
        as = actions(depth, A)
        chance_weights = pweights([chance(depth, a, g, gs) for a in as])
        a = as[sample(chance_weights)]
        ha = setindex(h, a, depth)
        # return updatetree!(ha, player, πi, πo, r, s, g, gs)
        return updatetree!(ha, player, πi, πo, r, s, g, gs, depth + 1)
    end
    I = infoset(h, depth, g, gs)
    u = 0.0
    σ = matchregret(I, r, gs)
    m, u = updateutility!(player, I, h, σ, u, πi, πo, r, s, g, gs, depth)
    updateregret!(player, I, h, σ, u, m, πi, r, s, gs, depth)
    return u
end

function cfr(T::Int64, g::AirDefenseGame, gs::GameSet)
    ni, ni_stage, na_stage, players = gs.ni, gs.ni_stage, gs.na_stage, gs.players
    r = [zeros(numactions(i, ni_stage, na_stage)) for i in 1:ni] # regret
    s = [zeros(numactions(i, ni_stage, na_stage)) for i in 1:ni] # cumulative strategies
    u = 0.0
    u1 = Vector{Float64}(undef, T)
    println("MCCFR: starting $T iterations...")
    # @progress for t in 1:T
    for t in 1:T
        for player in players
            h0 = SVector(0, 0, 0, 0, 0, 0)
            # u = updatetree!(h0, player, 1.0, 1.0, r, s, g, gs)
            u = updatetree!(h0, player, 1.0, 1.0, r, s, g, gs, 1)
            (player == 1) && (u1[t] = u)
        end
        t % 10 == 0 && print("\rIteration $t complete. Utility: $u")
    end
    println("...$T iterations complete.")
    denoms = [sum(i) for i in s]
    σ = Vector{Vector{Float64}}(undef, gs.ni)
    for i in eachindex(s)
        if denoms[i] == 0
            σ[i] = fill(1 / length(s[i][:]), length(s[i][:]))
        else
            σ[i] = s[i][:] ./ denoms[i]
        end
    end
    return u1, r, s, σ
end

function strat_cfr(s::Vector{Vector{Float64}})
    denoms = [sum(i) for i in s]
    σ = Vector{Vector{Float64}}(undef, length(s))
    for i in eachindex(s)
        if denoms[i] == 0
            σ[i] = fill(1 / length(s[i][:]), length(s[i][:]))
        else
            σ[i] = s[i][:] ./ denoms[i]
        end
    end
    return σ
end

function cfr_exploit(T::Int64, g::AirDefenseGame, gs::GameSet, gos::GameOptSet)
    ni, ni_stage, na_stage, players = gs.ni, gs.ni_stage, gs.na_stage, gs.players
    r = [zeros(numactions(i, ni_stage, na_stage)) for i in 1:ni] # regret
    s = [zeros(numactions(i, ni_stage, na_stage)) for i in 1:ni] # cumulative strategies
    u = 0.0
    Tx = T ÷ 10 # iterations for best response calculations
    tx = 1 # index for best respoinse calculations
    u1 = Vector{Float64}(undef, T)
    u_br = Vector{Float64}(undef, 10)
    println("MCCFR: starting $T iterations...")
    # @progress for t in 1:T
    for t in 1:T
        for player in players
            h0 = SVector(0, 0, 0, 0, 0, 0)
            # u = updatetree!(h0, player, 1.0, 1.0, r, s, g, gs)
            u = updatetree!(h0, player, 1.0, 1.0, r, s, g, gs, 1)
            (player == 1) && (u1[t] = u)
        end
        if t % Tx == 0 # calculate best response
            println("")
            σ = strat_cfr(s)
            rc1, _ = real_strat(σ, gs, gos)
            status, u_br[tx], _ = lp_best_response(gos, rc1, fixedplayer = 1)
            tx += 1
            status != :Optimal && @warn "Best response LP not optimal in iter $t."
        end
        t % 10 == 0 && print("\rIteration $t complete. Utility: $u")
    end
    σ = strat_cfr(s)
    println("...$T iterations complete.")
    return u1, r, s, σ, -u_br
end

function cfr_timebound(T::Int64, g::AirDefenseGame, gs::GameSet)
    ni, ni_stage, na_stage, players = gs.ni, gs.ni_stage, gs.na_stage, gs.players
    r = [zeros(numactions(i, ni_stage, na_stage)) for i in 1:ni] # regret
    s = [zeros(numactions(i, ni_stage, na_stage)) for i in 1:ni] # cumulative strategies
    u = 0.0
    nbr = 10
    Tx = round(Int64, T * 60 / nbr) # iterations for best response calculations
    tx = 1 # index for sigma calculations
    u1 = Float64[]
    σ1 = Vector{Vector{Vector{Float64}}}(undef, nbr)
    println("MCCFR: starting $T minutes...")
    tstart, tcurr, tbr = now(), now(), now()
    while tominute(tcurr - tstart) < T
        for player in players
            h0 = SVector(0, 0, 0, 0, 0, 0)
            u = updatetree!(h0, player, 1.0, 1.0, r, s, g, gs, 1)
            (player == 1) && push!(u1, u)
        end
        if tosecond(now() - tbr) > Tx # calculate sigma every Tx seconds
            σ1[tx] = strat_cfr(s)
            tx += 1
            tbr = now()
        end
        tcurr = now()
    end
    σ = strat_cfr(s)
    return u1, r, s, σ, σ1
end

function cfr_stoprule(T::Int64, g::AirDefenseGame, gs::GameSet; tolerance::Float64 = 5e-5)
    ni, ni_stage, na_stage, players = gs.ni, gs.ni_stage, gs.na_stage, gs.players
    r = [zeros(numactions(i, ni_stage, na_stage)) for i in 1:ni] # regret
    s = [zeros(numactions(i, ni_stage, na_stage)) for i in 1:ni] # cumulative strategies
    u = 0.0
    n_hist = 10 # number of intermediate strategies to record (for best responses)
    Tx = round(Int64, T * 60 / n_hist) # frequency of recording history
    tx = 1 # index for sigma calculations
    t = 1 # index for number of iter
    u1deriv = 0
    u1prev = 0
    converged = false
    u1 = Float64[]
    σ1 = Vector{Vector{Vector{Float64}}}(undef, n_hist)
    println("MCCFR: starting $T minutes at $(now())...")
    tstart, t_hist = now(), now(), now()
    while !converged && tominute(now() - tstart) < T
        for player in players
            h0 = SVector(0, 0, 0, 0, 0, 0)
            u = updatetree!(h0, player, 1.0, 1.0, r, s, g, gs, 1)
            (player == 1) && push!(u1, u)
        end
        if tosecond(now() - t_hist) > Tx # calculate sigma every Tx seconds
            σ1[tx] = strat_cfr(s)
            tx += 1
            t_hist = now()
        end
        u1deriv = u1deriv * (1 - 1 / t) + (u1[t] - u1prev) * (1 / t)
        converged = abs(u1deriv) < tolerance
        u1prev = u1[t]
        t += 1
        t % 10 == 0 && print("\rIter: $t")
    end
    println("")
    σ = strat_cfr(s)
    return u1, r, s, σ, σ1, converged
end

max_stars_repo_name: ajkeith/Cyber-Air-Defense

max_issues_repo_name: ajkeith/Cyber-Air-Defense


text:
#!/usr/bin/env julia

import Celeste.ParallelRun: one_node_joint_infer, infer_init, BoundingBox, find_neighbors
import Celeste.SDSSIO: RunCamcolField, load_field_images, PlainFITSStrategy

const datadir = joinpath(dirname(@__FILE__), "..", "..", "test", "data")
wd = pwd()
cd(datadir)
run(`make`)
cd(wd)

"""
This benchmark operates on a box of the sky that contains
seven light sources. During the optimization, some pixel is visited
254,771 times (``pixel visits'').
"""
function benchmark_seven_light_sources()
    # very small patch of sky that turns out to have 4 sources.
    # We checked that this patch is in the given field.
    box = BoundingBox(164.39, 164.41, 39.11, 39.13)
    rcfs = [RunCamcolField(3900, 6, 269),]

    catalog, target_sources, neighbor_map, images =
                        infer_init(rcfs, PlainFITSStrategy(datadir); box=box)
    ctni = (catalog, target_sources, neighbor_map, images)

    # Warm up---this compiles the code
    ctni2 = (catalog, target_sources[1:1], neighbor_map[1:1], images[1:1])
    one_node_joint_infer(ctni2...)

    # clear allocations in case julia is running with --track-allocations=user
    Profile.clear_malloc_data()

    if isempty(ARGS)
        @time one_node_joint_infer(ctni...)
    elseif ARGS[1] == "--profile"
        Profile.init(delay=0.01)
        @profile one_node_joint_infer(ctni...)
        Profile.print(format=:flat, sortedby=:count)
    end
end

benchmark_seven_light_sources()


text:
test_str = "nop +0
acc +1
jmp +4
acc +3
jmp -3
acc -99
acc +1
jmp -4
acc +6"

function acc_till_loop(list)
    acc = 0
    i = 1
    while list[i] != "∞"
        inst, num = split(list[i], ' ')
        list[i] = "∞"
        if inst == "acc"
            acc += parse(Int, num)
            i += 1
        elseif inst == "jmp"
            i += parse(Int, num)
        elseif inst == "nop"
            i += 1
        else
            throw(ErrorException("Invalid instruction `$inst $num` at line $i"))
        end
    end
    return acc, i
end

function acc_till_loop_alt(list)
    function parse_line(line)
        inst, num_str = split(line, ' ')
        return [inst, parse(Int, num_str)]
    end
    parsed_list = parse_line.(list)
    acc = 0
    i = 1
    while parsed_list[i][1] != "∞"
        inst, num = parsed_list[i]
        parsed_list[i][1] = "∞"
        if inst == "acc"
            acc += num
            i += 1
        elseif inst == "jmp"
            i += num
        elseif inst == "nop"
            i += 1
        else
            throw(ErrorException("Invalid instruction $line at line $i"))
        end
    end
    return acc, i
end

test_list = split(test_str, '\n')
@show acc_till_loop(test_list)[1] == 5

list = readlines("day8_input.txt")
@show acc_till_loop(list)[1] == 1489

test_list = split(test_str, '\n')
@show acc_till_loop_alt(test_list)[1] == 5

list = readlines("day8_input.txt")
@show acc_till_loop_alt(list)[1] == 1489

#=
# @btime doesn't work if function modifies input, so we use `copy(list)`, but
# in this case is okay, because `@btime copy(list)` yields about 633 ns, which
# does not harm the estimate below of tens of microseconds.
```julia
julia> @btime acc_till_loop!(copy(list))
  84.806 μs (1012 allocations: 75.45 KiB)
(1489, 164)

julia> @btime acc_till_loop2(list)
  673.271 μs (5918 allocations: 346.25 KiB)
(1489, 164)
```
 =#

# Part 2

function fix_code(list)
    function parse_line(line)
        inst, num_str = split(line, ' ')
        return [inst, parse(Int, num_str)]
    end
    len = length(list)
    attempt = 0
    parsed_list = parse_line.(list)
    visited = fill(false, len)
    i = 1
    j = 1
    acc = 0
    while i ≤ len
        inst, num = parsed_list[i]
        if visited[i]
            attempt += 1
            parsed_list = parse_line.(list)
            i = 1
            j = 1
            visited = fill(false, len)
            acc = 0
        else
            visited[i] = true
            if inst == "acc"
                acc += num
                i += 1
            elseif ((j == attempt && inst == "nop") || (j != attempt && inst == "jmp"))
                i += num
                j += 1
            else
                i += 1
                j += 1
            end
        end
    end
    return acc, i
end

test_list = split(test_str, '\n')
@show fix_code(test_list)[1] == 8

list = readlines("day8_input.txt")
@show fix_code(list)[1] == 1539

# As expected, part 2 is much slower than part 1.
# If we guess that they would make the "error" closer to the end, we could code backwards,
# but it would be much more work and memory intensive, to record where we have been and
# how much backwards we should go each time.
#=
```julia
julia> @btime fix_code(list)
  48.314 ms (383943 allocations: 22.19 MiB)
(1539, 602)
```
 =#

nothing


text:
# See interface docs
@inline inbounds(xs::Tuple, data::Union{GridData,SimData}) =
    inbounds(xs, gridsize(data), boundary(data))
@inline function inbounds(xs::Tuple, maxs::Tuple, boundary)
    a, inbounds_a = inbounds(xs[1], maxs[1], boundary)
    b, inbounds_b = inbounds(xs[2], maxs[2], boundary)
    (a, b), inbounds_a & inbounds_b
end
@inline function inbounds(x::Number, max::Number, boundary::Remove)
    x, isinbounds(x, max)
end
@inline function inbounds(x::Number, max::Number, boundary::Wrap)
    if x < oneunit(x)
        max + rem(x, max), true
    elseif x > max
        rem(x, max), true
    else
        x, true
    end
end

@inline isinbounds(x::Tuple, data::Union{SimData,GridData}) =
    isinbounds(x::Tuple, gridsize(data))
@inline isinbounds(xs::Tuple, maxs::Tuple) = all(isinbounds.(xs, maxs))

@inline isinbounds(x::Number, max::Number) = x >= one(x) && x <= max

#= Wrap boundary where required. This optimisation allows us to ignore
bounds checks on neighborhoods and still use a wraparound grid. =#
_handleboundary!(grids::Tuple) = map(_handleboundary!, grids)
_handleboundary!(griddata::GridData) = _handleboundary!(griddata, boundary(griddata))
_handleboundary!(griddata::GridData, ::Remove) = griddata
function _handleboundary!(griddata::GridData, ::Wrap)
    r = radius(griddata)
    r < 1 && return griddata

    # TODO optimise this. Its mostly a placeholder so wrapping still works in GOL tests.
    src = parent(source(griddata))
    nrows, ncols = gridsize(griddata)
    startpadrow = startpadcol = 1:r
    endpadrow = nrows+r+1:nrows+2r
    endpadcol = ncols+r+1:ncols+2r
    startrow = startcol = 1+r:2r
    endrow = nrows+1:nrows+r
    endcol = ncols+1:ncols+r
    rows = 1+r:nrows+r
    cols = 1+r:ncols+r

    # Left
    @inbounds copyto!(src, CartesianIndices((rows, startpadcol)),
                      src, CartesianIndices((rows, endcol)))
    # Right
    @inbounds copyto!(src, CartesianIndices((rows, endpadcol)),
                      src, CartesianIndices((rows, startcol)))
    # Top
    @inbounds copyto!(src, CartesianIndices((startpadrow, cols)),
                      src, CartesianIndices((endrow, cols)))
    # Bottom
    @inbounds copyto!(src, CartesianIndices((endpadrow, cols)),
                      src, CartesianIndices((startrow, cols)))

    # Copy four corners
    # Top Left
    @inbounds copyto!(src, CartesianIndices((startpadrow, startpadcol)),
                      src, CartesianIndices((endrow, endcol)))
    # Top Right
    @inbounds copyto!(src, CartesianIndices((startpadrow, endpadcol)),
                      src, CartesianIndices((endrow, startcol)))
    # Botom Left
    @inbounds copyto!(src, CartesianIndices((endpadrow, startpadcol)),
                      src, CartesianIndices((startrow, endcol)))
    # Botom Right
    @inbounds copyto!(src, CartesianIndices((endpadrow, endpadcol)),
                      src, CartesianIndices((startrow, startcol)))

    _wrapstatus!(sourcestatus(griddata))
end

_wrapstatus!(status::Nothing) = nothing
function _wrapstatus!(status::AbstractArray)
    # This could be further optimised.
    status[end-1, :] .|= status[1, :]
    status[:, end-1] .|= status[:, 1]
    #status[end-2, :] .|= status[1, :] .|= status[2, :]
    status[end-2, :] .= true
    status[:, end-2] .|= status[:, 1] .|= status[:, 2]
    #status[1, :] .|= status[end-2, :] .|= status[end-1, :]
    status[1, :] .= true
    status[:, 1] .|= status[:, end-2] .|= status[:, end-1]
    status .= true
end


text:
module Tensors

using Base.Cartesian
using TensorOperations

# Scalartype

export scalartype
scalartype(::Type{Any}) = throw(MethodError(scalartype, (Any,)))
scalartype(t::DataType) = scalartype(super(t))
scalartype(x) = scalartype(typeof(x))
scalartype{T <: Number}(::Type{T}) = T

# Typedefs

export Mode
immutable Mode
    mlabel::Any
    msize::Int
end

export Tag
immutable Tag{T,U} mlabel::U end

immutable Virtual{T}
    mlabel::T
    scale::Int
end

export Index
const Index = Dict{Any,Int}

export Tensor
type Tensor{T}
    modes::Vector{Mode}
    data::Vector{T}
end
Tensor{T}(mlabel::AbstractVector, x::AbstractArray{T}) = Tensor([Mode(k,n) for (k,n) in zip(mlabel,size(x))], vec(x))

# Output formatting

Base.show(io::IO, k::Mode   ) = (print(io,"Mode("); show(io,mlabel(k)); print(io,",",msize(k),")"))
Base.show{T}(io::IO, k::Tag{T}) = (print(io, T,"("); show(io,k.mlabel); print(io,")"))
Base.show(io::IO, k::Virtual) = (print(io,"Virtual("); show(io,k.mlabel); print(io,",",k.scale,")"))
Base.show(io::IO, t::Tensor ) = print(io,"Tensor{",eltype(t),"}([",join(map(string, mode(t)), ", "),"])")

# Mode

export msize,mlabel
msize(k::Mode) = k.msize
msize(K::AbstractVector{Mode}) = length(K) > 0 ? prod(map(msize,K)) : 1
mlabel(k::Mode) = k.mlabel
mlabel(K::AbstractVector{Mode}) = Any[k.mlabel for k in K]
import Base.==; ==(k::Mode, l::Mode) = mlabel(k) == mlabel(l) && msize(k) == msize(l)
multiplies(k::Any,l::Any) = k == l
replacemodes(k::Any, l::Any) = Pair[]

# Mode tags

export tag, untag
tag(T, k) = Tag{T,typeof(k)}(k)
tag(T, k::Mode) = Mode(tag(T,mlabel(k)),msize(k))
tag(T, K::AbstractVector) = Any[tag(T,k) for k in K]
tag(T, K::AbstractVector{Mode}) = Mode[tag(T,k) for k in K]
untag(k) = k
untag(k::Tag) = k.mlabel
untag(k::Mode) = Mode(untag(mlabel(k)), msize(k))
untag(K::AbstractVector) = Any[untag(k) for k in K]
untag(K::AbstractVector{Mode}) = Mode[untag(k) for k in K]
import Base.==; =={T}(k::Tag{T}, l::Tag{T}) = untag(k) == untag(l)

export istagged
istagged{T}(k::Tag{T}, TT) = T == TT
istagged(k::Any, TT) = false
istagged(k::Mode, TT) = istagged(k.mlabel,TT)

export tag!
function tag!(t::Tensor, T, modes)
    for k in modes
        i = findfirst(mlabel(t), k)
        if i == 0 error("Could not find mode $(k)!"); end
        t.modes[i] = tag(T, t.modes[i])
    end
    return t
end
tag(t::Tensor, T, modes) = tag!(copy(t), T, modes)

export untag!
function untag!(t::Tensor, modes::AbstractVector)
    for k in modes
        i = findfirst(untag(mlabel(t)), k)
        if i == 0 error("Could not find mode $k!"); end
        t.modes[i] = untag(t.modes[i])
    end
    return t
end
untag(t::Tensor, modes::AbstractVector) = untag!(copy(t), modes)

function untag!(t::Tensor, T)
    for (i,k) in enumerate(t.modes)
        if istagged(k,T)
            t.modes[i] = untag(k)
        end
    end
    return t
end
untag(t::Tensor, T) = untag!(copy(t), T)

export retag!, retag
function retag!(t::Tensor, p::Pair...)
    for (OldT,NewT) in p
        for (i,k) in enumerate(t.modes)
            if istagged(k,OldT)
                t.modes[i] = tag(NewT,untag(k))
            end
        end
    end
    return t
end
retag(t::Tensor, args...) = retag!(copy(t), args...)

export filtertags
function filtertags(K::AbstractVector{Any}, T)
    KK = Vector{Any}()
    for k in K
        if istagged(k,T)
            push!(KK,untag(k))
        end
    end
    return KK
end

# Row / column mode tags

export square
square(k) = [tag(:R,k); tag(:C,k)]
square(k::Mode) = Mode[tag(:R,k); tag(:C,k)]
square(K::AbstractVector{Any}) = Any[tag(:R,K); tag(:C,K)]
square(K::AbstractVector{Mode}) = Mode[tag(:R,K); tag(:C,K)]

multiplies(k::Tag{:C}, l::Tag{:R}) = multiplies(k.mlabel, l.mlabel)
multiplies(k::Any    , l::Tag{:R}) = multiplies(k       , l.mlabel)
multiplies(k::Tag{:C}, l::Any    ) = multiplies(k.mlabel, l       )

replacemodes(k::Tag{:C}, l::Tag{:R}) = Pair[]
replacemodes(k::Tag{:C}, l::Any    ) = Pair[tag(:R,untag(k)) => untag(k)]
replacemodes(k::Any    , l::Tag{:R}) = Pair[tag(:C,untag(k)) => untag(k)]

# Basic functions

export mode,msize
scalartype{T}(::Type{Tensor{T}}) = T
Base.eltype{T}(::Type{Tensor{T}}) = T
Base.length(t::Tensor) = length(t.data)
Base.start(t::Tensor) = start(t.data)
Base.next(t::Tensor,state) = next(t.data, state)
Base.done(t::Tensor,state) = done(t.data, state)
Base.ndims(t::Tensor) = length(t.modes)
Base.size(t::Tensor) = map(msize,t.modes)
Base.size(t::Tensor,k) = msize(t,k)
mode(t::Tensor) = t.modes
mode(t::Tensor, k::Any) = (i = findfirst(mlabel(t),k); i > 0 ? t.modes[i] : throw(ArgumentError("Tensor has no mode $k")))
mode(t::Tensor, K::AbstractVector) = Mode[mode(t,k) for k in K]
mlabel(t::Tensor) = mlabel(t.modes)
msize(t::Tensor) = size(t)
msize(t::Tensor, k) = msize(mode(t,k))
Base.copy(t::Tensor) = Tensor(copy(t.modes), copy(t.data))
Base.convert{T}(::Type{Tensor{T}}, t::Tensor) = Tensor(t.modes, convert(Vector{T}, t.data))

Base.isapprox(
    x::Tensor,y::Tensor,
    rtol::Real=Base.rtoldefault(scalartype(x),scalartype(y)),
    atol::Real=0
) = norm(x-y) <= atol + rtol*max(norm(x), norm(y))

export scalar
function scalar(t::Tensor)
    @assert ndims(t) == 0
    return t.data[1]
end

# Indexing

function linearindex(D::AbstractVector{Mode}, i::Index)
    n = 1; j = 0
    for k in D
        if !( 1 <= i[mlabel(k)] <= msize(k))
            throw(BoundsError(k,i[mlabel(k)]))
        end
        j = j + n*(i[mlabel(k)]-1)
        n *= msize(k)
    end
    return j+1
end
Base.getindex( t::Tensor, i::Index) = t.data[linearindex(t.modes,i)]
Base.setindex!(t::Tensor, ti, i::Index) = (t.data[linearindex(t.modes,i)] = ti; t)
Base.getindex( t::Tensor, i::Pair...) = t[Index(i...)]
Base.setindex!(t::Tensor, ti, i::Pair...) = t[Index(i...)] = ti

# Index iteration

immutable IndexSet D::Vector{Mode} end
index(D::AbstractVector{Mode}) = IndexSet(D)
index(t::Tensor) = index(mode(t))

Base.start(idx::IndexSet) = Index(map(m -> (mlabel(m),1), idx.D))
function Base.next(idx::IndexSet, i)
    ip = copy(i)
    done = true
    for k in idx.D
        if ip[mlabel(k)] < msize(k)
            ip[mlabel(k)] += 1
            done = false
            break
        end
        ip[mlabel(k)] = 1
    end
    if done return (i,Index()) end
    return (i,ip)
end
Base.done(idx::IndexSet, i) = isempty(i)
Base.eltype(::Type{IndexSet}) = Index
Base.length(idx::IndexSet) = msize(idx.D)

# Construction and Initialization

export empty
empty{T}(::Type{T}, D::AbstractVector{Mode}) = Tensor(D, Array{T}(msize(D)))
empty(D::AbstractVector{Mode}) = empty(Float64, D)
for f in (:zeros, :ones, :rand, :randn)
    @eval begin
        Base.$f{T}(::Type{T}, D::AbstractVector{Mode}) = Tensor{T}(D, $f(T, msize(D)))
        Base.$f(D::AbstractVector{Mode}) = $f(Float64, D)
    end
end
Base.eye{T}(::Type{T}, D::AbstractVector{Mode}) = Tensor(square(D), vec(eye(T, msize(D))))
Base.eye(D::AbstractVector{Mode}) = eye(Float64, D)
Base.one{T}(::Type{Tensor{T}}) = Tensor(Mode[],ones(T,1))
Base.one(t::Tensor) = one(scalartype(t))

export init
function init{T}(f::Function, ::Type{T}, D::AbstractVector{Mode})
    x = empty(T,D)
    for (il,i) in enumerate(index(D))
        x.data[il] = f(i)
    end
    return x
end
init(f::Function, D::AbstractVector{Mode}) = init(f, Float64, D)

# Storage order permutations

function permutedimsI!(t::Tensor, perm::AbstractVector{Int})
    if length(perm) != 0 && perm != collect(1:length(perm))
        data = similar(t.data)
        t.data = vec(TensorOperations.add!(
            one(scalartype(t)), reshape(t.data,size(t)...), Val{:N},
            zero(scalartype(t)), similar(t.data, size(t)[perm]...),
            perm
        ))
        t.modes = t.modes[perm]
    end
    return t
end
Base.permutedims!(t::Tensor, newmodes::AbstractVector) = permutedimsI!(t, findpermutation(mlabel(t.modes), newmodes))

# Unfoldings

function Base.getindex(t::Tensor, Ms...)
    permutedims!(t,vcat(Ms...))
    return reshape(t.data, [msize(t,M) for M in Ms]...)
end

# Reshaping

export splitm, splitm!
splitm(t::Tensor, s#=::Dict{Any,Vector{Mode}}=#) = splitm!(copy(t),s)
function splitm!(t::Tensor, s#=::Dict{Any,Vector{Mode}}=#)
    for k in mode(t)
        if msize(k) != msize(s[mlabel(k)])
            throw(ArgumentError("Invalid splitting for mode ", mlabel(k), "! (", msize(k), " != ", join("*",[msize(kk) for kk in s[mlabel(k)]]),")"))
        end
    end
    t.modes = vcat([s[mlabel(k)] for k in mode(t)]...)
    return t
end

export mergem, mergem!
mergem(t::Tensor, m#=::Dict{Vector, Any}=#) = mergem!(copy(t),m)
function mergem!(t::Tensor, m#=::Dict{Vector, Any}=#)
    K = vcat([K for K in keys(m)]...)
    DmK = setdiff(mlabel(t),K)
    permutedims!(t, vcat(K, DmK))
    n = (k = 1; [msize(t.modes[k:(k+=length(K))-1]) for K in keys(m)])
    t.modes = vcat(Mode[Mode(k,nk) for (k,nk) in zip(values(m),n)], mode(t,DmK))
    return t
end

export resize
function resize(t::Tensor, n)
    modes = [Mode(mlabel(k), get(n, mlabel(k), msize(k))) for k in t.modes]
    return Tensor(modes, vec(reshape(t.data, [nk for nk in size(t)]...)[[1:msize(k) for k in modes]...]))
end

# Concatenation and padding (for TN sum)

export padcat
@generated function padcat_(x::Array, y::Array, extendmode::AbstractVector{Bool})
    N = ndims(x)
    quote
        @assert ndims(x) == ndims(y) == length(extendmode)
        z = zeros(promote_type(eltype(x),eltype(y)), (collect(size(x)) + extendmode.*collect(size(y)))...)

        stridez_1 = 1; @nexprs $N d->(stridez_{d+1} = stridez_d*size(z, d))
        stridex_1 = 1; @nexprs $N d->(stridex_{d+1} = stridex_d*size(x, d))
        stridey_1 = 1; @nexprs $N d->(stridey_{d+1} = stridey_d*size(y, d))

        $(Symbol(:offsetz_, N)) = 1
        $(Symbol(:offsetx_, N)) = 1
        $(Symbol(:offsety_, N)) = 1

        @nexprs $N d->(nx_d = extendmode[d] ? size(x,d) : 0)

        @nloops $N i x d->(
            offsetz_{d-1} = offsetz_d + (i_d-1)*stridez_d;
            offsetx_{d-1} = offsetx_d + (i_d-1)*stridex_d
        ) begin
            @inbounds z[offsetz_0] = x[offsetx_0]
        end

        @nloops $N i y d->(
            offsetz_{d-1} = offsetz_d + (nx_d+i_d-1)*stridez_d;
            offsety_{d-1} = offsety_d + (i_d-1)*stridey_d
        ) begin
            @inbounds z[offsetz_0] = y[offsety_0]
        end

        return z
    end
end

function padcat(x::Tensor, y::Tensor, extendmodes)
    @assert length(extendmodes) > 0 "Must extend over at least one mode!"
    permutedims!(x, mlabel(y))
    zdata = padcat_(
        reshape(x.data,size(x)...), reshape(y.data,size(y)...),
        [k in extendmodes for k in mlabel(x)]
    )
    return Tensor(
        [Mode(k,n) for (k,n) in zip(mlabel(x), size(zdata))],
        vec(zdata)
    )
end

# Transposition and conjugation

Base.conj!(t::Tensor) = (conj!(t.data); return t)
Base.conj(t::Tensor) = conj!(copy(t))
function Base.transpose!(t::Tensor)
    for (i,k) in enumerate(t.modes)
        if istagged(k,:R) t.modes[i] = tag(:C,untag(k))
        elseif istagged(k,:C) t.modes[i] = tag(:R,untag(k))
        end
    end
    return t
end
Base.transpose(t::Tensor) = transpose!(copy(t))
Base.ctranspose!(t::Tensor) = transpose!(conj!(t))
Base.ctranspose(t::Tensor) = transpose!(conj(t))

# Vector arithmetic

for f in (:+,:-)
    @eval begin
        function Base.$f(t1::Tensor, t2::Tensor)
            permutedims!(t2, mlabel(t1))
            return Tensor(copy(t1.modes), $f(t1.data, t2.data))
        end
        Base.$f(a::Number, t::Tensor) = Tensor(copy(t.modes), $f(a, t.data))
        Base.$f(t::Tensor, a::Number) = Tensor(copy(t.modes), $f(t.data, a))
        Base.$f(t::Tensor) = Tensor(copy(t.modes), $f(t.data))
    end
end

import Base: *, /
*(a::Number, t::Tensor) = Tensor(copy(t.modes), a*t.data)
*(t::Tensor, a::Number) = a*t
/(t::Tensor, a::Number) = Tensor(copy(t.modes), t.data/a)

function Base.dot(t1::Tensor,t2::Tensor)
    permutedims!(t2, mlabel(t1))
    return dot(t1.data, t2.data)
end
Base.norm(t::Tensor) = vecnorm(t.data)

import Base: \
function \(A::Tensor, x::Tensor)
    KA = filtertags(mlabel(A), :R)
    Kx = [x.modes[findfirst(l -> multiplies(k,l), mlabel(x))] for k in KA]
    Kcx = setdiff(x.modes,Kx)
    return Tensor([Kx;Kcx], vec(A[tag(:R,KA), tag(:C,KA)]\x[mlabel(Kx),mlabel(Kcx)]))
end

function Base.inv(A::Tensor)
    D = filtertags(mlabel(A), :R)
    return Tensor(mode(A, square(D)), vec(inv(A[tag(:R,D),tag(:C,D)])))
end

# Unary elementwise arithmetic

for f in (:abs, :sin, :cos, :tan, :sqrt, :exp, :log)
    @eval (Base.$f)(t::Tensor, args...) = Tensor(copy(t.modes), $f(t.data, args...))
end
for f in (:mod, :mod1)
    @eval (Base.$f)(t::Tensor, m) = Tensor(copy(t.modes), $f(t.data, m))
end
for f in (:floor, :ceil, :base)
    @eval (Base.$f)(t::Tensor, opt::Integer...) = Tensor(copy(t.modes), $f(t.data, opt...))
end

# Binary elementwise arithmetic

# Unsupported for the time being, as overloading broadcast(!) correctly is complicated
# and it is not clear how much this functionality is used.

# for f in (:.*,:./)
#     @eval begin
#         Base.$f(t1::Tensor, t2::Tensor) = Tensor(copy(t1.modes), $f(t1.data, t2[mlabel(t1)]))
#         Base.$f(a::Number, t::Tensor) = Tensor(copy(t.modes), $f(a, t.data))
#         Base.$f(t::Tensor, a::Number) = Tensor(copy(t.modes), $f(t.data, a))
#     end
# end

# Mode product

function matchmodes(M1::AbstractVector{Mode},M2::AbstractVector{Mode})
    deleteindex! = (v::Vector, i) -> (v[i] = v[end]; pop!(v))
    M = Array{Mode}(0)
    K1 = Array{Mode}(0)
    K2 = Array{Mode}(0)
    N = copy(M2)
    for m1 in M1
        i2 = findfirst(m2->multiplies(mlabel(m1),mlabel(m2)), N)
        if i2 > 0
            push!(K1, m1)
            push!(K2, N[i2])
            deleteindex!(N, i2)
        else
            push!(M, m1)
        end
    end
    return M,K1,K2,N
end

function mergemodes(M, K1, K2, N)
    D = [M;N]
    for (k1,k2) in zip(K1,K2)
        for (old,new) in replacemodes(mlabel(k1),mlabel(k2))
            i = findfirst(mlabel(D), old)
            if i > 0
                D[i] = Mode(new, msize(D[i]))
            end
        end
    end
    return D
end

import Base.*
function *(t1::Tensor, t2::Tensor)
    M,K1,K2,N = matchmodes(t1.modes,t2.modes)
    D = mergemodes(M,K1,K2,N)
    return Tensor(D, vec(t1[mlabel(M),mlabel(K1)] * t2[mlabel(K2),mlabel(N)]))
end

# Scaling by diagonal matrices

@generated function _scale!(R::Array, A::Array, b::AbstractVector...)
    N = length(b)
    quote
        @assert size(R) == size(A)

        stride_1 = 1
        @nexprs $N d->(stride_{d+1} = stride_d*size(A, d))
        $(Symbol(:offset_, N)) = 1
        $(Symbol(:b_, N)) = 1

        @nloops $N i A d->(
            offset_{d-1} = offset_d + (i_d-1)*stride_d;
            b_{d-1} = b[d][i_d]*b_d
        ) begin
            @inbounds R[offset_0] = b_0*A[offset_0]
        end
        return A
    end
end
_scale_inplace!(A::Array, b::AbstractVector...) = _scale!(A,A,b...)

function splitAb(t::Tensor...)
    i = max(findfirst(t -> ndims(t) > 1, t), 1)
    A = t[i]
    b = [
        (j = findfirst(t -> multiplies(mlabel(t)[1], mlabel(k)), t[1:i-1])) > 0 ? t[j].data :
            (j = findfirst(t -> multiplies(mlabel(k), mlabel(t)[1]), t[i+1:end])) > 0 ? t[i+j].data :
                ones(msize(k))
        for k in mode(A)
    ]
    return A,b
end

function Base.scale!(t::Tensor...)
    A,b = splitAb(t...)
    _scale_inplace!(reshape(A.data, size(A)...), b...)
    return A
end
export scale
function scale(t::Tensor...)
    A,b = splitAb(t...)
    R = Tensor(copy(A.modes), Array{eltype(A)}(length(A)))
    _scale!(reshape(R.data, size(A)...), reshape(A.data, size(A)...), b...)
    return R
end

# QR decomposition

function Base.qr(t::Tensor, k::Any, r = maxrank(t,k))
    Dmk = setdiff(mlabel(t),[k])
    F = qrfact(t[Dmk,k])
    kk = Mode(k,r)
    return Tensor([mode(t,Dmk);kk], vec(F[:Q]*eye(eltype(t), msize(t,Dmk),r))),
           Tensor([tag(:R,kk),tag(:C,mode(t,k))], vec(F[:R][1:r,:]))
end
function Base.qr(t::Tensor, K::AbstractVector, k::Any, r = maxrank(t,K))
    DmK = setdiff(mlabel(t),K)
    F = qrfact(t[DmK,K])
    k = Mode(k,r)
    return Tensor([mode(t,DmK);k], vec(F[:Q]*eye(eltype(t), msize(t,DmK),r))),
           Tensor([k;mode(t,K)], vec(F[:R][1:r,:]))
end

# SVD decomposition

function Base.svd(t::Tensor, k::Any, rank = maxrank())
    Kc = setdiff(mlabel(t), [k])
    F = svdfact(t[Kc,k])
    kk = Mode(k, rank(F[:S]))
    return Tensor([mode(t,Kc); kk], vec(F[:U][:,1:msize(kk)])),
           Tensor([kk], F[:S][1:msize(kk)]),
           Tensor([tag(:R,kk);tag(:C,mode(t,k))], vec(F[:Vt][1:msize(kk), :]))
end
function Base.svd(t::Tensor, K::AbstractVector, k::Any, rank = maxrank())
    Kc = setdiff(mlabel(t), K)
    F = svdfact(t[Kc,K])
    k = Mode(k, rank(F[:S]))
    return Tensor([mode(t,Kc); k], vec(F[:U][:,1:msize(k)])),
           Tensor([k], F[:S][1:msize(k)]),
           Tensor([k;mode(t,K)], vec(F[:Vt][1:msize(k), :]))
end

# Common rank functions

export adaptive, fixed, maxrank
adaptive(eps; rel = true) = (s) -> max(length(s) - searchsortedlast(cumsum(reverse(s).^2), (eps*(rel ? vecnorm(s) : 1))^2),1)
fixed(r) = (s) -> r
maxrank() = (s) -> length(s)

# Helpers

# Determine perm such that after == before[perm]
function findpermutation(before, after)
    if length(before) != length(after)
        throw(ArgumentError("The number of modes must be preserved during permutation!"))
    end

    perm = Array{Int}(length(after))
    for i in 1:length(after)
        perm[i] = findfirst(before, after[i])
        if perm[i] == 0
            throw(ArgumentError("Could not find mode $(k)!"))
        end
    end
    return perm
end

maxrank(t::Tensor, K) = min(div(length(t),msize(t,K)), msize(t,K))

end # module


text:
#
# Define output path string:
#
# called by a driver right before run()
#
# ex.
#    #Create unique output path directory:
#    OUTPATH = IOstrings_outpath_name(problem_name, grid_resolution)
#
#
function IOstrings_outpath_name(problem_name, grid_resolution)
    #
    # Arguments:
    # problem_name = "NAME"
    # grid_resolution[1]   = Δx
    # grid_resolution[2]   = Δy
    # grid_resolution[end] = Δz
    #
    ndim = length(grid_resolution)

    outpath_string = string(grid_resolution[1], "mx")
    for i = 2:ndim
        ds = grid_resolution[i]
        outpath_string = string(outpath_string, ds, "mx")
    end
    current_time = string(Dates.format(convert(Dates.DateTime, Dates.now()), Dates.dateformat"yyyymmdd_HHMMSS"))
    OUTPATH = string("./output/",problem_name, "/", outpath_string,"_", current_time)

    mkpath(OUTPATH)

    return OUTPATH
end

max_stars_repo_name: akshaysridhar/CLIMA

max_issues_repo_name: akshaysridhar/CLIMA

max_forks_repo_name: akshaysridhar/CLIMA


text:
using ConicNonlinearBridge
using MathProgBase
using Ipopt

using Compat.Test
using Compat.LinearAlgebra
using Compat.SparseArrays

if VERSION < v"0.7.0-"
    mpb_path = Pkg.dir("MathProgBase")
end

if VERSION > v"0.7.0-"
    mpb_path = joinpath(dirname(pathof(MathProgBase)), "..")
end

include(joinpath(mpb_path, "test", "conicinterface.jl"))

@testset "ConicNonlinearBridge Tests" begin

@testset "linear and exponential cone tests (remove_single_rows is $remove_single_rows)" for remove_single_rows in [true, false]
    solver = ConicNLPWrapper(nlp_solver=IpoptSolver(print_level=0), remove_single_rows=remove_single_rows)
    coniclineartest(solver; duals=false, tol=1e-6)
    conicEXPtest(solver; duals=false, tol=1e-6)
end

@testset "second-order cone tests (soc_as_quadratic is $soc_as_quadratic)" for soc_as_quadratic in [true, false]
    solver = ConicNLPWrapper(nlp_solver=IpoptSolver(print_level=0), soc_as_quadratic=soc_as_quadratic, disaggregate_soc=false)
    conicSOCRotatedtest(solver; duals=false, tol=1e-6)

    @testset "disaggregation tests (disaggregate_soc is $disaggregate_soc)" for disaggregate_soc in [true, false]
        solver = ConicNLPWrapper(nlp_solver=IpoptSolver(print_level=0), soc_as_quadratic=soc_as_quadratic, disaggregate_soc=disaggregate_soc)
        conicSOCtest(solver; duals=false, tol=1e-6)
    end
end

@testset "miscellaneous tests" begin
    # min 0x - 2y - 1z
    #  st  x == 1,
    #      x >= norm(y, z) as constraint cone
    #      x continuous, y binary, z integer
    solver = ConicNLPWrapper(nlp_solver=IpoptSolver(print_level=0))
    @test issubset([:Free, :Zero, :NonNeg, :NonPos, :SOC, :SOCRotated, :ExpPrimal], MathProgBase.supportedcones(solver))
    m = MathProgBase.ConicModel(solver)
    MathProgBase.loadproblem!(m, [0, -2, -1], [1 0 0; Matrix(-1.0I, 3, 3)], [1, 0, 0, 0], [(:Zero, 1), (:SOC, 2:4)], [(:Free, 1:3)])
    MathProgBase.setvartype!(m, [:Cont, :Bin, :Int])
    MathProgBase.setwarmstart!(m, [3, 4, 2])
    @test m.nlp_model.colCat[1:3] == [:Cont, :Bin, :Int]
    @test m.solution[1:3] == [3, 4, 2]
    @test MathProgBase.numvar(m) == 3
    @test MathProgBase.numconstr(m) == 4

    # make variables continuous and solve
    MathProgBase.setvartype!(m, [:Cont, :Cont, :Cont])
    MathProgBase.optimize!(m)
    @test MathProgBase.status(m) == :Optimal
    @test MathProgBase.getsolution(m) ≈ [1.0, 0.894427, 0.447214] atol=1e-4 rtol=1e-4
    @test MathProgBase.getobjval(m) ≈ -2.236067 atol=1e-4 rtol=1e-4
    MathProgBase.freemodel!(m)
end

end


text:
"""
    ErrorDiffusion(filter::AbstractMatrix; clamp_error=true)

Generalized error diffusion algorithm. When calling `dither` using a color palette `cs`,
this will iterate over pixels and round them to the closest color in `cs`.
The rounding error is then "diffused" over the neighborhood defined by the matrix `filter`.
This diffused error can additionally be clamped to ``[0, 1]`` by setting
`clamp_error = true` (default).

When calling `dither` on a grayscale image without specifying a palette, `ErrorDiffusion`
algorithms will default to settings for binary dithering: `clamp_error=false` and the metric
`BinaryDitherMetric()`, which simply rounds to the closest binary number.

# Example
```julia-repl
julia> alg = FloydSteinberg() # returns ErrorDiffusion instance
DitherPunk.ErrorDiffusion{OffsetArrays.OffsetMatrix{Rational{Int64}, Matrix{Rational{Int64}}}}(Rational{Int64}[0//1 0//1 7//16; 3//16 5//16 1//16], true)

julia> cs = ColorSchemes.PuOr_7.colors; # using ColorSchemes.jl for color palette presets

julia> dither!(img, alg, cs);
```
"""
struct ErrorDiffusion{T<:AbstractMatrix} <: AbstractDither
    filter::T
    clamp_error::Bool
end
ErrorDiffusion(filter; clamp_error=true) = ErrorDiffusion(filter, clamp_error)

function binarydither!(alg::ErrorDiffusion, out::GenericGrayImage, img::GenericGrayImage)
    # this function does not yet support OffsetArray
    require_one_based_indexing(img)

    # Change from normalized intensities to Float as error will get added!
    # eagerly promote to the same type to make loop run faster
    FT = floattype(eltype(img))
    img = FT.(img)
    filter = eltype(FT).(alg.filter)

    drs = axes(alg.filter, 1)
    dcs = axes(alg.filter, 2)

    FT0, FT1 = FT(0), FT(1)

    @inbounds for r in axes(img, 1)
        for c in axes(img, 2)
            px = img[r, c]
            alg.clamp_error && (px = clamp01(px))

            px >= 0.5 ? (col = FT1) : (col = FT0) # round to closest color
            out[r, c] = col # apply pixel to dither
            err = px - col  # diffuse "error" to neighborhood in filter

            for dr in drs
                for dc in dcs
                    if (r + dr) in axes(img, 1) && (c + dc) in axes(img, 2)
                        img[r + dr, c + dc] += err * filter[dr, dc]
                    end
                end
            end
        end
    end

    return out
end

function colordither(
    alg::ErrorDiffusion,
    img::GenericImage,
    cs::AbstractVector{<:Pixel},
    metric::DifferenceMetric,
)
    # this function does not yet support OffsetArray
    require_one_based_indexing(img)

    index = Matrix{UInt8}(undef, size(img)...) # allocate matrix of color indices

    # Change from normalized intensities to Float as error will get added!
    # Eagerly promote to the same type to make loop run faster.
    FT = floattype(eltype(eltype(img))) # type of Float
    CT = floattype(eltype(img)) # type of colorant

    img = CT.(img)
    cs = CT.(cs)
    labcs = Lab{FT}.(cs) # otherwise each call to colordiff converts cs to Lab
    filter = FT.(alg.filter)

    drs = axes(alg.filter, 1)
    dcs = axes(alg.filter, 2)

    @inbounds for r in axes(img, 1)
        for c in axes(img, 2)
            px = img[r, c]
            alg.clamp_error && (px = clamp01(px))

            colorindex = _closest_color_idx(px, labcs, metric)
            index[r, c] = colorindex
            err = px - cs[colorindex]  # diffuse "error" to neighborhood in filter

            for dr in drs
                for dc in dcs
                    if (r + dr) in axes(img, 1) && (c + dc) in axes(img, 2)
                        img[r + dr, c + dc] += err * filter[dr, dc]
                    end
                end
            end
        end
    end
    return index
end

"""
    SimpleErrorDiffusion()

Error diffusion algorithm using the filter
```
*   1
1   0         (1//2)
```

[1]  Floyd, R.W. and L. Steinberg, "An Adaptive Algorithm for Spatial Gray
     Scale."  SID 1975, International Symposium Digest of Technical Papers,
     vol 1975m, pp. 36-37.
"""
SimpleErrorDiffusion() = ErrorDiffusion(OffsetMatrix([0 1; 1 0]//2, 0:1, 0:1))

"""
    FloydSteinberg()

Error diffusion algorithm using the filter
```
    *   7
3   5   1     (1//16)
```

[1]  Floyd, R.W. and L. Steinberg, "An Adaptive Algorithm for Spatial Gray
     Scale."  SID 1975, International Symposium Digest of Technical Papers,
     vol 1975m, pp. 36-37.
"""
function FloydSteinberg(; kwargs...)
    return ErrorDiffusion(OffsetMatrix([0 0 7; 3 5 1]//16, 0:1, -1:1); kwargs...)
end

"""
    JarvisJudice()

Error diffusion algorithm using the filter
```
        *   7   5
3   5   7   5   3
1   3   5   3   1   (1//48)
```
Also known as the Jarvis, Judice, and Ninke filter.

[1]  Jarvis, J.F., C.N. Judice, and W.H. Ninke, "A Survey of Techniques for
     the Display of Continuous Tone Pictures on Bi-Level Displays," Computer
     Graphics and Image Processing, vol. 5, pp. 13-40, 1976.
"""
function JarvisJudice(; kwargs...)
    return ErrorDiffusion(
        OffsetMatrix([0 0 0 7 5; 3 5 7 5 3; 1 3 5 3 1]//48, 0:2, -2:2); kwargs...
    )
end

"""
    Stucki()

Error diffusion algorithm using the filter
```
        *   8   4
2   4   8   4   2
1   2   4   2   1   (1//42)
```

[1]  Stucki, P., "MECCA - a multiple-error correcting computation algorithm
     for bilevel image hardcopy reproduction."  Research Report RZ1060, IBM
     Research Laboratory, Zurich, Switzerland, 1981.
"""
function Stucki(; kwargs...)
    return ErrorDiffusion(
        OffsetMatrix([0 0 0 8 4; 2 4 8 4 2; 1 2 4 2 1]//42, 0:2, -2:2); kwargs...
    )
end

"""
    Burkes()

Error diffusion algorithm using the filter
```
        *   8   4
2   4   8   4   2
1   2   4   2   1   (1//42)
```

[1] Burkes, D., "Presentation of the Burkes error filter for use in preparing
    continuous-tone images for presentation on bi-level devices." Unpublished, 1988.
"""
function Burkes(; kwargs...)
    return ErrorDiffusion(OffsetMatrix([0 0 0 8 4; 2 4 8 4 2]//32, 0:1, -2:2); kwargs...)
end

"""
    Sierra()

Error diffusion algorithm using the filter
```
        *   5   3
2   4   5   4   2
    2   3   2       (1//32)
```
Also known as Sierra3 or three-row Sierra due to the filter shape.
"""
function Sierra(; kwargs...)
    return ErrorDiffusion(
        OffsetMatrix([0 0 0 5 3; 2 4 5 4 2; 0 2 3 2 0]//32, 0:2, -2:2); kwargs...
    )
end

"""
    TwoRowSierra()

Error diffusion algorithm using the filter
```
        *   4   3
1   2   3   2   1   (1//16)
```
Also known as Sierra2.
"""
function TwoRowSierra(; kwargs...)
    return ErrorDiffusion(OffsetMatrix([0 0 0 4 3; 1 2 3 2 1]//16, 0:1, -2:2); kwargs...)
end

"""
    SierraLite()

Error diffusion algorithm using the filter
```
    *   2
1   1               (1//4)
```
Also known as Sierra-2-4A filter.
"""
function SierraLite(; kwargs...)
    return ErrorDiffusion(OffsetMatrix([0 0 2; 1 1 0]//4, 0:1, -1:1); kwargs...)
end

"""
    Atkinson()

Error diffusion algorithm using the filter
```
    *   1   1
1   1   1
    1               (1//8)
```
"""
function Atkinson(; kwargs...)
    return ErrorDiffusion(
        OffsetMatrix([0 0 1 1; 1 1 1 0; 0 1 0 0]//8, 0:2, -1:2); kwargs...
    )
end

"""
    Fan93()

Error diffusion algorithm using the filter
```
      *  7
1  3  5            (1//16)
```
A modification of the weights used in the Floyd-Steinberg algorithm.

[1] Z. Fan, "A Simple Modification of Error Diffusion Weights",
    IS&T's 46th Annual Conference, May 9-14, 1993, Final Program and Advanced Printing of
    Paper Summaries, pp 113-115 (1993).
"""
function Fan93(; kwargs...)
    return ErrorDiffusion(OffsetMatrix([0 0 0 7; 1 3 5 0]//16, 0:1, -2:1); kwargs...)
end

"""
    ShiauFan()

Error diffusion algorithm using the filter
```
        *   4
1   1   2           (1//8)
```

[1]  J. N. Shiau and Z. Fan. "Method for quantization gray level pixel data with extended
     distribution set", US 5353127A, United States Patent and Trademark Office, Oct. 4, 1993
"""
function ShiauFan(; kwargs...)
    return ErrorDiffusion(OffsetMatrix([0 0 0 4; 1 1 2 0]//8, 0:1, -2:1); kwargs...)
end

"""
    ShiauFan2()

Error diffusion algorithm using the filter
```
            *   8
1   1   2   4       (1//16)
```

[1]  J. N. Shiau and Z. Fan. "Method for quantization gray level pixel data with extended
     distribution set", US 5353127A, United States Patent and Trademark Office, Oct. 4, 1993
[2]  J. N. Shiau and Z. Fan. "A set of easily implementable coefficients in error-diffusion
     with reduced worm artifacts" Color Imaging: Device-Independent Color, Color Hard Copy,
     and Graphics Arts, volume 2658, pages 222–225. SPIE, March 1996.
"""
function ShiauFan2(; kwargs...)
    return ErrorDiffusion(OffsetMatrix([0 0 0 0 8; 1 1 2 4 0]//16, 0:1, -3:1); kwargs...)
end

"""
    FalseFloydSteinberg()

Error diffusion algorithm using the filter
```
*   3
3   2         (1//8)
```
Occasionally, you will see this filter erroneously called the Floyd-Steinberg filter.

# Note
There is no reason to use this algorithm, which is why DitherPunk doesn't export it.

[1] http://www.efg2.com/Lab/Library/ImageProcessing/DHALF.TXT
"""
function FalseFloydSteinberg(; kwargs...)
    return ErrorDiffusion(OffsetMatrix([0 3; 3 2]//8, 0:1, 0:1); kwargs...)
end


text:
using SMTPClient, JSON, Dates

# 設定ファイルの作成
function initsmtpsetting()
    setting = Dict(
        "url" => "smtps://smtp.gmail.com:465",
        "from" => "You",
        "username" => "you@gmail.com",
        "passwd" => "yourgmailpassword"
    )
    json = JSON.json(setting)
    open("./setting/smtp.json", "w") do f
        print(f, json)
    end
end

# 設定ファイルの読み込み
function readsmtpsetting()
    str = open(f -> read(f, String), "./setting/smtp.json")
    setting = JSON.parse(str)
    return setting
end

# タイトルと本文を指定すると，自分のメールアカウントを使って自分にメールを送る
function sendmail(subject::String, message::String)
    setting = readsmtpsetting()
    opt = SendOptions(
        isSSL=true,
        username=setting["username"],
        passwd=setting["passwd"])
    # Provide the message body as RFC5322 within an IO
    dtstr = Dates.format(Dates.now(), "e, dd u YYYY HH:MM:SS +0900")
    body = IOBuffer(
        "Date: " * dtstr * "\r\n" *
        "From: " * setting["from"] * " <" * setting["username"] * ">\r\n" *
        "To: " * setting["username"] * "\r\n" *
        "Subject: " * subject * "\r\n" *
        "\r\n" *
        message *
        "\r\n")
    url = setting["url"]
    rcpt = ["<" * setting["username"] * ">"]
    from = "<" * setting["username"] * ">"
    resp = send(url, rcpt, from, body, opt)
    return resp
end