text:
#######################################################################
#
# An example of creating Excel Doughnut charts with Python and XlsxWriter.
#
# The demo also shows how to set segment colors. It is possible to
# define chart colors for most types of XlsxWriter charts
# via the add_series() method. However, Pie/Doughnut charts are a special
# case since each segment is represented as a point so it is necessary to
# assign formatting to each point in the series.

# Original Python Copyright 2013-2016, John McNamara, jmcnamara@cpan.org
# https://github.com/jmcnamara/XlsxWriter

using XlsxWriter

function test()

	wb = Workbook("chart_doughnut.xlsx")
	ws = add_worksheet!(wb)
	bold = add_format!(wb, Dict("bold"=> 1))

	# Add the worksheet data that the charts will refer to.
	headings = ["Category" "Values"]
	data = [
		["Glazed" "Chocolate" "Cream"];
		[50 35 15];
	]

	write_row!(ws, "A1", headings, fmt=bold)
	write_column!(ws, "A2", data[1,:])
	write_column!(ws, "B2", data[2,:])

	#######################################################################
	#
	# Create a new chart object.
	#
	chart1 = add_chart!(wb, Dict("type"=> "doughnut"))

	# Configure the series. Note the use of the list syntax to define ranges:
	add_series!(chart1, Dict(
		"name"=>       "Doughnut sales data",
		"categories"=> ["Sheet1", 1, 0, 3, 0],
		"values"=>     ["Sheet1", 1, 1, 3, 1]
	))

	# Add a title.
	set_title!(chart1, Dict("name"=> "Popular Doughnut Types"))

	# Set an Excel chart style. Colors with white outline and shadow.
	set_style!(chart1, 10)

	# Insert the chart into the worksheet (with an offset).
	insert_chart!(ws, "C2", chart1, Dict("x_offset"=> 25, "y_offset"=> 10))

	#######################################################################
	#
	# Create a Doughnut chart with user defined segment colors.
	#

	# Create an example Doughnut chart like above.
	chart2 = add_chart!(wb, Dict("type"=> "doughnut"))

	# Configure the series and add user defined segment colors.
	add_series!(chart2, Dict(
		"name"=> "Doughnut sales data",
		"categories"=> "=Sheet1!\$A\$2:\$A\$4",
		"values"=>     "=Sheet1!\$B\$2:\$B\$4",
		"points"=> [
			Dict("fill"=> Dict("color"=> "#FA58D0")),
			Dict("fill"=> Dict("color"=> "#61210B")),
			Dict("fill"=> Dict("color"=> "#F5F6CE")),
		];
	))

	# Add a title.
	set_title!(chart2, Dict("name"=> "Doughnut Chart with user defined colors"))

	# Insert the chart into the worksheet (with an offset).
	insert_chart!(ws, "C18", chart2, Dict("x_offset"=> 25, "y_offset"=> 10))

	#######################################################################
	#
	# Create a Doughnut chart with rotation of the segments.
	#

	# Create an example Doughnut chart like above.
	chart3 = add_chart!(wb, Dict("type"=> "doughnut"))

	# Configure the series.
	add_series!(chart3, Dict(
		"name"=> "Doughnut sales data",
		"categories"=> "=Sheet1!\$A\$2:\$A\$4",
		"values"=>     "=Sheet1!\$B\$2:\$B\$4",
	))

	# Add a title.
	set_title!(chart3, Dict("name"=> "Doughnut Chart with segment rotation"))

	# Change the angle/rotation of the first segment.
	set_rotation!(chart3, 90)

	# Insert the chart into the worksheet (with an offset).
	insert_chart!(ws, "C34", chart3, Dict("x_offset"=> 25, "y_offset"=> 10))

	#######################################################################
	#
	# Create a Doughnut chart with user defined hole size and other options.
	#

	# Create an example Doughnut chart like above.
	chart4 = add_chart!(wb, Dict("type"=> "doughnut"))

	# Configure the series.
	add_series!(chart4, Dict(
		"name"=> "Doughnut sales data",
		"categories"=> "=Sheet1!\$A\$2:\$A\$4",
		"values"=>     "=Sheet1!\$B\$2:\$B\$4",
		"points"=> [
			Dict("fill"=> Dict("color"=> "#FA58D0")),
			Dict("fill"=> Dict("color"=> "#61210B")),
			Dict("fill"=> Dict("color"=> "#F5F6CE")),
		];
	))

	# Set a 3D style.
	set_style!(chart4, 26)

	# Add a title.
	set_title!(chart4, Dict("name"=> "Doughnut Chart with options applied"))

	# Change the angle/rotation of the first segment.
	set_rotation!(chart4, 28)

	# Change the hole size.
	set_hole_size!(chart4, 33)

	# Insert the chart into the worksheet (with an offset).
	insert_chart!(ws, "C50", chart4, Dict("x_offset"=> 25, "y_offset"=> 10))

	close(wb)
	isfile("chart_doughnut.xlsx")
end

test()


text:
import AlphaZero.GI

using Printf
using Crayons
using StaticArrays

const NUM_HOUSES_PER_PLAYER = 6

const NUM_INITIAL_SEEDS_PER_HOUSE = 3

const NUM_HOUSES = 2 * NUM_HOUSES_PER_PLAYER

const NUM_SEEDS =  NUM_HOUSES * NUM_INITIAL_SEEDS_PER_HOUSE

const Player = Int
const WHITE = 1
const BLACK = 2

other(p::Player) = 3 - p

struct Board
  stores :: SVector{2, UInt8}
  houses :: SMatrix{2, NUM_HOUSES_PER_PLAYER, UInt8, NUM_HOUSES}
end

const INITIAL_BOARD =
  let stores = SVector(0, 0)
  let houses = @SMatrix[
      NUM_INITIAL_SEEDS_PER_HOUSE
      for i in 1:2, j in 1:NUM_HOUSES_PER_PLAYER]
    Board(stores, houses)
  end end

mutable struct Game
  board :: Board
  curplayer :: Player
  function Game(board::Board=INITIAL_BOARD, white_playing::Bool=true)
    new(board, white_playing ? WHITE : BLACK)
  end
end

GI.Board(::Type{Game}) = Board
GI.Action(::Type{Game}) = Int

Base.copy(g::Game) = Game(g.board, g.curplayer == WHITE)

#####
##### Position system to ease board manipulation
#####

struct HousePos
  player :: Int
  num :: Int
end

struct StorePos
  player :: Int
end

const Pos = Union{HousePos, StorePos}

function next_pos(pos, player)
  if isa(pos, HousePos)
    if pos.num > 1
      return HousePos(pos.player, pos.num - 1)
    elseif pos.num == 1
      if pos.player == player
        return StorePos(player)
      else
        return HousePos(player, NUM_HOUSES_PER_PLAYER)
      end
    else
      @assert false
    end
  elseif isa(pos, StorePos)
    @assert pos.player == player
    return HousePos(other(player), NUM_HOUSES_PER_PLAYER)
  else
    @assert false
  end
end

function read_pos(b::Board, pos)
  isa(pos, HousePos) ?
    b.houses[pos.player, pos.num] :
    b.stores[pos.player]
end

function write_pos(b::Board, pos, v)
  if isa(pos, HousePos)
    houses = setindex(b.houses, v, pos.player, pos.num)
    return Board(b.stores, houses)
  else
    stores = setindex(b.stores, v, pos.player)
    return Board(stores, b.houses)
  end
end

#####
##### Defining game rules
#####

GI.available_actions(g::Game, player) = findall(>(0), g.board.houses[player,:])

GI.available_actions(g::Game) = GI.available_actions(g, g.curplayer)

sum_houses(b::Board, player) = sum(b.houses[player,:])

function GI.play!(g::Game, a)
  pos = HousePos(g.curplayer, a)
  nseeds = read_pos(g.board, pos)
  @assert nseeds > 0
  g.board = write_pos(g.board, pos, 0)
  for i in 1:nseeds
    pos = next_pos(pos, g.curplayer)
    g.board = write_pos(g.board, pos, read_pos(g.board, pos) + 1)
  end
  # Check endgame
  if sum_houses(g.board, g.curplayer) == 0
    other_sum = sum_houses(g.board, other(g.curplayer))
    new_store = g.board.stores[g.curplayer] + other_sum
    stores = setindex(g.board.stores, new_store, g.curplayer)
    g.board = Board(stores, zero(typeof(g.board.houses)))
  # Free turn is last seed was put in a store
  elseif isa(pos, HousePos)
    g.curplayer = other(g.curplayer)
  end
  return
end

GI.board(g::Game) = g.board

function GI.board_symmetric(g::Game) :: Board
  b = g.board
  stores = @SVector[b.stores[2], b.stores[1]]
  houses = vcat(b.houses[2,:]', b.houses[1,:]')
  Board(stores, houses)
end

GI.white_playing(g::Game) = g.curplayer == WHITE

game_terminated(g::Game) = all(==(0), g.board.houses)

#####
##### Reward shaping
#####

function zero_one_reward(nw, nb)
  nw > nb && return 1.
  nw < nb && return -1.
  return 0.
end

linear_reward(nw, nb) = (nw - nb) / NUM_SEEDS

function GI.white_reward(g::Game)
  if game_terminated(g)
    nw, nb = g.board.stores
    return zero_one_reward(nw, nb) # + linear_reward(nw, nb)
  else
    return nothing
  end
end

#####
##### ML interface
#####

GI.num_actions(::Type{Game}) = NUM_HOUSES_PER_PLAYER

GI.action_id(::Type{Game}, a) = a

GI.action(::Type{Game}, id) = id

function GI.vectorize_board(::Type{Game}, board)
  function cell(pos, chan)
    if chan == :nstones
      return read_pos(board, pos)
    elseif chan == :whouse
      return isa(pos, HousePos) && pos.player == WHITE
    elseif chan == :wstore
      return isa(pos, StorePos) && pos.player == WHITE
    elseif chan == :bhouse
      return isa(pos, HousePos) && pos.player == BLACK
    elseif chan == :bstore
      return isa(pos, StorePos) && pos.player == BLACK
    end
    @assert false
  end
  positions = [
    [HousePos(WHITE, i) for i in NUM_HOUSES_PER_PLAYER:-1:1];
    [StorePos(WHITE)];
    [HousePos(BLACK, i) for i in NUM_HOUSES_PER_PLAYER:-1:1];
    [StorePos(BLACK)]]
  return Float32[
    cell(p, c)
    for p in positions,
        y in 1:1,
        c in [:nstones, :whouse, :wstore, :bhouse, :bstore]]
end

#####
##### User interface
#####

GI.action_string(::Type{Game}, a) = string(a)

function GI.parse_action(g::Game, str)
  try
    p = parse(Int, str)
    1 <= p <= NUM_HOUSES_PER_PLAYER ? p : nothing
  catch
    nothing
  end
end

#          1    2    3    4    5    6
#  +----+----+----+----+----+----+----+----+
#  |    |  3 |  3 |  3 |  3 |  3 |  3 |    |
#  |    +----+----+----+----+----+----+    |
#  |    |  3 |  3 |  3 |  3 |  3 |  3 |    |
#  +----+----+----+----+----+----+----+----+
#          6    5    4    3    2    1

function GI.print_state(g::Game, with_position_names=true)
  b = g.board
  gray(s...) = print(crayon"blue", s..., crayon"reset")
  white(s...) = print(crayon"white", s..., crayon"reset")
  bold(s...) = print(crayon"yellow", s..., crayon"reset")
  blank_cell() = repeat(" ", 4)
  num_cell(n) = n == 0 ? blank_cell() : @sprintf(" %2d ", n)
  function show_labels(labels)
    print(" "); print(blank_cell()); print(" ")
    for i in labels
      gray(num_cell(i)); print(" ")
    end
    print("\n")
  end
  function hline(len; newline=true)
    for i in 1:len
      gray("+----")
    end
    newline && gray("+\n")
  end
  function print_houses(houses, self)
    printing = self ? bold : white
    for v in houses
      gray("|"); printing(num_cell(v))
    end
  end
  function print_half(houses, self=false)
    gray("|"); print(blank_cell())
    print_houses(houses, self)
    gray("|"); print(blank_cell())
    gray("|\n")
  end
  function print_middle(stores)
    gray("|"); white(num_cell(stores[2]))
    hline(NUM_HOUSES_PER_PLAYER, newline=false)
    gray("+"); white(num_cell(stores[1]))
    gray("|\n")
  end
  show_labels(1:NUM_HOUSES_PER_PLAYER)
  hline(NUM_HOUSES_PER_PLAYER + 2)
  print_half(b.houses[2,:], g.curplayer == BLACK)
  print_middle(b.stores)
  print_half(reverse(b.houses[1,:]), g.curplayer == WHITE)
  hline(NUM_HOUSES_PER_PLAYER + 2)
  show_labels(reverse(1:NUM_HOUSES_PER_PLAYER))
end

function GI.read_state(::Type{Game})
  try
    function read_houses(player)
      print("Player $(player) houses: ")
      houses = [parse(Int, s) for s in split(readline())]
      @assert all(>=(0), houses)
      @assert length(houses) == NUM_HOUSES_PER_PLAYER
      return houses
    end
    function read_store(player)
      print("Player $(player) store: ")
      store = parse(Int, readline())
      @assert store >= 0
      return store
    end
    # Read board
    h1 = read_houses(1)
    h2 = read_houses(2)
    s1 = read_store(1)
    #s2 = read_store(2)
    nsofar = sum(h1) + sum(h2) + s1
    @assert nsofar <= NUM_SEEDS
    s2 = NUM_SEEDS - nsofar
    board = Board([s1, s2], [h1'; h2'])
    # Read current player
    print("Current player (1/2): ")
    curplayer = parse(Int, readline())
    @assert 1 <= curplayer <= 2
    return Game(board, curplayer == WHITE)
  catch
    return nothing
  end
end

# Example of creating game objects manually:
# Game(Board([10, 20], [[1 0 0 0 0 0]; [1 0 2 0 2 0]]), true)

ext: jl


text:
abstract type AbstractBoundary end

abstract type AbstractBoundaryCondition <: AbstractBoundary end
abstract type AbstractBoundaryLayer <: AbstractBoundary end

include("./BoundaryConditions/PEC.jl")
include("./BoundaryConditions/PMC.jl")
include("./BoundaryConditions/Periodic.jl")
include("./BoundaryConditions/Bloch.jl")
include("./BoundaryConditions/ABC.jl")


text:
# Floyd Warshall algorithm to find shortest paths between all pairs of vertices

function floyd_warshall!(dists::AbstractMatrix{W}) where {W} # dists: minimum distance matrix (initialized to edge distances)

    # argument checking

    n = size(dists, 1)
    if size(dists, 2) != n
        throw(ArgumentError("dists should be a square matrix."))
    end

    # initialize

    for i = 1 : n
        dists[i,i] = 0
    end

    # main loop

    for k = 1 : n, i = 1 : n, j = 1 : n
        d = dists[i,k] + dists[k,j]
        if d < dists[i,j]
            dists[i,j] = d
        end
    end

    dists
end

function floyd_warshall!(
    dists::AbstractMatrix{W},                # minimum distance matrix (initialized to edge distances)
    nexts::AbstractMatrix{Int})  where {W}   # nexts(i,j) = the next hop from i when traveling from i to j via shortest path

    # argument checking

    n = size(dists, 1)
    if size(dists, 2) != n
        throw(ArgumentError("dists should be a square matrix."))
    end

    if size(nexts) != (n, n)
        throw(ArgumentError("nexts should be an n-by-n matrix."))
    end

    # initialize

    for i = 1 : n
        dists[i,i] = 0
    end

    for j = 1 : n, i = 1 : n
        nexts[i,j] = isfinite(dists[i,j]) ? j : 0
    end

    # main loop

    for k = 1 : n, i = 1 : n, j = 1 : n
        d = dists[i,k] + dists[k,j]
        if d < dists[i,j]
            dists[i,j] = d
            nexts[i,j] = nexts[i,k]
        end
    end

    dists
end

floyd_warshall(weights::AbstractMatrix) = floyd_warshall!(copy(weights))


text:
@testset "877.stone-game.jl" begin
    @test stone_game([5, 3, 4, 5]) == true
end


text:
export generate_test_problems, generate_noisy_data!,
    generate_noisy_data, generate_clustered_noisy_data,
    generate_clustered_noisy_data!, generate_image_noisy_data,
    generate_circle, generate_ncircle, generate_image_circle

"""

This dictionary represents the list of models used in the generation of random tests.
Return the tuple `(n, model, model_str)`, where

  - `n` is the number of parameters of the model
  - `model` is the model of the form `m(x, θ)`, where `x` are the
    variables and `θ` are the parameters
  - `model_str` is the string representing the model, used to build random generated problems

"""
const model_list = Dict(
    "linear" => (2, (x, θ) -> θ[1] * x[1] + θ[2],
                 "(x, θ) -> θ[1] * x[1] + θ[2]"),
    "cubic" => (4, (x, θ) -> θ[1] * x[1]^3 + θ[2] * x[1]^2 + θ[3] * x[1] + θ[4],
                "(x, θ) -> θ[1] * x[1]^3 + θ[2] * x[1]^2 + θ[3] * x[1] + θ[4]"),
    "expon" => (3, (x, θ) -> θ[1] + θ[2] * exp(- θ[3] * x[1]),
                "(x, θ) -> θ[1] + θ[2] * exp(- θ[3] * x[1])"),
    "logistic" => (4, (x, θ) -> θ[1] + θ[2] / (1.0 + exp(- θ[3] * x[1] + θ[4])),
                   "(x, θ) -> θ[1] + θ[2] / (1.0 + exp(- θ[3] * x[1] + θ[4]))"),
    "circle" => (3, (x, θ) -> (x[1] - θ[1])^2 + (x[2] - θ[2])^2 - θ[3]^2,
                 "(x, θ) -> (x[1] - θ[1])^2 + (x[2] - θ[2])^2 - θ[3]^2"),
    "ellipse" => (6, (x, θ) -> θ[1] * x[1]^2 + θ[2] * x[1] * x[2] + θ[3] * x[2]^2 +
                  θ[4] * x[1] + θ[5] * x[2] + θ[6],
                  "(x, θ) -> θ[1] * x[1]^2 + θ[2] * x[1] * x[2] + θ[3] * x[2]^2 + θ[4] * x[1] + θ[5] * x[2] + θ[6]")
)

"""

    interval_rand!(x::Vector{Float64},
        intervals::Vector{Tuple{Float64, Float64}})

Fill a vector `x` with uniformly distributed random numbers generated
in the interval given by `intervals`. It is assumed that `length(x) ==
length(intervals)`.

Throws an `ErrorException` if the dimension of `x` is smaller the
dimension of `intervals` or if the intervals are invalid.

"""
function interval_rand!(x::Vector{Float64},
    intervals::Vector{Tuple{Float64, Float64}})

    (length(x) < length(intervals)) &&
        error("Length of vector smaller than length of intervals.")

    for v in intervals

        (v[1] > v[2]) && error("Bad interval $(v[1])>$(v[2]).")

    end

    map!((v) -> v[1] + rand() * (v[2] - v[1]), x, intervals)

end

"""

    generate_test_problems(datFilename::String, solFilename::String,
        model::Function, modelStr::String, n::Int, np::Int, p::Int;
        x_interval::Tuple{Float64, Float64}=(-10.0, 10.0),
        θSol::Vector{Float64}=10.0 * randn(n), std::Float64=200.0,
        out_times::Float64=7.0)

    generate_test_problems(datFilename::String, solFilename::String,
        model::Function, modelStr::String, n::Int, np::Int, p::Int,
        cluster_interval::Tuple{Float64, Float64};
        x_interval::Tuple{Float64, Float64}=(-10.0, 10.0),
        θSol::Vector{Float64}=10.0 * randn(n), std::Float64=200.0,
        out_times::Float64=7.0)

Generate random data files for testing fitting problems.

  - `datFilename` and `solFilename` are strings with the name of the
    files for storing the random data and solution, respectively.
  - `model` is the model function and `modelStr` is a string
    representing this model function, e.g.

         model = (x, θ) -> θ[1] * x[1] + θ[2]
         modelStr = "(x, θ) -> θ[1] * x[1] + θ[2]"

    where vector `θ` represents the parameters (to be found) of the
    model and vector `x` are the variables of the model.
  - `n` is the number of parameters
  - `np` is the number of points to be generated.
  - `p` is the number of trusted points to be used in the LOVO
    approach.

If `cluster_interval` is provided, then generates outliers only in
this interval.

Additional parameters:

  - `xMin`, `xMax`: interval for generating points in one dimensional
    tests *Deprecated*
  - `x_interval`: interval for generating points in one dimensional
    tests
  - `θSol`: true solution, used for generating perturbed points
  - `std`: standard deviation
  - `out_times`: deviation for outliers will be `out_times * std`.

"""
function generate_test_problems(datFilename::String,
    solFilename::String, model::Function, modelStr::String, n::Int,
    np::Int, p::Int; gn_kwargs...)

    # Generate data file

    θSol = nothing

    open(datFilename, "w") do data

        vdata, θSol, outliers = generate_noisy_data(model, n, np, p;
                                                    gn_kwargs...)

        # Dimension of the domain of the function to fit
        @printf(data, "%d\n", 1)

        for k = 1:np

            @printf(data, "%20.15f %20.15f %1d\n",
                    vdata[k, 1], vdata[k, 2], Int(k in outliers))

        end

    end

    # Generate solution file

    open(solFilename, "w") do sol

        println(sol, n) # number of variables

        println(sol, θSol) # parameters

        println(sol, modelStr) # function expression

    end

end

function generate_test_problems(datFilename::String,
    solFilename::String, model::Function, modelStr::String, n::Int,
    np::Int, p::Int, x_interval::Tuple{Float64, Float64},
    cluster_interval::Tuple{Float64, Float64}; gn_kwargs...)

    # Generate data file

    θSol = nothing

    open(datFilename, "w") do data

        vdata, θSol, outliers = generate_clustered_noisy_data(model,
            n, np, p, x_interval, cluster_interval; gn_kwargs...)

        # Dimension of the domain of the function to fit
        @printf(data, "%d\n", 1)

        for k = 1:np

            @printf(data, "%20.15f %20.15f %1d\n",
                    vdata[k, 1], vdata[k, 2], Int(k in outliers))

        end

    end

    # Generate solution file

    open(solFilename, "w") do sol

        println(sol, n) # number of variables

        println(sol, θSol) # parameters

        println(sol, modelStr) # function expression

    end

end

"""

    generate_circle(dat_filename::String, np::Int, p::Int;
        std::Float64=0.1, θSol::Vector{Float64}=1.0*randn(Float64, 3),
        outTimes::Float64=3.0, interval=(rand(i)*2.0*π for i = 1:np))

Generate perturbed points in a circle given by `θSol` and save to
`dat_filename` in RAFF format. Return the np x 4 matrix with data (the
4th column is 0 if the point is "correct") and a `np - p` integer
vector containing the points selected to be outliers.

  - `dat_filename` is a String with the name of the file to store
    generated data.
  - `np` is the number of points to be generated.
  - `p` is the number of *trusted points* to be used in the LOVO
    approach.

Additional configuration parameters are

  - `std`: standard deviation.
  - `θSol`: true solution, used for generating perturbed points.
  - `out_times`: deviation for outliers will be `out_times * std`.
  - `interval`: any iterable object containing `np` numbers between 0
    and 2π.

"""
function generate_circle(dat_filename::String, np::Int, p::Int;
    std::Float64=0.1, θSol::Vector{Float64}=1.0*randn(Float64, 3),
    outTimes::Float64=3.0, interval=(rand()*2.0*π for i = 1:np))

    (length(interval) != np) &&
        error("Size of interval different from given value of np")

    ρ = (α, ρ) -> [ρ * cos(α) + θSol[1], ρ * sin(α) + θSol[2]]

    data = Array{Float64, 2}(undef, np, 4)

    # Points selected to be outliers
    v = RAFF.get_unique_random_points(np, np - p)

    for (i, α) in enumerate(interval)

        pt = ρ(α, θSol[3] + std * randn())

        data[i, 3:4] .= 0.0

        # Follow the noise idea of J. Yu, H. Zheng, S. R. Kulkarni,
        # and H. V. Poor, "Two-Stage Outlier Elimination for Robust
        # Curve and Surface Fitting," EURASIP J. Adv. Signal Process.,
        # vol. 2010, no. 1, p. 154891, Dec. 2010.
        if i in v

            pt[1] += outTimes * std * randn()
            pt[2] += outTimes * std * randn()

            data[i, 4] = 1.0

        end

        data[i, 1:2] = pt

    end

    open(dat_filename, "w") do fp

        # Dimension of the domain of the function to fit
        @printf(fp, "%d\n", 2)

        for k = 1:np

            @printf(fp, "%20.15f %20.15f %20.15f %1d\n",
                    data[k, 1], data[k, 2], data[k, 3], Int(k in v))

        end

    end

    return data, v

end

"""

    generate_ncircle(dat_filename::String,np::Int, p::Int;
      std::Float64=0.1, θSol::Vector{Float64}=10.0*randn(Float64, 3),
      interval=(rand()*2.0*π for i = 1:np))

Generate perturbed points and uniform noise in a square containing the
circle given by `θSol` and save data to `dat_filename` in RAFF
format. Return the np x 4 matrix with data (the 4th column is 0 if the
point is "correct") and a `np - p` integer vector containing the
points selected to be outliers.

  - `dat_filename` is a String with the name of the file to store
    generated data.
  - `np` is the number of points to be generated.
  - `p` is the number of *trusted points* to be used in the LOVO
    approach.

Additional configuration parameters are

  - `std`: standard deviation.
  - `θSol`: true solution, used for generating perturbed points.
  - `interval`: any iterable object containing `np` numbers between 0
    and 2π.
  - `leftd`: number of times the radius of the circle that will be
    used for computing the **lower left corner** of the square for
    generation of the random noise
  - `lngth`: number of times the radius of the circle that will be
    used for computing the **side of the square** for generation of
    the random noise

"""
function generate_ncircle(dat_filename::String, np::Int, p::Int;
    std::Float64=0.1, θSol::Vector{Float64}=10.0*randn(Float64, 3),
    interval=(rand()*2.0*π for i = 1:np),leftd=2, lngth=4)

    (length(interval) != np) &&
        error("Size of interval different from given value of np")

    ρ = (α, ρ) -> [ρ * cos(α) + θSol[1], ρ * sin(α) + θSol[2]]

    data = Array{Float64, 2}(undef, np, 4)

    for (i, α) in enumerate(interval)

        pt = ρ(α, θSol[3] + std * randn())

        data[i, 1:2]  = pt
        data[i, 3:4] .= 0.0

    end

    # Just random noise

    v = Vector{Int}(undef, np - p)

    for i = p + 1:np

        data[i, 1] = θSol[1] - leftd * θSol[3] + rand() * lngth * θSol[3]
        data[i, 2] = θSol[2] - leftd * θSol[3] + rand() * lngth * θSol[3]
        data[i, 3] = 0.0
        data[i, 4] = 1.0

        v[i - p] = i

    end

    open(dat_filename, "w") do fp

        # Dimension of the domain of the function to fit
        @printf(fp, "%d\n", 2)

        for k = 1:np

            @printf(fp, "%20.15f %20.15f %20.15f %1d\n",
                    data[k, 1], data[k, 2], data[k, 3], Int(k in v))

        end

    end

    return data, v

end

"""

    generate_image_circle(dat_filename::String, w::Int, h::Int,
        np::Int, p::Int; std=0.1,
        θSol::Vector{Float64}=10.0*randn(Float64, 3),
        interval=(rand()*2.0*π for i = 1:p), thck::Int=2,
        funcsize=min(w, h))

Generate perturbed points and uniform noise in a `w`x`h` image
containing the circle given by `θSol` and save data to `dat_filename`
in RAFF format. Return the 0-1 matrix representing the black and white
image generate.

  - `dat_filename` is a String with the name of the file to store
    generated data.
  - `w` and `h` are the dimensions of the image
  - `np` is the number of points to be generated.
  - `p` is the number of *trusted points* to be used in the LOVO
    approach.

Additional configuration parameters are

  - `std`: standard deviation.
  - `θSol`: true solution, used for generating perturbed points.
  - `interval`: any iterable object containing `np` numbers between 0
    and 2π.
  - `thck`: thickness of the point in the image
  - `funcsize`: size (in pixels) that the function will use in the
    image.

"""
function generate_image_circle(dat_filename::String, w::Int, h::Int,
    np::Int, p::Int; std=0.1,
    θSol::Vector{Float64}=10.0*randn(Float64, 3),
    interval=(rand()*2.0*π for i = 1:p), thck::Int=2, funcsize=min(w, h))

    (length(interval) != p) &&
        error("Size of interval different from given value of p")

    img = zeros(Int, h, w)

    # Random uniform noise
    for k = 1:np - p

        j = 1 + Int(round(rand() * w))
        i = 1 + Int(round(rand() * h))

        img[max(1, i - thck):min(i + thck, h),
            max(1, j - thck):min(j + thck, w)] .= 1

    end

    # Circle function
    ρ = (α, ρ) -> [ρ * cos(α) + θSol[1], ρ * sin(α) + θSol[2]]

    # Store 'true' perturbed points
    x = Array{Float64}(undef, p)
    y = Array{Float64}(undef, p)

    for (i, α) in enumerate(interval)
        x[i], y[i] = ρ(α, θSol[3] + std * randn())
    end

    dx = minimum(x)
    dy = minimum(y)

    scale = funcsize / max(maximum(x) - dx, maximum(y) - dy)
    xshift = Int(round(rand() * (w - funcsize)))
    yshift = Int(round(rand() * (h - funcsize)))

    # Correct points are black

    for (tx, ty) in zip(x, y)

        j = xshift + Int(round((- dx + tx) * scale))
        i = h - yshift - Int(round((- dy + ty) * scale))

        ((i < 1) || (i > h) || (j < 1) || (j > w)) && continue

        img[max(1, i - thck):min(i + thck, h),
            max(1, j - thck):min(j + thck, w)] .= 1
    end

    open(dat_filename, "w") do fp

        # Dimension of the domain of the function to fit
        @printf(fp, "%d\n", 2)

        # Only consider the white points in the image
        for j = 1:w
            for i = 1:h

                (img[i, j] == 1) &&
                    @printf(fp, "%6d %6d %1d %1d\n",
                            j, i, 0, 0)

            end
        end

    end

    return img

end

"""

    get_unique_random_points(np::Int, npp::Int)

Choose exactly `npp` unique random points from a set containing `np`
points. This function is similar to `rand(vector)`, but does not allow
repetitions.

If `npp` < `np`, returns all the `np` points. Note that this function
is not very memory efficient, since the process of selecting unique
elements involves creating several temporary vectors.

Return a vector with the selected points.

"""
function get_unique_random_points(np::Int, npp::Int)

    # Check invalid arguments
    ((np <= 0) || (npp <=0)) && (return Vector{Int}())

    ntp = min(npp, np)

    v = Vector{Int}(undef, ntp)

    return get_unique_random_points!(v, np, npp)

end

"""

    get_unique_random_points!(v::Vector{Int}, np::Int, npp::Int)

Choose exactly `npp` unique random points from a set containing `np`
points. This function is similar to `rand(vector)`, but does not allow
repetitions.

If `npp` < `np`, returns all the `np` points. Note that this function
is not very memory efficient, since the process of selecting unique
elements involves creating several temporary vectors.

Return the vector `v` provided as argument filled with the selected
points.

"""
function get_unique_random_points!(v::Vector{Int}, np::Int, npp::Int)

    # Check invalid arguments
    ((np <= 0) || (npp <=0)) && return v

    ntp = min(np, npp)

    # Check invalid size
    (length(v) < ntp) && throw(DimensionMismatch("Incorrect size for vector."))

    # Check small np for efficiency
    if np == ntp

        for i = 1:np
            v[i] = i
        end

        return v

    end

    points = [1:np;]

    while ntp > 0

        u = rand(points, ntp)

        unique!(u)

        for i in u
            v[ntp] = i
            ntp -= 1
        end

        setdiff!(points, u)

    end

    return v

end

"""

    generate_noisy_data(model::Function, n::Int, np::Int, p::Int;
        x_interval::Tuple{Float64, Float64}=(-10.0, 10.0),
        θSol::Vector{Float64}=10.0 * randn(Float64, n),
        std::Float64=200.0, out_times::Float64=7.0)

    generate_noisy_data(model::Function, n::Int, np::Int, p::Int,
        x_interval::Tuple{Float64, Float64})

    generate_noisy_data(model::Function, n::Int, np::Int, p::Int,
        θSol::Vector{Float64}, x_interval::Tuple{Float64, Float64})

Random generate a fitting one-dimensional data problem.

This function receives a `model(x, θ)` function, the number of parameters
`n`, the number of points `np` to be generated and the number of
trusted points `p`.

If the `n`-dimensional vector `θSol` is provided, then the exact
solution will not be random generated. The interval `[xMin, xMax]`
(*deprecated*) or `x_interval` for generating the values to evaluate
`model` can also be provided.

It returns a tuple `(data, θSol, outliers)` where

  - `data`: (`np` x `3`) array, where each row contains `x` and
    `model(x, θSol)`.
  - `θSol`: `n`-dimensional vector with the exact solution.
  - `outliers`: the outliers of this data set

"""
function generate_noisy_data(model::Function, n::Int, np::Int,
    p::Int; gn_kwargs...)

    # Data matrix
    data = Array{Float64}(undef, np, 3)

    # Points selected to be outliers
    v = Vector{Int}(undef, np - p)

    return generate_noisy_data!(data, v, model, n, np, p;
                                gn_kwargs...)

end

# Deprecated

@deprecate(generate_noisy_data(model::Function, n, np, p,
    xMin::Float64, xMax::Float64), generate_noisy_data(model, n, np,
    p, (xMin, xMax)))

@deprecate(generate_noisy_data(model::Function, n::Int, np::Int,
    p::Int, θSol::Vector{Float64}, xMin::Float64, xMax::Float64),
    generate_noisy_data(model, n, np, p, θSol, (xMin, xMax)))

generate_noisy_data(model::Function, n::Int, np::Int, p::Int,
    x_interval::Tuple{Float64, Float64}) = generate_noisy_data(model,
    n, np, p; x_interval=x_interval)

generate_noisy_data(model::Function, n::Int, np::Int, p::Int,
    θSol::Vector{Float64}, x_interval::Tuple{Float64, Float64}) =
    generate_noisy_data(model, n, np, p; x_interval=x_interval,
    θSol=θSol)

"""

    generate_noisy_data!(data::AbstractArray{Float64, 2},
        v::Vector{Int}, model::Function, n::Int, np::Int, p::Int;
        x_interval::Tuple{Float64, Float64}=(-10.0, 10.0),
        θSol::Vector{Float64}=10.0 * randn(Float64, n),
        std::Float64=200.0, out_times::Float64=7.0)

Random generate a fitting one-dimensional data problem, storing the
data in matrix `data` and the outliers in vector `v`.

This function receives a `model(x, θ)` function, the number of parameters
`n`, the number of points `np` to be generated and the number of
trusted points `p`.

If the `n`-dimensional vector `θSol` is provided, then the exact
solution will not be random generated. The interval `[xMin, xMax]`
(*deprecated*) or `x_interval` for generating the values to evaluate
`model` can also be provided.

It returns a tuple `(data, θSol, outliers)` where

  - `data`: (`np` x `3`) array, where each row contains `x` and
    `model(x, θSol)`.
  - `θSol`: `n`-dimensional vector with the exact solution.
  - `outliers`: the outliers of this data set

"""
function generate_noisy_data!(data::AbstractArray{Float64, 2},
    v::Vector{Int}, model::Function, n::Int, np::Int, p::Int;
    x_interval::Tuple{Float64, Float64}=(-10.0, 10.0),
    θSol::Vector{Float64}=10.0 * randn(Float64, n),
    std::Float64=200.0, out_times::Float64=7.0)

    @assert(x_interval[1] <= x_interval[2],
            "Invalid interval for random number generation.")
    @assert(size(data) == (np, 3),
            "Invalid size of data matrix. $(size(data)) != $((np, 3)).")
    @assert(length(v) >= np - p,
            "Invalid size for vector of outliers.")

    # Generate (x_i) where x_interval[1] <= x_i <= x_interval[2] (data)
    # Fix the problem of large interval with 1 element.
    x = (np == 1) ? sum(x_interval) / 2.0 : LinRange(x_interval[1], x_interval[2], np)

    # Points selected to be outliers
    get_unique_random_points!(v, np, np - p)

    # Add noise to some random points
    sgn = sign(randn())

    for k = 1:np

        y = model(x[k], θSol) + randn() * std

        noise = 0.0

        if k in v
            y = model(x[k], θSol)
            noise = (1.0 + 2 * rand()) * out_times * std * sgn
        end

        data[k, 1] = x[k]
        data[k, 2] = y + noise
        data[k, 3] = noise

    end

    return data, θSol, v

end

"""

    generate_clustered_noisy_data(model::Function, n::Int, np::Int,
        p::Int, x_interval::Tuple{Float64,Float64},
        cluster_interval::Tuple{Float64, Float64}; kwargs...)

    generate_clustered_noisy_data(model::Function, n::Int,
        np::Int, p::Int, θSol::Vector{Float64},
        x_interval::Tuple{Float64,Float64},
        cluster_interval::Tuple{Float64, Float64}; kwargs...)

Generate a test set with clustered outliers.

The arguments and optional arguments are the same for
[`generate_noisy_data!`](@ref), with exception of tuple
`cluster_interval` which is the interval to generate the clustered
outliers.

It returns a tuple `(data, θSol, outliers)` where

  - `data`: (`np` x `3`) array, where each row contains `x` and
    `model(x, θSol)`. The same array given as argument
  - `θSol`: `n`-dimensional vector with the exact solution.
  - `outliers`: the outliers of this data set. The same vector given
    as argument.

"""
function generate_clustered_noisy_data(model::Function, n::Int,
    np::Int, p::Int, x_interval::Tuple{Float64,Float64},
    cluster_interval::Tuple{Float64, Float64}; kwargs...)

    data = Array{Float64, 2}(undef, np, 3)

    v = Vector{Int}(undef, np - p)

    return generate_clustered_noisy_data!(data, v, model, n, np, p,
               x_interval, cluster_interval; kwargs...)

end

generate_clustered_noisy_data(model::Function, n::Int, np::Int,
    p::Int, θSol::Vector{Float64}, x_interval::Tuple{Float64,Float64},
    cluster_interval::Tuple{Float64, Float64}; kwargs...) =
    generate_clustered_noisy_data(model, n, np, p, x_interval,
    cluster_interval, θSol=θSol; kwargs...)

"""

    generate_clustered_noisy_data!(data::Array{Float64, 2},
        v::Vector{Int}, model::Function, n::Int, np::Int, p::Int,
        x_interval::Tuple{Float64,Float64},
        cluster_interval::Tuple{Float64, Float64}; kwargs...)

Generate a test set with clustered outliers. This version overwrites
the content of (`np` x `3`) matrix `data` and vector `v` with integer
indices to the position of outliers in `data`.

The arguments and optional arguments are the same for
[`generate_noisy_data!`](@ref), with exception of tuple
`cluster_interval` which is the interval to generate the clustered
outliers.

It returns a tuple `(data, θSol, outliers)` where

  - `data`: (`np` x `3`) array, where each row contains `x` and
    `model(x, θSol)`. The same array given as argument
  - `θSol`: `n`-dimensional vector with the exact solution.
  - `outliers`: the outliers of this data set. The same vector given
    as argument.

"""
function generate_clustered_noisy_data!(data::Array{Float64, 2},
            v::Vector{Int}, model::Function, n::Int, np::Int,
            p::Int, x_interval::Tuple{Float64,Float64},
            cluster_interval::Tuple{Float64, Float64}; kwargs...)

    if (np - p > 0) &&
        !(x_interval[1] <= cluster_interval[1] <
          cluster_interval[2] <= x_interval[2])

        error("Bad interval for clustered data generation.")

    end

    interval_len = x_interval[2] - x_interval[1]

    fr1 = (cluster_interval[1] - x_interval[1]) / interval_len
    fr2 = (cluster_interval[2] - cluster_interval[1]) / interval_len
    fr3 = (x_interval[2] - cluster_interval[2]) / interval_len

    # Distribute data in a proportional way

    # Interval 2 will contain the clustered outliers
    np2 = max(np - p, Int(round(fr2 * np)))
    np1 = min(np - np2, Int(round(fr1 * np)))
    # Interval 3 will contain the remaining points
    np3 = max(0, np - np1 - np2)

    @debug("Clustered points: $(np1), $(np2), $(np3).")

    # Avoid repetition in the borders of the intervals

    δ_c = (cluster_interval[2] - cluster_interval[1]) / (np2 + 2)

    # Generate data

    tmpv = Vector{Int}()

    tmpd, θSol, tmpv = generate_noisy_data!(@view(data[1:np1, :]),
        tmpv, model, n, np1, np1; x_interval=(x_interval[1],
        cluster_interval[1]), kwargs...)

    generate_noisy_data!(@view(data[np1 + 1:np1 + np2, :]), v, model,
        n, np2, np2 - (np - p); x_interval=(cluster_interval[1] + δ_c,
        cluster_interval[2] - δ_c), θSol=θSol, kwargs...)

    # Update the outlier number with the correct number
    map!((x) -> x + np1, v, v)

    generate_noisy_data!(@view(data[np1 + np2 + 1:np, :]), tmpv,
        model, n, np3, np3; x_interval=(cluster_interval[2],
        x_interval[2]), θSol=θSol, kwargs...)

    return data, θSol, v

end

"""
    function generate_image_noisy_data(dat_filename::String,
    w::Int, h::Int, model::Function, n::Int, np::Int, p::Int;
    x_interval::Tuple{Number, Number}=(-10.0, 10.0),
    θSol::Vector{Float64}=10.0 * randn(Float64, n), std=2,
    thck::Int=2, funcsize=min(w, h))

Create a file `dat_filename` with data information to detect `model`
in a `w`x`h` image containing random uniform noise. **Attention**:
this function only works with `1`-dimensional models.

Return a black and white matrix representing the image.

The parameters are

  - `dat_filename`: name of the file to save data
  - `w` and `h`: dimension of the image
  - `model`: real-valued model given by a function `model(x, θ)`
  - `n`: dimension of the parameters of the model
  - `np`: number of points to be generated
  - `p`: number of trusted points that will define the correct points
    in the model

The function also accepts the following optional arguments:

  - `x_interval`: tuple representing the interval for the `x` variable
  - `θSol`: vector with the 'exact' parameters of the solution
  - `std`: error that will be added to the simulated 'correct' points
  - `thck`: thickness of the point in the image
  - `funcsize`: size (in pixels) that the function will use in the image.

"""
function generate_image_noisy_data(dat_filename::String,
    w::Int, h::Int, model::Function, n::Int, np::Int, p::Int;
    x_interval::Tuple{Number, Number}=(-10.0, 10.0),
    θSol::Vector{Float64}=10.0 * randn(Float64, n), std=2,
    thck::Int=2, funcsize=min(w, h))

    @assert(x_interval[1] <= x_interval[2],
            "Invalid interval for random number generation.")

    # Generate (x_i) where x_interval[1] <= x_i <= x_interval[2] (data)
    # Fix the problem of large interval with 1 element.
    x = (np == 1) ? sum(x_interval) / 2.0 : LinRange(x_interval[1], x_interval[2], p)

    img = zeros(Int, h, w)

    # Correct points are black
    y = map(x -> (model(x, θSol) + randn() * std), x)
    dx = minimum(x)
    dy = minimum(y)

    scale = funcsize / max(maximum(x) - dx, maximum(y) - dy)
    xshift = Int(round(rand() * (w - funcsize)))
    yshift = Int(round(rand() * (h - funcsize)))

    # Noise is gray
    for k = 1:np - p

        j = 1 + Int(round(rand() * w))
        i = 1 + Int(round(rand() * h))

        img[max(1, i - thck):min(i + thck, h),
            max(1, j - thck):min(j + thck, w)] .= 1

    end

    for (tx, ty) in zip(x, y)

        j = xshift + Int(round((- dx + tx) * scale))
        i = h - yshift - Int(round((- dy + ty) * scale))

        ((i < 1) || (i > h) || (j < 1) || (j > w)) && continue

        img[max(1, i - thck):min(i + thck, h),
            max(1, j - thck):min(j + thck, w)] .= 1
    end

    # Create data filename
    open(dat_filename, "w") do fp

        # Dimension of the domain of the function to fit
        @printf(fp, "%d\n", 2)

        # Only consider the white points in the image
        for j = 1:w
            for i = 1:h

                (img[i, j] == 1) &&
                    @printf(fp, "%6d %6d %1d %1d\n",
                            j, i, 0, 0)

            end
        end

    end

    return img

end


text:
"""
    mbgetdata(cmd0::String=""; kwargs...)

Extract bathymetry, sidescan or amplitude data from datafiles.

Parameters
----------

- **A** | **flagged** :: [Type => Number]	``Arg = value``

    Replace flagged beans with NaN. Use -A<val> to assign a constant value to the flagged beans.
- $(GMT.opt_J)
- **C** | **datatype** | **data_type** :: [Type => Number | Str | Tuple]	``Arg = 0 or "a"``

    Output SideScan, or amplitude, instead of bathymetry. This case ignores **A**
- **D** | **scaling** :: [Type => Str | Tuple]	``Arg = <mode>/<ampscale>/<ampmin>/<ampmax>``

    Sets scaling of beam amplitude or sidescan pixel values which can be applied before plotting.
- **F** | **format** :: [Type => Int]

    Sets the format for the input swath sonar data using MBIO integer format identifiers.
- **S** | **speed** :: [Type => Number]

    Sets the parameters controlng simulated illumination of bathymetry.
- **T** | **timegap** :: [Type => number]

    Sets the maximum time gap in minutes between adjacent pings before being considered a gap.
- $(GMT.opt_R)
- $(GMT.opt_V)
- $(GMT.opt_n)
- $(GMT.opt_t)
"""
function mbgetdata(cmd0::String=""; kwargs...)

	length(kwargs) == 0 && occursin(" -", cmd0) && return monolitic("mbgetdata", cmd0)

	d = init_module(false, kwargs...)[1]		# Also checks if the user wants ONLY the HELP mode

	cmd, opt_R = parse_R(d, "")
	cmd, = parse_common_opts(d, cmd, [:UVXY :params], true)

	cmd  = parse_these_opts(cmd, d, [[:A :flagged], [:D :scaling], [:F :format],
	                                 [:S :speed], [:T :timegap], [:b :star_time], [:e :end_time]])
	cmd = add_opt(d, cmd, 'C', [:C :datatype :data_type], (sidescan="_0", amplitude="_a"))

	cmd = "mbgetdata -I" * cmd0 * cmd				# In any case we need this
	finish_PS_module(d, cmd, "", true, false, false)
end


text:
function check_return_type(f::Function, M::Module, print_only_bad::Bool = false)
  color(x) = contains_any(x) ? :red : :green

  met = methods(f)

  for m in met
    if (m.func.code.module != M)
      continue
    end

    return_type_bad = true

    a = Base.return_types(f, m.sig)
    for ret in a
      return_type_bad = contains_any(ret) && return_type_bad
    end

    if (print_only_bad && return_type_bad) || !print_only_bad
      print("name:\n$(m.func.code.name)\n")
      print("signature:\n")
      print_with_color(color(m.sig), "$(m.sig)\n")
      a = Base.return_types(f, m.sig)
      print("return values:\n")
      for ret in a
        print_with_color(color(ret), "$ret ")
      end
      print("\n")
      print("location:\n")
      print("$(m.func.code.file):$(m.func.code.line)\n")
      print("--------------------------------------------------------------------------------\n")
    end
  end
end

function check_return_type(f::DataType, M::Module, print_only_bad::Bool = false)
  color(x) = contains_any(x) ? :red : :green

  met = methods(f)

  for m in met
    if (m.func.code.module != M)
      continue
    end

    return_type_bad = true

    a = Any[]
    try
      a = Base.return_types(call, m.sig)
    catch
      continue
    end

    for ret in a
      return_type_bad = contains_any(ret) && return_type_bad
    end

    if (print_only_bad && return_type_bad) || !print_only_bad
      print("name:\n$(m.func.code.name)\n")
      print("signature:\n")
      print_with_color(color(m.sig), "$(m.sig)\n")
      a = Base.return_types(f, m.sig)
      print("return values:\n")
      for ret in a
        print_with_color(color(ret), "$ret ")
      end
      print("\n")
      print("location:\n")
      print("$(m.func.code.file):$(m.func.code.line)\n")
      print("--------------------------------------------------------------------------------\n")
    end
  end
end

function check_bad_return_type(M::Module)
  for x in names(M)
    try y = eval(x)
      if isa(y, Function) || isa(y, DataType)
        check_return_type(y, M, true)
      end
    catch
      nothing
    end
  end
end

function contains_any(T)
  if typeof(T) == TypeVar
    #return contains_any(T.ub)
    return false
  elseif typeof(T) == Union
    return false
  elseif isleaftype(T)
    return false
  elseif length(T.parameters) == 0
    return T == Any
  else
    for P in T.parameters
      if contains_any(P)
        return true
      end
    end
  end
  return false
end


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:light
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.3'
#       jupytext_version: 1.0.4
#   kernelspec:
#     display_name: Julia 1.0.3
#     language: julia
#     name: julia-1.0
# ---

using SemanticModels
using Random
using LightGraphs
using MetaGraphs
using DataFrames
using GraphPlot
using DifferentialEquations
using Plots

# ## Dynamic Analysis and Dataflow Graph Construction for Type Disambiguation
#
# Our objective in this task is to ingest the main method of a program, and recursively enumerate: (1) the set of subroutines, S = {s | s is called by main} that main calls; (2) the set of variables that main interacts with (as args, or in the form of locally scoped assignment statements), V = {v | v is passed to main, or is a local variable inside of main}; (3) the set of sub-subroutines (e.g., subroutines called by each sub-routine called by main), $S_{sr}$ = {s.sr | s.sr called by s; $s \in S$}; and (4) the set of variables that the set of methods, $S$, interacts with, which may be global (w/respect to main) or local (w/respect to S) in scope, $V_{sr}$ = {s.sr.v | v is passed to an s.sr $\in S_{sr}$, or is a local variable inside of an s.sr $\in S_{sr}$}.
#
# Our hypothesis is that if we `extract method -> subroutine calls`, and `method -> variable interactions`, and use the union of these tuples and their associated metadata to construct a (directed) dataflow graph, we can use this graph to identify patterns in method call stacks, variable co-location, and different variables whose types are the same in the program as it is currently written, and as such, are candidates for type disambiguation vis-a-vis struct creation.
#
# A basic workflow that this type of formalization can support would be:
#
# Let $C$ be a class of models and $f$ be an oracle that determines $f(m) = true \iff m \in C$. Define a function `swap(m, a, b)` that takes the model $m$ and swaps all occurances of the value $a$ with the value $b$. If `f(swap(m,a,b)) = f(m)` then $a$ and $b$ are interchangeable.
#
# 1. For a given program $p_1$, identify/detect potentially ambiguous types.
# 2. For all elements of a given data type, perform pairwise swaps, and observe whether oracle still returns true.
#     - Note that we assume here that we're able to filter the set of candidate data types, to avoid (in most cases) primitive-level swaps.
# 3. Transform $p_1$ to $p_2$ by introducing structs and/or more precise data types to resolve ambiguities.
#     - Most likely a human (expert)-in-the-loop process, where we programatically generate viable stubs, and the scientist fills in these stubs with implementations.
# 4. If we iteratively apply step one to $p_2$, the ambiguities we initially detected will have been resolved, and will not be identified as candidates for type disambiguation.
#
# This workflow touches on validation-related tasks in that we should expect some swaps to result in alterations to programs that do not halt execution, but do meaningfully alter results. Ideally, we'd be able to flag inappropriate swaps of this nature by observing cases where the results represent marked deviations from previously observed values/outputs.
#
# For demonstration purposes, we begin by defining a set of extraction, graph construction, and graph visualization methods. We then provide example programs to illustrate how to apply this analytical pipeline.
#
# ### 1. Define Structs and Code Annotation Methods
# The data structures and methods defined below are intended to facilitate the annotation of our input "program" via injection (e.g., we insert data structures into the toy methods defined above for the purpose of extracting information that we can use to build the dataflow graph, including directed relationships of the types: `method -calls-> method`; `method -interacts_with-> variable`; `method -returns-> variable`; and `variable -takes-> value`).

# +
"""
Struct to hold a single method, its args, and its caller method, when recovered from parsing a program's AST.
"""
mutable struct Method{name, args, expression, called_by}
    name::name
    args::args
    expression::expression
    called_by::called_by
end

"""
Struct to hold declared methods extracted from parsing a program's AST.
"""
mutable struct DeclaredMethodsCollector{declared_methods}
    declared_methods::declared_methods
end

"""
Struct to hold called methods extracted from parsing a program's AST.
"""
mutable struct CalledMethodsCollector{called_methods}
    called_methods::called_methods
end

"""
Struct to hold a single variable with its current (updateable) value(s), along with the function performing i/o related operation with this variable.
Note that a single variable can be interact with (e.g., be passed as args to, be instantiated/modified by, and/or be returned by, multiple functions)
"""
mutable struct Variable{name, cur_val, read_write_by}
    name::name
    cur_val::cur_val
    read_write_by::read_write_by
end

"""
    set_val!(v::Variable, new_val::Any)
A method to update the "current value" of a Variable struct as a value modification (e.g., assignment or math operation) is encountered in the code.
"""
set_val!(v::Variable, new_val::Any) = (v.cur_val = new_val)

"""
Struct to hold called variables extracted from parsing and evaling a program's AST.
"""
mutable struct VariableCollector{vars}
    vars::vars
end

# +
"""
    get_all_top_level_signatures(expr_to_parse::Expr, expr_name::String, decl_methods_col::DeclaredMethodsCollector, called_methods_col::CalledMethodsCollector)
Take a Julia program (wrapped in quote), parses the AST, and recover all top level method signatures from declared and called methods.
"""
function get_all_top_level_signatures(expr_to_parse::Expr, expr_name::String, decl_methods_col::DeclaredMethodsCollector, called_methods_col::CalledMethodsCollector)

    for ex in expr_to_parse.args

        if isa(ex, LineNumberNode)
            continue
        end

        if ex.head == :function

            # the first arg of each function expression will be the function call
            signature = ex.args[1]

            # the remainder of the expression args will be the args that are passed into the function
            args = signature.args[2:length(signature.args)]

            decl_method = Method(signature, args, ex, "$expr_name")
            push!(decl_methods_col.declared_methods, decl_method)
        end

        if (ex.head == :macrocall)

            # arg 1 is the :macrocall head; arg 2 is the line number
            signature = ex.args[3]
            args = signature.args[2:length(signature.args)]
            called_method = Method(signature, args, Expr(:call, signature, args...), "$expr_name")
            push!(called_methods_col.called_methods, called_method)
        end

    end

    return decl_methods_col, called_methods_col

end
# -

"""
    collect_method_info(expr_to_parse::Expr)
Helper function that ingests an expression to parse, and outputs collector structs containing all declared and called methods
"""
function collect_method_info(expr_to_parse::Expr)

    delc_methods = DeclaredMethodsCollector(Method[])
    called_methods = CalledMethodsCollector(Method[])
    delc_methods, called_methods = get_all_top_level_signatures(expr_to_parse, "toy_expr",
        delc_methods, called_methods)
end

# +
"""
    nametype(ex::Expr)
Returns the value and type of the expression
"""
nametype(ex::Expr) = begin
    ex.head == :(::) || error("$ex is not a type assertion")
    avar = ex.args[1]
    atyp = ex.args[2]
    return avar, atyp
end

nametype(s::Symbol) = (s, :Any)

# +
"""
    collect_subroutine_calls(collector, subtree_root, caller_func)
Collects all first-level subroutine calls made by the caller_func. Note that this is (by design) not currently recursively applied.
"""

function collect_subroutine_calls(collector, subtree_root, caller_func)

    if subtree_root.head == :call

        subroutine_name = subtree_root.args[1]
        sr_args = subtree_root.args[2:end]
        subroutine = Method(subroutine_name, sr_args, subtree_root, caller_func)
        push!(collector, subroutine)
    end
end
# -

"""
    describe_args(fu::Expr, collector::Array)
A function to inject variable collection commands into the Expr representation of the program so that when the program is evaled, variables (and their values) can be collected.
"""
function describe_args(fu::Expr, collector::Array)
    f = fu
#     @show argl = f.args[1].args[2:end]
#     @show argl
#     @show body = f.args[2].args

    fname = string(f.args[1])

    for a in reverse(argl)
        avar, atyp = nametype(a)
        varname = string(avar)
        type_avar = typeof(avar)
        pushfirst!(body,
            :(println("F: ", $fname, "; ",
                $varname, " = ", $avar, " :: ", $type_avar, " <: ", $atyp))
        )

        variable_obj = Variable(varname, avar, fname)
        pushfirst!(body, :(push!($collector, $variable_obj)))

    end
    return f
end

# +
"""
annotate_program(input_expr::Expr)
Takes a program (represented as an expression), and injects annotations such that first and second-level method calls, along with variable assignment statements, can be parsed/collected when the expression is evaled.
"""
function annotate_program(input_expr::Expr)

    expr_to_annotate = copy(input_expr)
    declared, called = collect_method_info(expr_to_annotate)

    first_and_second_order_method_calls = CalledMethodsCollector([])
    program_level_variables = VariableCollector([])

    for method in declared.declared_methods

        f = method.expression

        describe_args(f, program_level_variables.vars)

        fname = string(f.args[1])
        func_obj = f.args[1]
        body = f.args[2].args

        # We want to put the top-level function call with its args into our called methods collector
        method_obj = Method(func_obj, body, func_obj, func_obj)
        pushfirst!(body, :(push!($first_and_second_order_method_calls.called_methods, $method_obj)))

        for ex in body

            if (isa(ex,LineNumberNode)) || (ex.head == :push!) || (ex.head == :println)
                continue
            end

            if ex.head == :(=)
                a = ex.args[1]
                b = ex.args[2]
                varname = string(a)

                insert!(body, length(body)-1,:(println("A: ", $fname, "; " , $varname, " = ", $b)))
                insert!(body, length(body)-1,:(push!($program_level_variables.vars, Variable($varname, $b, $fname))))

                # Collect second-level method calls
                if isa(b, Expr)
                    collect_subroutine_calls(first_and_second_order_method_calls.called_methods, b, func_obj)
                end

            end

           if ex.head in (:(*=), :(+=), :(-=), :(/=))

                a = ex.args[1]
                varname = string(a)
                b = ex.args[2]

                # The math op is the symbol right before the equals sign in the ex.head
                math_op = string(ex.head)[1]

                insert!(body, length(body)-1,:(println("A: ", $fname, "; " , $varname, " = ", $varname, " ", $math_op, " ", $b, "; == ", $a)))
                insert!(body, length(body)-1,:(push!($program_level_variables.vars, Variable($varname, $a, $fname))))

                # Collect second-level method calls
                if isa(b, Expr)
                    collect_subroutine_calls(first_and_second_order_method_calls.called_methods, b, func_obj)
                end

            end

            # The last line is a return statement; we want to capture the method -returns-> variable relationship
            if ex.head == :return && isa(body[length(body)-1],LineNumberNode)
                return_vals = ex.args[1:end]

                # can catch all assgmts take lhs and print out values of lhs
                insert!(body, length(body), :(println("R: ", $fname, " returns ", $return_vals)))

                # When we're about to return, update the value attribute of each variable
                for v in program_level_variables.vars
                    insert!(body, length(body), :(set_val!(v, v.name)))
                end
            end

            # We've reached the last line and this program doesn't contain a return statement; capture var values.
            if (ex == body[length(body)])

                if !(isa(body[length(body)-2], LineNumberNode)) && (body[length(body)-2].args[1] == :push!)

                    for v in program_level_variables.vars
                        insert!(body, length(body), :(set_val!(v, v.name)))
                    end

                    break

                end
            end
        end
    end
    return expr_to_annotate, first_and_second_order_method_calls, program_level_variables

end

# +
# Graph construction helper methods
function create_vertex_metadata_dict(vertex_struct)

    meta_data = Dict()

    for f in fieldnames(typeof(vertex_struct))
        meta_data[f] = getfield(vertex_struct, f)
    end

    return meta_data
end

function get_method_calls_method_edges(called_methods)

    methods_df = DataFrame(caller=Expr[], callee=Expr[])

    for m in called_methods
        push!(methods_df, [m.called_by, m.expression])
    end

#     @show methods_df
    return methods_df
end

function get_method_rw_variable_edges(prog_variables)

    vars_df = DataFrame(method=String[], variable=String[])

    for v in prog_variables
        push!(vars_df, [v.read_write_by, v.name])
    end

#     @show vars_df
    return vars_df
end

function add_vertices!(vertices_to_add, g, counter)

    for x in vertices_to_add
        v_name = string(x.name)

        try
            g[string(v_name), :v_name]

        catch error

            if isa(error, ErrorException)

                meta_data = create_vertex_metadata_dict(x)
                v_name = string(meta_data[:name])
                vertex_props = Dict(:v_type=> typeof(x), :meta_data=> meta_data)

                add_vertex!(g)
                set_indexing_prop!(g, counter, :v_name, v_name)
                set_props!(g, counter, vertex_props)

                counter += 1
            end
        end

    end

    return g, counter

end

# need to be able to make this conditional on the array type (Eg methods vs variables)

function add_directed_edges!(g, edge_collector::CalledMethodsCollector, edge_type::String)

    for method_call in edge_collector.called_methods

        src_id = g[string(method_call.called_by.args[1]), :v_name]
        dst_id = g[string(method_call.name), :v_name]

        # Check if \e \in G; if Y, increment counter
        if has_edge(g, src_id, dst_id)
            cur_freq = get_prop(g, Edge(src_id, dst_id), :weight)
            updated_props = Dict(:weight => cur_freq + 1, :type => edge_type)
            set_props!(g, Edge(src_id, dst_id), updated_props)
        else

            e_props = Dict(:weight => 1, :type => edge_type)
            add_edge!(g, Edge(src_id, dst_id))
            set_props!(g, Edge(src_id, dst_id), e_props)
        end
    end
end

function add_directed_edges!(g, edge_collector::VariableCollector, edge_type::String)

    for var in edge_collector.vars

        src_name = split(string(var.read_write_by), "(")[1]
        src_id = g["$src_name", :v_name]
        dst_id = g[string(var.name), :v_name]

        # Check if \e \in G; if Y, increment counter
        if has_edge(g, src_id, dst_id)
            cur_freq = get_prop(g, Edge(src_id, dst_id), :weight)
            updated_props = Dict(:weight => cur_freq + 1, :type => edge_type)
            set_props!(g, Edge(src_id, dst_id), updated_props)
        else

            e_props = Dict(:weight => 1, :type => edge_type)
            add_edge!(g, Edge(src_id, dst_id))
            set_props!(g, Edge(src_id, dst_id), e_props)
        end
    end
end

function create_dataflow_graph(called_methods, prog_vars, m_df, v_df)

    g = MetaDiGraph()

    set_indexing_prop!(g, :v_name)

    # THIS HAS TO START AT ONE, NOT ZERO.
    counter = 1

    # Add the method vertices representing called methods
    g, counter = add_vertices!(called_methods.called_methods, g, counter)

    # Add the caller methods
    callers = [v.called_by for v in called_methods.called_methods]
    caller_methods = [Method(c.args[1], c.args[2:end], c, :nothing) for c in callers]

    g, counter = add_vertices!(caller_methods, g, counter)

    # Add the variable vertices
    g, counter = add_vertices!(prog_vars.vars, g, counter)

    # Add the method -calls-> method edges
    add_directed_edges!(g, called_methods, "calls")

    add_directed_edges!(g, prog_vars, "interacts_with")

    return g

end

# +
# Graph visualization helper method
function plot_dataflow_graph(g)
    println(ne(g))

    edge_labels = [(get_prop(g,e, :type),get_prop(g,e, :weight))  for e in edges(g)]

    vertex_labels = []

    for v in vertices(g)

        if get_prop(g,v, :v_type) <: Variable
            v_name = get_prop(g,v, :v_name)
            v_type = typeof(get_prop(g,v, :meta_data)[:cur_val])
            push!(vertex_labels, (string(v_name, "::", v_type)))
        else
            push!(vertex_labels, (get_prop(g,v, :v_name)))
        end
    end

    node_sizes = [LightGraphs.outdegree(g, v) + 3 for v in vertices(g)]

    # Our nodes either represent methods or variables
    node_colors = [get_prop(g,v, :v_type) <: Variable ? "#6682E0" : "#D5635C"  for v in vertices(g)]
    node_strokes = [get_prop(g,v, :v_type) <: Variable ? "#4063D8" : "#CB3C33"  for v in vertices(g)]

    node_label_sizes = [LightGraphs.outdegree(g, v) for v in vertices(g)]

    layout=(args...)->spring_layout(args...; C=40)

    gplot(g,
        edgelabel=edge_labels,
        nodelabel=vertex_labels,
        nodesize=node_sizes,
        edgelabelsize=2.0,
        nodelabelsize=2.5,
        #layout=stressmajorize_layout,
        layout=layout,
        nodefillc=node_colors,
        nodestrokec= node_strokes,
        nodestrokelw = 1.0,
        nodelabeldist=1.5,
        nodelabelangleoffset=π/4,
        NODELABELSIZE = 3.5,
        EDGELABELSIZE = 3.5)

end

# -

# ### 2. Define Example Programs
# To begin, we can define some toy methods with useful properties (e.g., sub-routine calls; recursive calls, empty args). We wrap our method declarations and calls in a `quote` so that we can parse the program (and its component pieces) as `Expr` objects; this allows us to traverse the ASTs.

sir_expr = quote
function sir_ode(du, u, p, t)
    #Infected per-Capita Rate
    β = p[1]
    #Recover per-capita rate
    γ = p[2]
    #Susceptible Individuals
    S = u[1]
    #Infected by Infected Individuals
    I = u[2]

    du[1] = -β * S * I
    du[2] = β * S * I - γ * I
    du[3] = γ * I
end
function sir_ode2(du,u,p,t)
    S,I,R = u
    b,g = p
    du[1] = -b*S*I
    du[2] = b*S*I-g*I
    du[3] = g*I
end
parms = [0.1,0.05]
init = [0.99,0.01,0.0]
tspan = (0.0,200.0)
sir_prob2 = ODEProblem(sir_ode2,init,tspan,parms)
@show sir_sol = solve(sir_prob2,saveat = 0.1)
end

# + {"endofcell": "--"}
macro_parasite_expr = quote
function main()
# # +
function macroParasiteModelFunction(dY,Y,p,t)
    #Host Birth Rate
    a = p[1]
    #Parasite Influence on host birth rate
    b = p[2]
    #Parasite induced host mortality
    α = p[3]
    #Parasite induced decrease in host reproduction
    β = p[4]
    #Intrinsic death rate of parasites
    μ = p[5]
    #dispersion aggregation parameter
    k = p[6]
    #Rate of production of new free-living stages
    λ = p[7]
    #Death Rate of Free-Living Stages
    γ = p[8]

    #Host Population
    H = Y[1]
    #Parasite Population
    P = Y[2]
    #Infective Stages
    W = Y[3]

dY[1] = (a-b)*H - α*P
dY[2] = β*H*W - (μ + α + b) * P - (α*((P^2)/H)*((k+1)/k))
dY[3] = λ*P - (γ*W) - (β*H*W)
end

# -

par=[1.4,1.05,0.0003,0.01,0.5,0.1,10.0,10.0]
init=[100.0,10.0,10.0]
tspan=(0.0,100.0)

macro_odeProblem = ODEProblem(macroParasiteModelFunction,init,tspan,par)

sol=solve(macro_odeProblem);
@show typeof(sol)

p = plot(sol,xlabel="Time",yscale=:log10)
end
@show main()
end
# --

# +
# Let's define some toy functions with useful properties (e.g., sub-routine calls; recursive calls, empty args):
example_expr = quote

function fib(n::Int64)
    if n <= 2
        return 1
    else
        return fib(n-1) + fib(n-2)
    end
end

function circumference(r)
    circum = 2*pi*r
    return circum
end

function exp_decay(init_value, decay_rate, t)
    y = init_value(1-decay_rate)^t
    return y
end

function compute_circle_metrics(r)
    c = circumference(r)
    a = pi*r^2
    return c,a
end

function add_some_noise(x)
    rng = MersenneTwister(1234);
    x *= Base.randn(rng,1)
    return x
end

function combine_fib_and_circle(x::Int64)
    fib_res = fib(x)
    circum, area = compute_circle_metrics(fib_res)
    return circum, area
end

@show combine_fib_and_circle(7)
end

# +
# Let's define some toy functions with useful properties (e.g., sub-routine calls; recursive calls, empty args):
example_expr_fib = quote

function fib(n::Int64)
    if n <= 2
        return 1
    else
        return fib(n-1) + fib(n-2)
    end
end
function main(n::Int64)
    out = fib(n)
    return out
end

@show main(5)
end
# -

# ### 3. Use code annotation script to collect metadata to build dataflow graph

# +
# First, let's run our annotation method on the example expression.
# When we do this, we'll get back a revised (annotated) expression that can be evaluated to yield the information
# we can use to build the dataflow graph.

annotated_expr, first_and_second_order_calls, prog_vars = annotate_program(sir_expr)
#annotated_expr, first_and_second_order_calls, prog_vars = annotate_program(macro_parasite_expr)
#annotated_expr, first_and_second_order_calls, prog_vars = annotate_program(example_expr)
# -

eval(annotated_expr)

methods_df = get_method_calls_method_edges(first_and_second_order_calls.called_methods)
vars_df = get_method_rw_variable_edges(prog_vars.vars)
g_out = create_dataflow_graph(first_and_second_order_calls, prog_vars, methods_df, vars_df)

plot_dataflow_graph(g_out)

# ## Conclusions
#
# Dynamic program analysis is able to capture variables, values, functions, and types. Knowledge graphs extracted from scientific model programs include rich information that can be tied to the concept graphs from the text, through the intermediary of variables.
#
# We can collect this dynamic analysis information from programs that use highly optimized complex libraries such as DifferentialEquations.jl, which we did not write. This shows promise for the ability of model metaprogramming and semantically annotated syntax trees to provide information about the structure and function of scientific models.


text:
# In theory, for load-based MMF routing, using any monotone function of the load
# should give the same result. We check that fact. We proceed as follows:
# for each topology, for several objective functions (including the simple
# load), compute a routing; compare all those routings.

# For simplicity, we will focus on load functions that can be solved with
# a standard SOCP solver (like CPLEX, Gurobi, Mosek).

# Theoretical reference:
# https://ieeexplore.ieee.org/document/4738463 : Section 4.1.2.

include("topologies.jl")

using PrettyTables, LightGraphs, DelimitedFiles
using JuMP, CPLEX, Gurobi, Mosek, MosekTools

# opt = CPLEX.Optimizer
opt = optimizer_with_attributes(Gurobi.Optimizer, "NonConvex" => 2)
# opt = Mosek.Optimizer

ftopologies = [t for t in topologies if ! (t.name in ["att", "abilene"])]
ε = 1.0e-8

objs = [Load(), KleinrockLoad(true), FortzThorupLoad()]#, AlphaFairness(0.5)]#, AlphaFairness(1.0)]#, AlphaFairness(1.5), AlphaFairness(2.0)]
table = Matrix{Any}(undef, 2 * length(ftopologies), 2 + length(objs))
csv_cols = ["Instance", "Number of edges", "Number of nodes", "Number of demands",
            "Objective function for MMF", "MMF iterations", "Time to create the master problem [ms]",
            "Computation time (per iteration) [ms]", "Maximum iteration time [ms]", "Total iteration time [ms]",
            "Total time [ms]", "Maximum load of the obtained routing", "MMF status code"]
csv = Matrix{String}(undef, 1 + length(ftopologies) * length(objs), length(csv_cols))
csv[1, :] = csv_cols
csv_row = 2

for (i, t) in enumerate(ftopologies)
    println("=> $(t.name)")
    g, k, dm = topology_to_graphs(t)

    loads = Any[]
    n_paths = Any[]

    for obj in objs
        println("  -> $(obj)")

        if typeof(obj) != AlphaFairness
            mt = ModelType(obj, MinMaxFair(ε))
        else
            mt = ModelType(obj, MaxMinFair(ε))
        end

        rd = RoutingData(g, k, opt, mt, name=t.name, traffic_matrix=dm)
        sol = compute_routing(rd)

        r = sol.routings[sol.n_iter]
        max_load = compute_max_load(r)
        push!(loads, max_load)
        push!(n_paths, length(r.paths))

        global csv_row
        csv[csv_row, :] = String[t.name, string(ne(g)), string(nv(g)), string(ne(k)),
                                 string(obj), string(sol.n_iter), string(sol.time_create_master_model_ms),
                                 string(sol.time_solve_ms), string(maximum(values(sol.time_solve_ms))), string(sol.total_time_solve_ms),
                                 string(sol.total_time_ms), string(max_load), string(sol.result)]
        csv_row += 1
    end

    table[2 * i - 1, :] = vcat([t.name, "Max load"], loads)
    table[2 * i, :] = vcat(["", "Nb. paths"], n_paths)
end

@ptconfclean
@ptconf tf = unicode_rounded
@pt :header = ["", "", objs...] table

# println(csv)
writedlm("$(@__DIR__)/mmf_compare_objectives.csv", csv, ',')


text:
# TODO: Implement liblazy functions over standard arrays / other collections

walk(inner, outer, xs::Array) = @>> xs map(inner) outer

function Base.split{T}(xs::Vector{T}, x; keep = true)
  result = Vector{T}[]
  push!(result, T[])
  for i = 1:length(xs)
    if xs[i] == x
      (keep || !(isempty(result[end]) || i == length(xs))) &&
        push!(result, T[])
    else
      push!(result[end], xs[i])
    end
  end
  isempty(result[end]) && pop!(result)
  return result
end

export frequencies

function frequencies(xs)
  freqs = Dict{eltype(xs),Int}()
  for x in xs
    freqs[x] = get(freqs, x, 0) + 1
  end
  return freqs
end

export init

init(xs) = xs[1:end-1]

# function Base.getindex{K, V}(d::Dict{K, V}, ks::Vector{K})
#   vs = similar(ks, V, 0)
#   for k in ks
#     push!(vs, d[k])
#   end
#   return vs
# end

for f in (:takewhile, :splitby)
  @eval $(f)(f::Function, xs) = $(f)(f, seq(xs))
end

Base.get!(d, k, v) = (d[k] = get(d, k, v))


text:
__precompile__()

module SymEngine

import Base: show, convert, real, imag
import Compat: String, unsafe_string, @compat, denominator, numerator, invokelatest, Cvoid, Nothing, MathConstants.γ, MathConstants.e, MathConstants.φ, MathConstants.catalan, LinearAlgebra, finalizer, Libdl

export Basic, symbols, @vars, @funs, SymFunction
export free_symbols, get_args
export ascii_art
export subs, lambdify, N
export series

const deps_file = joinpath(dirname(@__FILE__), "..", "deps", "deps.jl")
const deps_in_file = joinpath(dirname(@__FILE__), "..", "deps", "deps.jl.in")

noinit = false

if !isfile(deps_file)
    cp(deps_in_file, deps_file)
end

include(deps_file)
include("utils.jl")
const have_mpfr = have_component("mpfr")
const have_mpc = have_component("mpc")
const libversion = get_libversion()

if VERSION > VersionNumber("0.7.0-DEV")
    _finalizer(f, o) = finalizer(f, o)
else
    _finalizer(f, o) = finalizer(o, f)
end

include("types.jl")
include("ctypes.jl")
include("display.jl")
include("mathops.jl")
include("mathfuns.jl")
include("subs.jl")
include("numerics.jl")
include("simplify.jl")
include("calculus.jl")
include("recipes.jl")
include("dense-matrix.jl")

if !noinit
    __init__() = init_constants()
end

end


text:
show_construction(sys::LTISystem; kwargs...) = show_construction(stdout, sys; kwargs...)

"""
    show_construction([io::IO,] sys::LTISystem; name = "temp", letb = true)

Print code to `io` that reconstructs `sys`.
- `letb`: If true, the code is surrounded by a let block.

```@example
julia> sys = ss(tf(1, [1, 1]))
StateSpace{Continuous, Float64}
A =
 -1.0
B =
 1.0
C =
 1.0
D =
 0.0

Continuous-time state-space model

julia> show_construction(sys, name="Jörgen")
sys = let
    JörgenA = [-1.0;;]
    JörgenB = [1.0;;]
    JörgenC = [1.0;;]
    JörgenD = [0.0;;]
    ss(JörgenA, JörgenB, JörgenC, JörgenD)
end
```
"""
function show_construction(io::IO, sys::LTISystem; name = "temp", letb = true)
    # sys = StateSpace(sys)
    letb && println(io, "sys = let")
    prestr = letb ? "    " : ""
    println(io, prestr*"$(name)A = ", sys.A)
    println(io, prestr*"$(name)B = ", sys.B)
    println(io, prestr*"$(name)C = ", sys.C)
    println(io, prestr*"$(name)D = ", sys.D)
    letb || print(io, "sys = ")
    if isdiscrete(sys)
        println(io, prestr*"ss($(name)A, $(name)B, $(name)C, $(name)D, $(sys.Ts))")
    else
        println(io, prestr*"ss($(name)A, $(name)B, $(name)C, $(name)D)")
    end
    letb && println(io, "end")
    nothing
end

function Base.vec(sys::LTISystem)
    [vec(sys.A); vec(sys.B); vec(sys.C); vec(sys.D)]
end

"""
    vec2sys(v::AbstractArray, ny::Int, nu::Int, ts = nothing)

Create a statespace system from the parameters
```julia
v = vec(sys) = [vec(sys.A); vec(sys.B); vec(sys.C); vec(sys.D)]
```
Use `vec(sys)` to create `v`.

This can be useful in order to convert to and from vectors for, e.g., optimization.

```@example
julia> sys  = ss(tf(1, [1, 1]))
StateSpace{Continuous, Float64}
A =
 -1.0
B =
 1.0
C =
 1.0
D =
 0.0

Continuous-time state-space model

julia> v    = vec(sys)
4-element Vector{Float64}:
 -1.0
  1.0
  1.0
  0.0

julia> sys2 = vec2sys(v, sys.ny, sys.nu)
StateSpace{Continuous, Float64}
A =
 -1.0
B =
 1.0
C =
 1.0
D =
 0.0

Continuous-time state-space model
```
"""
function vec2sys(v::AbstractArray, ny::Int, nu::Int, ts=nothing)
    n = length(v)
    p = (ny+nu)
    nx = Int(-p/2 + sqrt(p^2 - 4nu*ny + 4n)/2)
    @assert n == nx^2 + nx*nu + ny*nx + ny*nu
    ai = (1:nx^2)
    bi = (1:nx*nu) .+ ai[end]
    ci = (1:nx*ny) .+ bi[end]
    di = (1:nu*ny) .+ ci[end]
    A = reshape(v[ai], nx, nx)
    B = reshape(v[bi], nx, nu)
    C = reshape(v[ci], ny, nx)
    D = reshape(v[di], ny, nu)
    ts === nothing ? ss(A, B, C, D) : ss(A, B, C, D, ts)
end

# sys2vec = @(sys) [
#         size(sys.A,1);
#         size(sys.B,2);
#         size(sys.C,1);
#         sys.A(:);
#         sys.B(:);
#         sys.C(:);
#         sys.D(:)
#     ]
# function vec2sys(v::AbstractArray, ts = nothing)
#     nx = Int(v[1])
#     nu = Int(v[2])
#     ny = Int(v[3])
#     ai = (1:nx^2) .+ 3
#     bi = (1:nx*nu) .+ ai[end]
#     ci = (1:nx*ny) .+ bi[end]
#     di = (1:nu*ny) .+ ci[end]
#     A = reshape(v[ai], nx, nx)
#     B = reshape(v[bi], nx, nu)
#     C = reshape(v[ci], ny, nx)
#     D = reshape(v[di], ny, nu)
#     ts === nothing ? ss(A, B, C, D) : ss(A, B, C, D, ts)
# end


text:
module MimiPAGE2009

using Mimi

export getpage

include("utils/load_parameters.jl")
include("utils/mctools.jl")

include("mcs.jl")

include("compute_scc.jl")

include("components/CO2emissions.jl")
include("components/CO2cycle.jl")
include("components/CO2forcing.jl")
include("components/CH4emissions.jl")
include("components/CH4cycle.jl")
include("components/CH4forcing.jl")
include("components/N2Oemissions.jl")
include("components/N2Ocycle.jl")
include("components/N2Oforcing.jl")
include("components/LGemissions.jl")
include("components/LGcycle.jl")
include("components/LGforcing.jl")
include("components/SulphateForcing.jl")
include("components/TotalForcing.jl")
include("components/ClimateTemperature.jl")
include("components/SeaLevelRise.jl")
include("components/GDP.jl")
include("components/MarketDamages.jl")
include("components/NonMarketDamages.jl")
include("components/Discontinuity.jl")
include("components/AdaptationCosts.jl")
include("components/SLRDamages.jl")
include("components/AbatementCostParameters.jl")
include("components/AbatementCosts.jl")
include("components/TotalAbatementCosts.jl")
include("components/TotalAdaptationCosts.jl")
include("components/Population.jl")
include("components/EquityWeighting.jl")

function buildpage(m::Model, policy::String="policy-a")

    #add all the components
    add_comp!(m, co2emissions)
    add_comp!(m, co2cycle)
    add_comp!(m, co2forcing)
    add_comp!(m, ch4emissions)
    add_comp!(m, ch4cycle)
    add_comp!(m, ch4forcing)
    add_comp!(m, n2oemissions)
    add_comp!(m, n2ocycle)
    add_comp!(m, n2oforcing)
    add_comp!(m, LGemissions)
    add_comp!(m, LGcycle)
    add_comp!(m, LGforcing)
    add_comp!(m, SulphateForcing)
    add_comp!(m, TotalForcing)
    add_comp!(m, ClimateTemperature)
    add_comp!(m, SeaLevelRise)

    #Socio-Economics
    population = addpopulation(m)
    add_comp!(m, GDP)

    #Abatement Costs
    abatementcostparameters_CO2 = addabatementcostparameters(m, :CO2, policy)
    abatementcostparameters_CH4 = addabatementcostparameters(m, :CH4, policy)
    abatementcostparameters_N2O = addabatementcostparameters(m, :N2O, policy)
    abatementcostparameters_Lin = addabatementcostparameters(m, :Lin, policy)

    abatementcosts_CO2 = addabatementcosts(m, :CO2, policy)
    abatementcosts_CH4 = addabatementcosts(m, :CH4, policy)
    abatementcosts_N2O = addabatementcosts(m, :N2O, policy)
    abatementcosts_Lin = addabatementcosts(m, :Lin, policy)
    add_comp!(m, TotalAbatementCosts)

    #Adaptation Costs
    adaptationcosts_sealevel = addadaptationcosts_sealevel(m)
    adaptationcosts_economic = addadaptationcosts_economic(m)
    adaptationcosts_noneconomic = addadaptationcosts_noneconomic(m)
    add_comp!(m, TotalAdaptationCosts)

    # Impacts
    slrdamages = addslrdamages(m)
    marketdamages = addmarketdamages(m)
    nonmarketdamages = addnonmarketdamages(m)
    add_comp!(m, Discontinuity)

    #Equity weighting and Total Costs
    add_comp!(m, EquityWeighting)

    #connect parameters together
    connect_param!(m, :co2cycle => :e_globalCO2emissions, :co2emissions => :e_globalCO2emissions)
    connect_param!(m, :co2cycle => :rt_g0_baseglobaltemp, :ClimateTemperature => :rt_g0_baseglobaltemp)
    connect_param!(m, :co2cycle => :rt_g_globaltemperature, :ClimateTemperature => :rt_g_globaltemperature)

    connect_param!(m, :co2forcing => :c_CO2concentration, :co2cycle => :c_CO2concentration)

    connect_param!(m, :ch4cycle => :e_globalCH4emissions, :ch4emissions => :e_globalCH4emissions)
    connect_param!(m, :ch4cycle => :rtl_g0_baselandtemp, :ClimateTemperature => :rtl_g0_baselandtemp)
    connect_param!(m, :ch4cycle => :rtl_g_landtemperature, :ClimateTemperature => :rtl_g_landtemperature)

    connect_param!(m, :ch4forcing => :c_CH4concentration, :ch4cycle => :c_CH4concentration)
    connect_param!(m, :ch4forcing => :c_N2Oconcentration, :n2ocycle => :c_N2Oconcentration)

    connect_param!(m, :n2ocycle => :e_globalN2Oemissions, :n2oemissions => :e_globalN2Oemissions)
    connect_param!(m, :n2ocycle => :rtl_g0_baselandtemp, :ClimateTemperature => :rtl_g0_baselandtemp)
    connect_param!(m, :n2ocycle => :rtl_g_landtemperature, :ClimateTemperature => :rtl_g_landtemperature)

    connect_param!(m, :n2oforcing => :c_CH4concentration, :ch4cycle => :c_CH4concentration)
    connect_param!(m, :n2oforcing => :c_N2Oconcentration, :n2ocycle => :c_N2Oconcentration)

    connect_param!(m, :LGcycle => :e_globalLGemissions, :LGemissions => :e_globalLGemissions)
    connect_param!(m, :LGcycle => :rtl_g0_baselandtemp, :ClimateTemperature => :rtl_g0_baselandtemp)
    connect_param!(m, :LGcycle => :rtl_g_landtemperature, :ClimateTemperature => :rtl_g_landtemperature)

    connect_param!(m, :LGforcing => :c_LGconcentration, :LGcycle => :c_LGconcentration)

    connect_param!(m, :TotalForcing => :f_CO2forcing, :co2forcing => :f_CO2forcing)
    connect_param!(m, :TotalForcing => :f_CH4forcing, :ch4forcing => :f_CH4forcing)
    connect_param!(m, :TotalForcing => :f_N2Oforcing, :n2oforcing => :f_N2Oforcing)
    connect_param!(m, :TotalForcing => :f_lineargasforcing, :LGforcing => :f_LGforcing)

    connect_param!(m, :ClimateTemperature => :ft_totalforcing, :TotalForcing => :ft_totalforcing)
    connect_param!(m, :ClimateTemperature => :fs_sulfateforcing, :SulphateForcing => :fs_sulphateforcing)

    connect_param!(m, :SeaLevelRise => :rt_g_globaltemperature, :ClimateTemperature => :rt_g_globaltemperature)

    connect_param!(m, :GDP => :pop_population, :Population => :pop_population)

    for allabatement in [
        (:AbatementCostParametersCO2, :AbatementCostsCO2),
        (:AbatementCostParametersCH4, :AbatementCostsCH4),
        (:AbatementCostParametersN2O, :AbatementCostsN2O),
        (:AbatementCostParametersLin, :AbatementCostsLin)]

        abatementcostparameters, abatementcosts = allabatement

        connect_param!(m, abatementcostparameters => :yagg, :GDP => :yagg_periodspan)
        connect_param!(m, abatementcostparameters => :cbe_absoluteemissionreductions, abatementcosts => :cbe_absoluteemissionreductions)

        connect_param!(m, abatementcosts => :zc_zerocostemissions, abatementcostparameters => :zc_zerocostemissions)
        connect_param!(m, abatementcosts => :q0_absolutecutbacksatnegativecost, abatementcostparameters => :q0_absolutecutbacksatnegativecost)
        connect_param!(m, abatementcosts => :blo, abatementcostparameters => :blo)
        connect_param!(m, abatementcosts => :alo, abatementcostparameters => :alo)
        connect_param!(m, abatementcosts => :bhi, abatementcostparameters => :bhi)
        connect_param!(m, abatementcosts => :ahi, abatementcostparameters => :ahi)

    end

    connect_param!(m, :TotalAbatementCosts => :tc_totalcosts_co2, :AbatementCostsCO2 => :tc_totalcost)
    connect_param!(m, :TotalAbatementCosts => :tc_totalcosts_n2o, :AbatementCostsN2O => :tc_totalcost)
    connect_param!(m, :TotalAbatementCosts => :tc_totalcosts_ch4, :AbatementCostsCH4 => :tc_totalcost)
    connect_param!(m, :TotalAbatementCosts => :tc_totalcosts_linear, :AbatementCostsLin => :tc_totalcost)
    connect_param!(m, :TotalAbatementCosts => :pop_population, :Population => :pop_population)

    connect_param!(m, :AdaptiveCostsEconomic => :gdp, :GDP => :gdp)
    connect_param!(m, :AdaptiveCostsNonEconomic => :gdp, :GDP => :gdp)
    connect_param!(m, :AdaptiveCostsSeaLevel => :gdp, :GDP => :gdp)

    connect_param!(m, :TotalAdaptationCosts => :ac_adaptationcosts_economic, :AdaptiveCostsEconomic => :ac_adaptivecosts)
    connect_param!(m, :TotalAdaptationCosts => :ac_adaptationcosts_noneconomic, :AdaptiveCostsNonEconomic => :ac_adaptivecosts)
    connect_param!(m, :TotalAdaptationCosts => :ac_adaptationcosts_sealevelrise, :AdaptiveCostsSeaLevel => :ac_adaptivecosts)
    connect_param!(m, :TotalAdaptationCosts => :pop_population, :Population => :pop_population)

    connect_param!(m, :SLRDamages => :s_sealevel, :SeaLevelRise => :s_sealevel)
    connect_param!(m, :SLRDamages => :cons_percap_consumption, :GDP => :cons_percap_consumption)
    connect_param!(m, :SLRDamages => :tct_per_cap_totalcostspercap, :TotalAbatementCosts => :tct_per_cap_totalcostspercap)
    connect_param!(m, :SLRDamages => :act_percap_adaptationcosts, :TotalAdaptationCosts => :act_percap_adaptationcosts)
    connect_param!(m, :SLRDamages => :atl_adjustedtolerablelevelofsealevelrise, :AdaptiveCostsSeaLevel => :atl_adjustedtolerablelevel, ignoreunits=true)
    connect_param!(m, :SLRDamages => :imp_actualreductionSLR, :AdaptiveCostsSeaLevel => :imp_adaptedimpacts)
    connect_param!(m, :SLRDamages => :isatg_impactfxnsaturation, :GDP => :isatg_impactfxnsaturation)

    connect_param!(m, :MarketDamages => :rtl_realizedtemperature, :ClimateTemperature => :rtl_realizedtemperature)
    connect_param!(m, :MarketDamages => :rgdp_per_cap_SLRRemainGDP, :SLRDamages => :rgdp_per_cap_SLRRemainGDP)
    connect_param!(m, :MarketDamages => :rcons_per_cap_SLRRemainConsumption, :SLRDamages => :rcons_per_cap_SLRRemainConsumption)
    connect_param!(m, :MarketDamages => :atl_adjustedtolerableleveloftemprise, :AdaptiveCostsEconomic => :atl_adjustedtolerablelevel, ignoreunits=true)
    connect_param!(m, :MarketDamages => :imp_actualreduction, :AdaptiveCostsEconomic => :imp_adaptedimpacts)
    connect_param!(m, :MarketDamages => :isatg_impactfxnsaturation, :GDP => :isatg_impactfxnsaturation)

    connect_param!(m, :NonMarketDamages => :rtl_realizedtemperature, :ClimateTemperature => :rtl_realizedtemperature)
    connect_param!(m, :NonMarketDamages => :rgdp_per_cap_MarketRemainGDP, :MarketDamages => :rgdp_per_cap_MarketRemainGDP)
    connect_param!(m, :NonMarketDamages => :rcons_per_cap_MarketRemainConsumption, :MarketDamages => :rcons_per_cap_MarketRemainConsumption)
    connect_param!(m, :NonMarketDamages =>:atl_adjustedtolerableleveloftemprise, :AdaptiveCostsNonEconomic =>:atl_adjustedtolerablelevel, ignoreunits=true)
    connect_param!(m, :NonMarketDamages => :imp_actualreduction, :AdaptiveCostsNonEconomic => :imp_adaptedimpacts)
    connect_param!(m, :NonMarketDamages => :isatg_impactfxnsaturation, :GDP => :isatg_impactfxnsaturation)

    connect_param!(m, :Discontinuity => :rgdp_per_cap_NonMarketRemainGDP, :NonMarketDamages => :rgdp_per_cap_NonMarketRemainGDP)
    connect_param!(m, :Discontinuity => :rt_g_globaltemperature, :ClimateTemperature => :rt_g_globaltemperature)
    connect_param!(m, :Discontinuity => :rgdp_per_cap_NonMarketRemainGDP, :NonMarketDamages => :rgdp_per_cap_NonMarketRemainGDP)
    connect_param!(m, :Discontinuity => :rcons_per_cap_NonMarketRemainConsumption, :NonMarketDamages => :rcons_per_cap_NonMarketRemainConsumption)
    connect_param!(m, :Discontinuity => :isatg_saturationmodification, :GDP => :isatg_impactfxnsaturation)

    connect_param!(m, :EquityWeighting => :pop_population, :Population => :pop_population)
    connect_param!(m, :EquityWeighting => :tct_percap_totalcosts_total, :TotalAbatementCosts => :tct_per_cap_totalcostspercap)
    connect_param!(m, :EquityWeighting => :act_adaptationcosts_total, :TotalAdaptationCosts => :act_adaptationcosts_total)
    connect_param!(m, :EquityWeighting => :act_percap_adaptationcosts, :TotalAdaptationCosts => :act_percap_adaptationcosts)
    connect_param!(m, :EquityWeighting => :cons_percap_consumption, :GDP => :cons_percap_consumption)
    connect_param!(m, :EquityWeighting => :cons_percap_consumption_0, :GDP => :cons_percap_consumption_0)
    connect_param!(m, :EquityWeighting => :cons_percap_aftercosts, :SLRDamages => :cons_percap_aftercosts)
    connect_param!(m, :EquityWeighting => :rcons_percap_dis, :Discontinuity => :rcons_per_cap_DiscRemainConsumption)
    connect_param!(m, :EquityWeighting => :yagg_periodspan, :GDP => :yagg_periodspan)

    return m
end

function initpage(m::Model, policy::String="policy-a")
    p = load_parameters(m, policy)
    p["y_year_0"] = 2008.
    p["y_year"] = Mimi.dim_keys(m.md, :time)
    set_leftover_params!(m, p)
end

function get_model(policy::String="policy-a")
    m = Model()
    set_dimension!(m, :time, [2009, 2010, 2020, 2030, 2040, 2050, 2075, 2100, 2150, 2200])
    set_dimension!(m, :region, ["EU", "USA", "OECD","USSR","China","SEAsia","Africa","LatAmerica"])

    buildpage(m, policy)

    # next: add vector and panel example
    initpage(m, policy)

    return m
end

getpage = get_model

end


text:
using Test
using DataFrames
using Random
using Dates

rng = MersenneTwister(123)

times = 150
df = DataFrame("timestamp"=>1:times,
               "A"=>randn(rng, Float16, times),
               "B"=>sin.(rand(rng, times)),
               "C"=>1:times)
train_prop = 0.8
train_end = Int(times*train_prop)
train_df = df[1:train_end,:]
valid_df = df[train_end+1:end,:]

include("window_generator.jl")
@testset "historical window of length 1, future window 1" begin
    h = 1
    f = 1
    wg = WindowGenerator(h, f, train_df, valid_df, "C")
    @test wg.target_idx == [4]
    @test size(wg.train[1][1]) == (length(names(df)),h)
    @test size(wg.valid[1][1]) == (length(names(df)),h)
    @test length(wg.train) == train_end - (h+f-1) # truncates necessary first and last sequence points
    @test length(wg.valid) == times - train_end - (h+f-1)
    @test (wg.train[1][1][1]) == (wg.train[1][2][1] - h)
    @test (wg.train[end][1][1]) == (wg.train[end][2][1] - h)
    @test wg.train[end][1][1,end] == (wg.train[end][2][1,begin] - 1) # 1st pred is immediately after last hist

end

@testset "historical window of length 7, future window 5" begin
    h = 7
    f = 5
    wg = WindowGenerator(h, f, train_df, valid_df, "C")
    @test wg.target_idx == [4]
    @test size(wg.train[1][1]) == (length(names(df)),h)
    @test size(wg.valid[1][1]) == (length(names(df)),h)
    @test length(wg.train) == train_end - (h+f-1)
    @test length(wg.valid) == times - train_end - (h+f-1)
    @test (wg.train[1][1][1]) == (wg.train[1][2][1] - h)
    @test (wg.train[end][1][1]) == (wg.train[end][2][1] - h)
    @test wg.train[end][1][1,end] == (wg.train[end][2][1,begin] - 1) # 1st pred is immediately after last hist
end

@testset "multiple label columns" begin
    h = 16
    f = 3
    wg = WindowGenerator(h, f, train_df, valid_df, "C")
    @test length(wg.target_idx) == 1
    wg = WindowGenerator(h, f, train_df, valid_df, ["A","C"])
    @test length(wg.target_idx) == 2
    wg = WindowGenerator(h, f, train_df, valid_df, ["A","C","B"])
    @test length(wg.target_idx) == 3
    wg = WindowGenerator(h, f, train_df, valid_df, ["timestamp","A","C","B"])
    @test length(wg.target_idx) == 4
end

@testset "ignores nonexistent/duplicate columns" begin
    h = 5
    f = 3
    wg = WindowGenerator(h, f, train_df, valid_df, ["timestamp","A","C","B","D"])
    @test length(wg.target_idx) == 4

    wg = WindowGenerator(h, f, train_df, valid_df, ["timestamp","timestamp"])
    @test length(wg.target_idx) == 1
end

@testset "labels_indices are correct" begin
    valid_end = times - train_end # limited by size of validation set
    for h in 1:(valid_end - 1)
        for f in 1:(valid_end - h)
            wg = WindowGenerator(h, f, train_df, valid_df, ["timestamp","A","C","B"])
            @test f == length(wg.label_indices)
            @test first(wg.label_indices) == h + 1
            @test last(wg.label_indices) == h + f
        end
    end
end

# batch_ts
include("batch_ts.jl")
@testset "" begin

end

max_stars_repo_name: adinhobl/model-zoo

max_issues_repo_name: adinhobl/model-zoo

max_forks_repo_name: adinhobl/model-zoo


text:
import PDESolver.solvePDE

#=
# Startup function for 1 dof advection equation
"""
  This function invokes the solver for the advection equation, using the
  specified input file

  Inputs:
    input_file: a string containing the path to an input file, or just a file
                name.  If it is just a file name, it is taken to be in the
                users pwd.

  Outputs:
    mesh: an AbstractMesh object
    sbp: the SBP operator used in solving the equation
    eqn: the AbstractSolutionData object used to solve the equation
    opts: the options dictonary

"""
function run_advection(input_file::AbstractString)

  mesh, sbp, eqn, opts, pmesh = createObjects(input_file)
  solve_advection(mesh, sbp, eqn, opts, pmesh)

  return mesh, sbp, eqn, opts
end
=#

"""
  This function creates and initializes the mesh, sbp, eqn, and opts objects

  **Inputs**

   * opts: options dictionary

  **Outputs**

   * mesh: an AbstractMesh.  The concrete type is determined by the options
          dictionary
   * sbp: an AbstractOperator.  The concrete type is determined by the options
         dictionary
   * eqn: an AdvectionData object
   * opts: the options dictionary
   * pmesh: mesh used for preconditioning, can be same object as mesh
"""
function createObjects(opts::Dict)

  dim = opts["dimensions"]
  dofpernode = 1

  sbp, mesh, pmesh, Tsol, Tres, Tmsh, mesh_time = createMeshAndOperator(opts, dofpernode)

  # Create advection equation object
  Tdim = dim
  eqn = AdvectionData_{Tsol, Tres, Tdim, Tmsh}(mesh, sbp, opts)
  eqn.params.time.t_meshinit = mesh_time  # save the mesh init time

  # Initialize the advection equation
  init(mesh, sbp, eqn, opts)

  return mesh, sbp, eqn, opts, pmesh
end

"""
  Given fully initialized mesh, sbp, eqn, opts, this function solves
  the advection equations.  The 4 object should be obtained from
  createObjects().

  Specifically, it applies an initial condition and invokes a nonlinear
  solver according to the options dictionary.

  Inputs:
    mesh: an AbstractMesh
    sbp: an AbstractOperator
    eqn: an AdvectionData
    opts: the options dictionary.  This must be the options dictionary returned
          by createObjects().  Changing values in the options dictionary after
          calling createObjects() results in undefined behavior.
    pmesh: mesh used for preconditioning, can be same object as mesh.
           default value of mesh

"""
function solvePDE(mesh::AbstractMesh, sbp::AbstractOperator, eqn::AdvectionData,
                  opts::Dict, pmesh::AbstractMesh=mesh)

  myrank = mesh.myrank

  fill!(eqn.res, 0.0)
  fill!(eqn.res_vec, 0.0)

  q_vec = eqn.q_vec

  # Populate with initial conditions
  @mpi_master println("\nEvaluating initial condition")
  ICfunc_name = opts["IC_name"]
  ICfunc = ICDict[ICfunc_name]
  @mpi_master println("ICfunc = ", ICfunc)
  ICfunc(mesh, sbp, eqn, opts, q_vec)

  #writeSolutionFiles(mesh, sbp, eqn, opts, "solution_ic")

  if opts["calc_error"]
    @mpi_master println("\ncalculating error of file ", opts["calc_error_infname"],
            " compared to initial condition")
    # read in this processors portion of the solution
    vals = copy(eqn.q_vec)
    readSolutionFiles(mesh, sbp, eqn, opts, opts["calc_error_infname"], vals)

    err_vec = vals - eqn.q_vec
    err_local = calcNorm(eqn, err_vec)
    if myrank == 0
      outname = opts["calc_error_outfname"]
      @mpi_master println("printed err = ", err, " to file ", outname)
      f = open(outname, "w")
      println(f, err)
      close(f)
    end
  end

  if opts["calc_trunc_error"]  # calculate truncation error
    @mpi_master println("\nCalculating residual for truncation error")
    tmp = physicsRhs(mesh, sbp, eqn, opts, eqn.res_vec, (evalResidual,))
    if myrank == 0
      f = open("error_trunc.dat", "w")
      println(f, tmp)
      close(f)
    end
  end

  if opts["calc_havg"]
    h_avg = calcMeshH(mesh, sbp, eqn, opts)
    if myrank == 0  # TODO: make this globally accurate
      rmfile("havg.dat")
      f = open("havg.dat", "w")
      println(f, h_avg)
      close(f)
    end
  end

  if opts["perturb_ic"]
    @mpi_master println("\nPerturbing initial condition")
    perturb_mag = opts["perturb_mag"]
    for i=1:mesh.numDof
      q_vec[i] += perturb_mag*rand()
    end
  end

  res_vec_exact = deepcopy(q_vec)

  saveSolutionToMesh(mesh, q_vec)
  writeVisFiles(mesh, "solution_ic")

  if opts["calc_dt"]
    alpha_net = sqrt(eqn.params.alpha_x^2 + eqn.params.alpha_y^2)
    opts["delta_t"] = opts["CFL"]*mesh.min_el_size/alpha_net
  end

  # this must be last, because all previous calculations need to use the
  # original initial condition
  if opts["is_restart"]
    loadRestartState(mesh, sbp, eqn, opts, pmesh)
  end

  MPI.Barrier( mesh.comm)

  call_nlsolver(mesh, sbp, eqn, opts, pmesh)
  postproc(mesh, sbp, eqn, opts)

  MPI.Barrier(mesh.comm)

  return mesh, sbp, eqn, opts
end  # end function

"""
  This function does post processing, if requested by the input options.
  Typical post processing includes calculation of errors, norms of important
  quantities, writing of files. etc.

  Inputs:
    mesh
    sbp
    eqn
    opts
"""
function postproc(mesh, sbp, eqn, opts)

  ##### Do postprocessing ######
  myrank = mesh.myrank
  @mpi_master println("\nDoing postprocessing")

  Tsol = eltype(eqn.q)
  Tres = eltype(eqn.res)
  Tmsh = eltype(mesh.dxidx)

  if opts["do_postproc"] && opts["solve"]
    @mpi_master println("final time = ", eqn.t)
    exfname = opts["exact_soln_func"]
    if haskey(ICDict, exfname)
      exfunc = ICDict[exfname]
      q_exact = zeros(Tsol, mesh.numDof)
      exfunc(mesh, sbp, eqn, opts, q_exact)
      q_diff = eqn.q_vec - q_exact
      diff_norm = calcNorm(eqn, q_diff)
      discrete_norm = norm(q_diff/length(q_diff))
      discrete_norm = MPI.Allreduce(discrete_norm*discrete_norm, MPI.SUM, mesh.comm)
      discrete_norm = sqrt(discrete_norm)

      if myrank == 0
        println("solution error norm = ", diff_norm)
        println("solution discrete L2 norm = ", discrete_norm)
      end

      saveSolutionToMesh(mesh, real(q_diff))
      writeVisFiles(mesh, "solution_error")

      sol_norm = calcNorm(eqn, eqn.q_vec)
      exact_norm = calcNorm(eqn, q_exact)
      @mpi_master println("numerical solution norm = ", sol_norm)
      @mpi_master println("exact solution norm = ", exact_norm)

      # calculate the average mesh size
      h_avg = calcMeshH(mesh, sbp, eqn, opts)
      if myrank == 0
#        println("mesh.min_node_distance = ", mesh.min_node_dist)
        # print to file
        outname = opts["calc_error_outfname"]
        @mpi_master println("printed err = ", diff_norm, " to file ", outname)
        f = open(outname, "w")
        println(f, diff_norm, " ", h_avg)
        close(f)
      end
#=
      outname = opts["calc_error_outfname"]
      f = open(outname, "w")
      println(f, mesh.numEl, " ", diff_norm, " ", discrete_norm)
      close(f)
=#
    end  # end if haskey(ICname)
  end  # end if do_postproc

  return nothing
end