text:
using Chain
using DataFrames
using DataFramesMeta
using Dates
using StatsBase

"""
Takes a DataFrame of Twitter data and a date manipulation function to group on and calculates volume

# Arguments
- `twit_data::DataFrame`: Twitter data
- `date_func::Function`: Date manipulation function

# Returns
-  `DataFrame`: Containing date function output and n volume
"""
function get_volume_by(twit_data::DataFrame, date_func::Function)
    volume_by = @chain twit_data begin
        @transform(:temp = date_func.(:date))
        combine(groupby(_, :temp), nrow)
    end

    rename!(volume_by, :nrow => "n")
    rename!(volume_by, :temp => Symbol(date_func))

    return volume_by
end

"""
Takes a DataFrame of Twitter data and a date manipulation function to group on and calculates cumulative volume

# Arguments
- `twit_data::DataFrame`: Twitter data
- `date_func::Function`: Date manipulation function

# Returns
-  `DataFrame`: Containing date function output and n cumulative volume
"""
function get_cumulative_volume_by(twit_data::DataFrame, date_func::Function)
    volume_by = get_volume_by(twit_data, date_func)
    volume_by = @transform(volume_by, :n = cumsum(:n))

    return volume_by
end

"""
Takes a DataFrame of Twitter data and a date manipulation function to group on and calculates
number of quotes and retweets

# Arguments
- `twit_data::DataFrame`: Twitter data
- `date_func::Function`: Date manipulation function
- `manip_func::Function`: Column manipulation function (default mean)

# Returns
-  `DataFrame`: Containing date function output, nquotes and nretweets
"""
function get_interactions_by(twit_data::DataFrame, date_func::Function, manip_func::Function = mean)
    interactions = @chain twit_data begin
        @transform(:temp = date_func.(:date))
        @by(:temp, :nquotes = manip_func(:quote_count), :nretweets = manip_func(:retweet_count), :nreplies = manip_func(:reply_count))
    end

    rename!(interactions, :temp => Symbol(date_func))

    return interactions
end

"""
Takes a DataFrame of Twitter data and extracts mention data

# Arguments
- `twit_data::DataFrame`: Twitter data

# Returns
-  `DataFrame`: Containing original twitter data plus mentions and Nmentions columns
"""
function add_mentions(twit_data::DataFrame)
    twit_data.mentions = get_mentions.(lowercase.(twit_data.text))
    twit_data.n_mentions = length.(twit_data.mentions)

    return twit_data
end

"""
Takes a DataFrame of Twitter data and a date manipulation function to group on and calculates
number of mentions

# Arguments
- `twit_data::DataFrame`: Twitter data
- `date_func::Function`: Date manipulation function
- `manip_func::Function`: Column manipulation function (default mean)

# Returns
-  `DataFrame`: Containing date function output and number of mentions
"""
function get_mentions_by(twit_data::DataFrame, date_func::Function, manip_func::Function = mean)
    twit_data = add_mentions(twit_data)

    mentions_by = @chain twit_data begin
        @transform(:temp = date_func.(:date))
        @by(:temp, :nmentions = manip_func(:n_mentions))
    end

    rename!(mentions_by, :temp => Symbol(date_func))

    return mentions_by
end

"""
Takes a DataFrame of Twitter data and plots the 'ratio' or number retweets against the number of quotes

# Arguments
- `twit_data::DataFrame`: Twitter data

# Output
A plot of retweets against quotes
"""
function generate_ratio_plot(twit_data::DataFrame)
    ratio_df = data(twit_data) *
        mapping(:quote_count => "N Quotes",:retweet_count => "N Retweets")
    return ratio_df
end

"""
Takes a DataFrame of Twitter data and date manipulation function and produces a plot
of the volume of tweets grouped by the date manipulation function

# Arguments
- `twit_data::DataFrame`: Twitter data
- `date_func::Function`: Date manipulation function

# Output
A plot of tweet volume by date function grouping
"""
function generate_volume_by_plot(twit_data::DataFrame, date_func::Function)
    volume = data(get_volume_by(twit_data, date_func)) *
        visual(Lines) *
        mapping(Symbol(date_func) => string(date_func), :n => "Number of posts")

    return volume
end

"""
Takes a DataFrame of Twitter data and date manipulation function and produces a plot
of the cumulative volume of tweets grouped by the date manipulation function

# Arguments
- `twit_data::DataFrame`: Twitter data
- `date_func::Function`: Date manipulation function

# Output
A plot of cumulative tweet volume by date function grouping
"""
function generate_cumulative_volume_by_plot(twit_data::DataFrame, date_func::Function)
    cumulative_volume = data(get_cumulative_volume_by(twit_data, date_func)) *
        visual(Lines) *
        mapping(Symbol(date_func) => string(date_func), :n => "Cumulative number of posts")

    return cumulative_volume
end

"""
Takes a DataFrame of Twitter data and date manipulation function and produces a plot
of the average number of mentions per tweet grouped by the date manipulation function

# Arguments
- `twit_data::DataFrame`: Twitter data
- `date_func::Function`: Date manipulation function

# Output
A plot of average mentions by date function grouping
"""
function generate_avg_mentions_by_plot(twit_data::DataFrame, date_func::Function)
    mentions = data(get_mentions_by(twit_data, date_func, mean)) *
        visual(Lines) *
        mapping(Symbol(date_func) => string(date_func), :nmentions => "Average number of mentions per tweet")

    return mentions
end

"""
Takes a DataFrame of Twitter data and date manipulation function and produces a plot
of the total number of mentions grouped by the date manipulation function

# Arguments
- `twit_data::DataFrame`: Twitter data
- `date_func::Function`: Date manipulation function

# Output
A plot of total mentions by date function grouping
"""
function generate_total_mentions_by_plot(twit_data::DataFrame, date_func::Function)
    mentions = data(get_mentions_by(twit_data, date_func, sum)) *
        visual(Lines) *
        mapping(Symbol(date_func) => string(date_func), :nmentions => "Total number of mentions")

    return mentions
end

"""
Takes a DataFrame of Twitter data and date manipulation function and produces a plot
of the average number of retweets per tweet grouped by the date manipulation function

# Arguments
- `twit_data::DataFrame`: Twitter data
- `date_func::Function`: Date manipulation function

# Output
A plot of average retweets by date function grouping
"""
function generate_avg_retweets_by_plot(twit_data::DataFrame, date_func::Function)
    df = select(get_interactions_by(twit_data, date_func, mean), Symbol(date_func), :nretweets)
    mentions = data(df) *
        visual(Lines) *
        mapping(Symbol(date_func) => string(date_func), :nretweets => "Average number of retweets per tweet")

    return mentions
end

"""
Takes a DataFrame of Twitter data and date manipulation function and produces a plot
of the total number of grouped by the date manipulation function

# Arguments
- `twit_data::DataFrame`: Twitter data
- `date_func::Function`: Date manipulation function

# Output
A plot of total retweets by date function grouping
"""
function generate_total_retweets_by_plot(twit_data::DataFrame, date_func::Function)
    df = select(get_interactions_by(twit_data, date_func, sum), Symbol(date_func), :nretweets)
    mentions = data(df) *
        visual(Lines) *
        mapping(Symbol(date_func) => string(date_func), :nretweets => "Total number of retweets")

    return mentions
end

"""
Takes a DataFrame of Twitter data and date manipulation function and produces a plot
of the average number of quotes per tweet grouped by the date manipulation function

# Arguments
- `twit_data::DataFrame`: Twitter data
- `date_func::Function`: Date manipulation function

# Output
A plot of average quote tweets by date function grouping
"""
function generate_avg_quotes_by_plot(twit_data::DataFrame, date_func::Function)
    df = select(get_interactions_by(twit_data, date_func, mean), Symbol(date_func), :nquotes)
    mentions = data(df) *
        visual(Lines) *
        mapping(Symbol(date_func) => string(date_func), :nquotes => "Average number of quotes per tweet")

    return mentions
end

"""
Takes a DataFrame of Twitter data and date manipulation function and produces a plot
of the total number of quotes grouped by the date manipulation function

# Arguments
- `twit_data::DataFrame`: Twitter data
- `date_func::Function`: Date manipulation function

# Output
A plot of total quote tweets by date function grouping
"""
function generate_total_quotes_by_plot(twit_data::DataFrame, date_func::Function)
    df = select(get_interactions_by(twit_data, date_func, sum), Symbol(date_func), :nquotes)
    mentions = data(df) *
        visual(Lines) *
        mapping(Symbol(date_func) => string(date_func), :nquotes => "Total number of quotes")

    return mentions
end

function get_unique_users(twit_data::DataFrame, date_func::Function)
    users = @chain twit_data begin
        @transform(:temp = date_func.(:date))
        @by(:temp, :nusers = length(unique(:author_id)))
    end

    rename!(users, :temp => Symbol(date_func))

    return users
end

function generate_unique_users(twit_data::DataFrame, date_func::Function)
    df = get_unique_users(twit_data, date_func)
    users = data(df) *
        visual(Lines) *
        mapping(Symbol(date_func) => string(date_func), :nusers => "N unique users")

    return users
end

function get_quote_volume(twit_data::DataFrame, date_func::Function)
    quotes = @chain twit_data begin
        @transform(:temp = date_func.(:date))
        @by(:temp, :nquotes = count(i -> (i=="quoted" || "quoted" in split(strip(i,['[', ']']),", ")), :sourcetweet_type))
    end

    rename!(quotes, :temp => Symbol(date_func))

    return quotes
end

function generate_quote_volume(twit_data::DataFrame, date_func::Function)
    df = select(get_interactions_by(twit_data, date_func, sum), Symbol(date_func), :nquotes)
    quotes = data(df) *
        visual(Lines) *
        mapping(Symbol(date_func) => string(date_func), :nquotes => "N quote tweets")

    return quotes
end

function get_retweet_volume(twit_data::DataFrame, date_func::Function)
    df = @chain twit_data begin
        @transform(:temp = date_func.(:date))
        @by(:temp, :nretweets = count(i -> (i=="retweeted" || "retweeted" in split(strip(i,['[', ']']),", ")), :sourcetweet_type))
    end

    rename!(df, :temp => Symbol(date_func))

    return df
end

function generate_retweet_volume(twit_data::DataFrame, date_func::Function)
    df = select(get_interactions_by(twit_data, date_func, sum), Symbol(date_func), :nretweets)
    quotes = data(df) *
        visual(Lines) *
        mapping(Symbol(date_func) => string(date_func), :nretweets => "N quote tweets")

    return quotes
end

function get_reply_volume(twit_data::DataFrame, date_func::Function)
    df = @chain twit_data begin
        @transform(:temp = date_func.(:date))
        @by(:temp, :nreplies = count(i -> (i!="NA"), :in_reply_to_user_id))
    end

    rename!(df, :temp => Symbol(date_func))

    return df
end

function generate_reply_volume(twit_data::DataFrame, date_func::Function)
    df = select(get_interactions_by(twit_data, date_func, sum), Symbol(date_func), :nreplies)
    quotes = data(df) *
        visual(Lines) *
        mapping(Symbol(date_func) => string(date_func), :nreplies => "N reply tweets")

    return quotes
end

function get_unique_user_rate(twit_data::DataFrame, date_func::Function)
    df = @chain twit_data begin
        @transform(:temp1 = minute.(:date), :temp2 = hour.(:date), :temp3 = dayofyear.(:date))
        @by([:temp1, :temp2, :temp3], :nunique = length(unique(:author_id)))
        @by(:temp3, :nuniqueAvg = mean(:nunique))
    end

    rename!(df, :temp3 => Symbol(date_func))

    return df
end

function generate_unique_user_rate(twit_data::DataFrame, date_func::Function)
    df = get_unique_user_rate(twit_data, date_func)
    quotes = data(df) *
        visual(Lines) *
        mapping(Symbol(date_func) => string(date_func), :nuniqueAvg => "Unique user tweet rate")

    return quotes
end

function get_quote_rate(twit_data::DataFrame, date_func::Function)
    is_quote(type) = type == "quoted" || "quoted" in split(strip(type,['[', ']']),", ")
    twit_data = filter(:sourcetweet_type => is_quote, twit_data)
    df = @chain twit_data begin
        @transform(:temp1 = minute.(:date), :temp2 = hour.(:date), :temp3 = dayofyear.(:date))
        @by([:temp1, :temp2, :temp3], :nunique = length(unique(:author_id)))
        @by(:temp3, :nuniqueAvg = mean(:nunique))
    end

    rename!(df, :temp3 => Symbol(date_func))

    return df
end

function generate_quote_rate(twit_data::DataFrame, date_func::Function)
    df = get_quote_rate(twit_data, date_func)
    quotes = data(df) *
        visual(Lines) *
        mapping(Symbol(date_func) => string(date_func), :nuniqueAvg => "Quote tweet rate")

    return quotes
end

function get_retweet_rate(twit_data::DataFrame, date_func::Function)
    is_retweet(type) = type == "retweeted" || "retweeted" in split(strip(type,['[', ']']),", ")
    twit_data = filter(:sourcetweet_type => is_retweet, twit_data)
    df = @chain twit_data begin
        @transform(:temp1 = minute.(:date), :temp2 = hour.(:date), :temp3 = dayofyear.(:date))
        @by([:temp1, :temp2, :temp3], :nunique = length(unique(:author_id)))
        @by(:temp3, :nuniqueAvg = mean(:nunique))
    end

    rename!(df, :temp3 => Symbol(date_func))

    return df
end

function generate_retweet_rate(twit_data::DataFrame, date_func::Function)
    df = get_retweet_rate(twit_data, date_func)
    retweets = data(df) *
        visual(Lines) *
        mapping(Symbol(date_func) => string(date_func), :nuniqueAvg => "Retweet rate")

    return retweets
end

function get_reply_rate(twit_data::DataFrame, date_func::Function)
    is_reply(type) = type != "NA"
    twit_data = filter(:in_reply_to_user_id => is_reply, twit_data)
    df = @chain twit_data begin
        @transform(:temp1 = minute.(:date), :temp2 = hour.(:date), :temp3 = dayofyear.(:date))
        @by([:temp1, :temp2, :temp3], :nunique = length(unique(:author_id)))
        @by(:temp3, :nuniqueAvg = mean(:nunique))
    end

    rename!(df, :temp3 => Symbol(date_func))

    return df
end

function generate_reply_rate(twit_data::DataFrame, date_func::Function)
    df = get_reply_rate(twit_data, date_func)
    replys = data(df) *
        visual(Lines) *
        mapping(Symbol(date_func) => string(date_func), :nuniqueAvg => "reply tweet rate")

    return replys
end

function get_quote_rate_2(twit_data::DataFrame, date_func::Function)
    df = select(get_interactions_by(twit_data, date_func, sum), Symbol(date_func), :nquotes)
    df = @chain df begin
        @transform(:nquotes = :nquotes/1440)
    end
    return df
end

function generate_quote_rate_2(twit_data::DataFrame, date_func::Function)
    df = get_quote_rate_2(twit_data, date_func)
    replys = data(df) *
        visual(Lines) *
        mapping(Symbol(date_func) => string(date_func), :nquotes => "quote tweet rate")

    return replys
end

function get_retweet_rate_2(twit_data::DataFrame, date_func::Function)
    df = select(get_interactions_by(twit_data, date_func, sum), Symbol(date_func), :nretweets)
    df = @chain df begin
        @transform(:nretweets = :nretweets/1440)
    end
    return df
end

function generate_retweet_rate_2(twit_data::DataFrame, date_func::Function)
    df = get_retweet_rate_2(twit_data, date_func)
    replys = data(df) *
        visual(Lines) *
        mapping(Symbol(date_func) => string(date_func), :nretweets => "retweet rate")

    return replys
end

function get_reply_rate_2(twit_data::DataFrame, date_func::Function)
    df = select(get_interactions_by(twit_data, date_func, sum), Symbol(date_func), :nreplies)
    df = @chain df begin
        @transform(:nreplies = :nreplies/1440)
    end
    print(df)
    return df
end

function generate_reply_rate_2(twit_data::DataFrame, date_func::Function)
    df = get_reply_rate_2(twit_data, date_func)
    replys = data(df) *
        visual(Lines) *
        mapping(Symbol(date_func) => string(date_func), :nreplies => "reply rate")

    return replys
end

function get_cum_quote_volume(twit_data::DataFrame, date_func::Function)
    df = select(get_interactions_by(twit_data, date_func, sum), Symbol(date_func), :nquotes)
    df = @transform(df, :n = cumsum(:nquotes))
    return df
end

function generate_cum_quote_volume(twit_data::DataFrame, date_func::Function)
    cumulative_volume = data(get_cum_quote_volume(twit_data, date_func)) *
        visual(Lines) *
        mapping(Symbol(date_func) => string(date_func), :n => "Cumulative number of quote tweets")

    return cumulative_volume
end

function get_cum_retweet_volume(twit_data::DataFrame, date_func::Function)
    df = select(get_interactions_by(twit_data, date_func, sum), Symbol(date_func), :nretweets)
    df = @transform(df, :n = cumsum(:nretweets))
    return df
end

function generate_cum_retweet_volume(twit_data::DataFrame, date_func::Function)
    cumulative_volume = data(get_cum_retweet_volume(twit_data, date_func)) *
        visual(Lines) *
        mapping(Symbol(date_func) => string(date_func), :n => "Cumulative number of retweets")

    return cumulative_volume
end

function get_cum_reply_volume(twit_data::DataFrame, date_func::Function)
    df = select(get_interactions_by(twit_data, date_func, sum), Symbol(date_func), :nreplies)
    df = @transform(df, :n = cumsum(:nreplies))
    return df
end

function generate_cum_reply_volume(twit_data::DataFrame, date_func::Function)
    cumulative_volume = data(get_cum_reply_volume(twit_data, date_func)) *
        visual(Lines) *
        mapping(Symbol(date_func) => string(date_func), :n => "Cumulative number of replies")

    return cumulative_volume
end

function wrap_in_makie(out_of_AoG, annotations, use_mad = false, fig = Figure())
    ax1 = Axis(fig[1,1],
        xlabel = annotations[1],
        ylabel = annotations[2],
        title = annotations[3]
        )

    draw!(ax1, out_of_AoG)

    if use_mad
        mean_data = mean(out_of_AoG.data[2])
        mad_data = mad(out_of_AoG.data[2])

        hlines!(ax1, mean_data)
        hspan!(ax1, mean_data-mad_data , mean_data+mad_data, color = (:blue, 0.2))
    end

    return fig
end

function cumulative_sorting(df, column)
    @chain df begin
        @subset(:true_type .!= "retweeted")
        groupby(:author_id)
        combine(column => sum => :tot)
        @orderby(-:tot)
        @transform(_,
        :cumulative = cumsum(:tot),
        :rownumber = rownumber.(eachrow(_)))
    end
end

function graph2(df, date_func)
    fig = Figure()

    g1 = wrap_in_makie(generate_quote_volume(df, date_func), ["Day of year", "N quotes", "Daily quote volume"], false, fig[1,1])
    g2 = wrap_in_makie(generate_retweet_volume(df, date_func), ["Day of year", "N retweets", "Daily retweet volume"], false, fig[1,2])
    g3 = wrap_in_makie(generate_reply_volume(df, date_func), ["Day of year", "N replies", "Daily reply volume"], false, fig[2,1])

    return fig
end

function graph4(df, date_func)
    fig = Figure()
    g1 = wrap_in_makie(generate_quote_rate_2(df, date_func), ["Day of year", "Avg n quotes per min", "Daily quote rate"], true, fig[1,1])
    g2 = wrap_in_makie(generate_retweet_rate_2(df, date_func), ["Day of year", "Avg n retweets per min", "Daily retweet rate"], true, fig[1,2])
    g3 = wrap_in_makie(generate_reply_rate_2(df, date_func), ["Day of year", "Avg n replies per min", "Daily reply rate"], true, fig[2,1])

    return fig
end

function graph5(df, date_func)
    fig = Figure()
    g1 = wrap_in_makie(generate_cum_quote_volume(df, date_func), ["Day of year", "Cumulative n quotes", "Daily quote volume"], false, fig[1,1])
    g2 = wrap_in_makie(generate_cum_retweet_volume(df, date_func), ["Day of year", "Cumulative n retweets", "Daily retweet volume"], false, fig[1,2])
    g3 = wrap_in_makie(generate_cum_reply_volume(df, date_func), ["Day of year", "Cumulative n replies", "Daily reply volume"], false, fig[2,1])

    return fig
end


text:
function write_problem_to_file(problem_name, g, ζ_vec, I_sols, f_sols)

outfile = joinpath(dirname(@__FILE__),
                    "problems",  problem_name * ".jl")
f = open(outfile, "w")

println(f, "# Automatically generated test file")
println(f, "function ", problem_name, "()\n")

print(f, "g = ")
if typeof(g) <: Vector
    println(f, "[\n", join([@sprintf "%0.15f" x for x in g], ",\n"), "]\n")
else
    println(f, g)
end

println(f, "ζ_vec = ", ζ_vec, "\n")

println(f, "I_sols = [\n", join([string(x) for x in I_sols], ",\n"), "]\n")

println(f, "f_sols = [\n", join([@sprintf "%0.15f" x for x in f_sols], ",\n"), "]\n")

println(f, "return g, ζ_vec, I_sols, f_sols\n")
println(f, "end")

close(f)

end

ext: jl


text:
# export readmseed, readmseed!

function seed_cleanup!(S::SeisData, BUF::SeisIOBuf)
  trunc_x!(S)
  fill!(getfield(BUF, :hdr_old), zero(UInt8))
  setfield!(BUF, :r1_old, zero(Int16))
  setfield!(BUF, :r2_old, zero(Int16))
  return nothing
end

function parsemseed!(S::SeisData, sid::IO, v::Int64, nx_new::Int64, nx_add::Int64)
  while !eof(sid)
    parserec!(S, BUF, sid, v, nx_new, nx_add)
  end
  seed_cleanup!(S, BUF)
  return S
end

function read_mseed_file!(S::SeisData, fname::String, v::Int64, nx_new::Int64, nx_add::Int64)
  io = open(fname, "r")
  skip(io, 6)
  c = read(io, UInt8)
  if c in (0x44, 0x52, 0x4d, 0x51)
    seekstart(io)
    parsemseed!(S, io, v, nx_new, nx_add)
    close(io)
  else
    error("Invalid file type!")
  end
  return nothing
end


text:
############################################################################################
# Define abstract types
"""
$(TYPEDEF)
Abstract super type for Tuning structs of differentiable functions.
"""
abstract type AbstractDifferentiableTune end

"""
    $(FUNCTIONNAME)(objective, θᵤ)
Compute log density of 'objective' at 'θᵤ'.

# Examples
```julia
```

"""
function log_density end

"""
    $(FUNCTIONNAME)(objective, θᵤ)
Compute log density and gradient of 'objective' at 'θᵤ'.

# Examples
```julia
```

"""
function log_density_and_gradient end

"""
    $(FUNCTIONNAME)(objective, θᵤ)
Compute log density, gradient and hessian of 'objective' at 'θᵤ'.

# Examples
```julia
```

"""
function log_density_and_hessian end

############################################################################################
# Include
include("tune_analytic.jl")
include("tune_automatic.jl")
include("diffobjective.jl")
include("results.jl")
include("checks.jl")

############################################################################################
# Export
export Objective, update!, log_density, log_density_and_gradient#,
#log_density_and_hessian


text:
module FranklinCLI

using Franklin, Comonicon
Comonicon.set_brief_length!(120)

"""
Generate a new folder for the website. It contains the
skeleton of a website that can be processed by Franklin.

# Arguments

- `path`: the path where you want to create the website.

# Options

- `-t,--template <name>`: template name, default is "basic"

# Flags

- `-v,--verbose`: add this flag to print verbose info.
"""
@cast function new(path; template="basic", verbose::Bool=false)
    Franklin.newsite(path; template=template, verbose=verbose, changedir=false)
    return
end

"""
Runs Franklin in the current directory.

# Options

- `--host <ip>`: the host to use for the local server, default is "127.0.0.1".
- `-p, --port <port>`: the port to use for the local server (should pick a number between 8000 and 9000).

# Flags

- `--clear`: whether to remove any existing output directory.
- `--verbose`: whether to display messages.
- `--single`: whether to run a single pass or run continuously.
- `--prerender`: whether to pre-render javascript (KaTeX and highlight.js).
- `--nomess`: suppresses all messages (internal use).
- `--eval_all`: whether to force re-evaluation of all code blocks.
- `--silent`: switch this on to suppress all output (including eval statements).
"""
@cast function serve(;
    host::String="127.0.0.1",
    port::Int=8000,
    clear::Bool=false, verbose::Bool=false,
    single::Bool=false,
    prerender::Bool=false,
    eval_all::Bool=false,
    silent::Bool=false,
    )

    Franklin.serve(;clear=clear, verb=verbose,
        port=port, single=single, prerender=prerender,
        eval_all=eval_all,
        silent=silent, host=host,
    )
end

@main name="franklin" doc="""
static site generator. Simple, fast, maths with KaTeX, code evaluation,
LaTeX-like commands, optional pre-rendering, in Julia.
"""

end


text:
import ..Rhmc:RHMC,get_order,get_α,get_β,get_α0,get_β_inverse,get_α_inverse, get_α0_inverse

struct StaggeredFermiAction{Dim,Dirac,fermion,gauge,Nf} <: FermiAction{Dim,Dirac,fermion,gauge}
    hascovnet::Bool
    covneuralnet::Union{Nothing,CovNeuralnet{Dim}}
    diracoperator::Dirac
    Nf::Int64
    rhmc_info_for_action::Union{Nothing,RHMC}
    rhmc_info_for_MD::Union{Nothing,RHMC}
    _temporary_fermionfields::Vector{fermion}
    _temporary_gaugefields::Vector{gauge}

    function StaggeredFermiAction(D::Dirac_operator{Dim},hascovnet,covneuralnet,parameters_action) where Dim
        @assert haskey(parameters_action,"Nf") "parameters for action should have the keyword Nf"
        Nf_in = parameters_action["Nf"]
        if Dim == 2 #Number of zero modes is different.
            Nf = 2*Nf_in
        elseif Dim == 4
            Nf = Nf_in
        end

        if Nf == 4 || Nf == 8 # 8 flavors if phi (Mdag M)^{-1} phi
            rhmc_info_for_action = nothing
            rhmc_info_for_MD = nothing
        else
            #for action: r_action
            #Nf = 8 -> alpha = 1 -> power x^{1/2}
            #Nf = 2 -> alpha = 1/4 -> power x^1/8
            #Nf = 1 -> alpha = 1/8  -> power x^1/16
            order = Nf //16

            rhmc_info_for_action = RHMC(order,n=15)

            #for MD: r_MD
            #Nf = 8 -> alpha = 1 -> power x^{1}
            #Nf = 2 -> alpha = 1/4 -> power x^1/4
            #Nf = 1 -> alpha = 1/8  -> power x^1/8
            order = Nf // 8
            #rhmcorder = 8 ÷ Nf
            rhmc_info_for_MD = RHMC(order,n=10)

            N_action = get_order(rhmc_info_for_action)
            N_MD = get_order(rhmc_info_for_MD)
        end

        if Nf == 4
            num = 1
        elseif Nf == 8
            num = 1
        else
            num = maximum((N_action,N_MD))+1
        end
        num += 4

        x = D._temporary_fermi[1]
        xtype = typeof(x)
        _temporary_fermionfields = Array{xtype,1}(undef,num)
        for i=1:num
            _temporary_fermionfields[i] = similar(x)
        end

        Utemp = D.U[1]
        Utype = typeof(Utemp)
        numU = 2
        _temporary_gaugefields = Array{Utype,1}(undef,numU)
        for i=1:numU
            _temporary_gaugefields[i] = similar(Utemp)
        end

        return new{Dim,typeof(D),xtype,Utype,Nf}(hascovnet,covneuralnet,D,Nf,rhmc_info_for_action,rhmc_info_for_MD,
                        _temporary_fermionfields,
                        _temporary_gaugefields)
    end
end

function evaluate_FermiAction(fermi_action::StaggeredFermiAction{Dim,Dirac,fermion,gauge,Nf},U,ϕ::AbstractFermionfields) where {Dim,Dirac,fermion,gauge,Nf}
    W = fermi_action.diracoperator(U)
    WdagW = DdagD_Staggered_operator(W)

    rhmc = fermi_action.rhmc_info_for_action

    N = get_order(rhmc)

    x = fermi_action._temporary_fermionfields[end-N]
    vec_x = fermi_action._temporary_fermionfields[end-N+1:end]
    for j=1:N
        clear_fermion!(vec_x[j])
    end

    vec_β = get_β_inverse(rhmc)
    vec_α = get_α_inverse(rhmc)
    α0 = get_α0_inverse(rhmc)
    shiftedcg(vec_x,vec_β,x,WdagW,ϕ,eps = W.eps_CG,maxsteps= W.MaxCGstep, verbose = W.verbose_print)
    clear_fermion!(x)
    add_fermion!(x,α0,ϕ)
    for j=1:N
        αk = vec_α[j]
        add_fermion!(x,αk,vec_x[j])
    end

    Sf = dot(x,x)
    return real(Sf)
end

function evaluate_FermiAction(fermi_action::StaggeredFermiAction{Dim,Dirac,fermion,gauge,8},U,ϕ::AbstractFermionfields) where {Dim,Dirac,fermion,gauge}
    W = fermi_action.diracoperator(U)
    η = fermi_action._temporary_fermionfields[1]
    solve_DinvX!(η,W',ϕ)
    Sf = dot(η,η)
    return real(Sf)
end

function evaluate_FermiAction(fermi_action::StaggeredFermiAction{Dim,Dirac,fermion,gauge,4},U,ϕ::AbstractFermionfields) where {Dim,Dirac,fermion,gauge}
    W = fermi_action.diracoperator(U)
    η = fermi_action._temporary_fermionfields[1]
    solve_DinvX!(η,W',ϕ)
    Sf = dot(η,η)
    return real(Sf)
end

function gauss_sampling_in_action!(η::AbstractFermionfields,U,fermi_action::StaggeredFermiAction{Dim,Dirac,fermion,gauge,Nf}) where {Dim,Dirac,fermion,gauge,Nf}
    gauss_distribution_fermion!(η)
end

function gauss_sampling_in_action!(η::AbstractFermionfields,U,fermi_action::StaggeredFermiAction{Dim,Dirac,fermion,gauge,4}) where {Dim,Dirac,fermion,gauge}
    evensite = false
    W = fermi_action.diracoperator(U)
    temp = fermi_action._temporary_fermionfields[1]
    gauss_distribution_fermion!(η)
    mul!(temp,W',η)
    clear_fermion!(temp,evensite)
    solve_DinvX!(η,W',temp)
end

function sample_pseudofermions!(ϕ::AbstractFermionfields,U,fermi_action::StaggeredFermiAction{Dim,Dirac,fermion,gauge,8},ξ::AbstractFermionfields) where {Dim,Dirac,fermion,gauge}
    W = fermi_action.diracoperator(U)
    mul!(ϕ,W',ξ)
    set_wing_fermion!(ϕ)
end

function sample_pseudofermions!(ϕ::AbstractFermionfields,U,fermi_action::StaggeredFermiAction{Dim,Dirac,fermion,gauge,4},ξ::AbstractFermionfields) where {Dim,Dirac,fermion,gauge}
    W = fermi_action.diracoperator(U)
    mul!(ϕ,W',ξ)
    set_wing_fermion!(ϕ)
end

function sample_pseudofermions!(ϕ::AbstractFermionfields,U,fermi_action::StaggeredFermiAction{Dim,Dirac,fermion,gauge,Nf},ξ::AbstractFermionfields) where {Dim,Dirac,Nf,fermion,gauge}
    W = fermi_action.diracoperator(U)
    WdagW = DdagD_Staggered_operator(W)

    rhmc = fermi_action.rhmc_info_for_action

    N = get_order(rhmc)

    x = ϕ #fermi_action._temporary_fermionfields[1]
    vec_x = fermi_action._temporary_fermionfields[end-N+1:end]
    for j=1:N
        clear_fermion!(vec_x[j])
    end

    vec_β = get_β(rhmc)
    vec_α = get_α(rhmc)
    α0 = get_α0(rhmc)
    #set_wing_fermion!(ξ)

    shiftedcg(vec_x,vec_β,x,WdagW,ξ,eps = W.eps_CG,maxsteps= W.MaxCGstep, verbose = W.verbose_print)
    clear_fermion!(ϕ)
    add_fermion!(ϕ,α0,ξ)
    for j=1:N
        αk = vec_α[j]
        add_fermion!(ϕ,αk,vec_x[j])
    end
    set_wing_fermion!(ϕ)

    #error("sample_pseudofermions!(ϕ,fermi_action) is not implemented in type ϕ:$(typeof(ϕ)), fermi_action:$(typeof(fermi_action))")
end

function calc_UdSfdU!(UdSfdU::Vector{<: AbstractGaugefields},fermi_action::StaggeredFermiAction{Dim,Dirac,fermion,gauge,Nf},U::Vector{<: AbstractGaugefields},ϕ::AbstractFermionfields) where {Dim,Dirac,fermion,Nf,gauge}
    @assert Nf != 4 && Nf != 8
    clear_U!(UdSfdU)

    W = fermi_action.diracoperator(U)
    WdagW = DdagD_Staggered_operator(W)

    rhmc = fermi_action.rhmc_info_for_MD
    N = get_order(rhmc)
    x = fermi_action._temporary_fermionfields[end-N]
    vec_x = fermi_action._temporary_fermionfields[end-N+1:end]
    Y = fermi_action._temporary_fermionfields[2]
    for j=1:N
        clear_fermion!(vec_x[j])
    end
    vec_β = get_β_inverse(rhmc)
    vec_α = get_α_inverse(rhmc)
    #α0 = get_α0_inverse(rhmc)

    shiftedcg(vec_x,vec_β,x,WdagW,ϕ,eps = W.eps_CG,maxsteps= W.MaxCGstep, verbose = W.verbose_print)

    for j=1:N
        set_wing_fermion!(vec_x[j])

        calc_UdSfdU_fromX!(UdSfdU,Y,fermi_action,U,vec_x[j],coeff = vec_α[j])
    end

    set_wing_U!(UdSfdU)

    #error("calc_UdSfdU!(UdSfdU,fermi_action,U) is not implemented in type dSfdU:$(typeof(UdSfdU)), fermi_action:$(typeof(fermi_action)), ϕ:$(typeof(ϕ))")
end

function calc_UdSfdU!(UdSfdU::Vector{<: AbstractGaugefields},fermi_action::StaggeredFermiAction{Dim,Dirac,fermion,gauge,4},U::Vector{<: AbstractGaugefields},ϕ::AbstractFermionfields) where {Dim,Dirac,fermion,gauge}
#function cald_UdSfdU!(UdSfdU::Vector{<: AbstractGaugefields},fermi_action::StaggeredFermiAction{Dim,Dirac,fermion,gauge,4},U::Vector{<: AbstractGaugefields},ϕ::AbstractFermionfields) where {Dim,Dirac,fermion,gauge}
    calc_UdSfdU_Nf_4_8!(UdSfdU,fermi_action,U,ϕ)
end

function calc_UdSfdU!(UdSfdU::Vector{<: AbstractGaugefields},fermi_action::StaggeredFermiAction{Dim,Dirac,fermion,gauge,8},U::Vector{<: AbstractGaugefields},ϕ::AbstractFermionfields) where {Dim,Dirac,fermion,gauge}
    calc_UdSfdU_Nf_4_8!(UdSfdU,fermi_action,U,ϕ)
end

function calc_UdSfdU_Nf_4_8!(UdSfdU::Vector{<: AbstractGaugefields},fermi_action::StaggeredFermiAction{Dim,Dirac,fermion,gauge,Nf},U::Vector{<: AbstractGaugefields},ϕ::AbstractFermionfields) where {Dim,Dirac,fermion,gauge,Nf}
    W = fermi_action.diracoperator(U)
    WdagW = DdagD_Staggered_operator(W)
    X = fermi_action._temporary_fermionfields[1]
    Y = fermi_action._temporary_fermionfields[2]
    #X = (D^dag D)^(-1) ϕ
    #
    solve_DinvX!(X,WdagW,ϕ)
    clear_U!(UdSfdU)

    calc_UdSfdU_fromX!(UdSfdU,Y,fermi_action,U,X)
    set_wing_U!(UdSfdU)
end

function calc_UdSfdU_fromX!(UdSfdU::Vector{<: AbstractGaugefields},Y,fermi_action::StaggeredFermiAction{Dim,Dirac,fermion,gauge,Nf},U,X;coeff = 1) where {Dim,Dirac,fermion,gauge,Nf}
    W = fermi_action.diracoperator(U)
    mul!(Y,W,X)

    temp0_f = fermi_action._temporary_fermionfields[3]
    temp0_g = fermi_action._temporary_gaugefields[1]

    #println(coeff)
    for μ=1:Dim
        #!  Construct U(x,mu)*P1

        # U_{k,μ} X_{k+μ}
        Xplus = shift_fermion(X,μ)

        Us = staggered_U(U[μ],μ)
        mul!(temp0_f,Us,Xplus)
        #println(temp0_f[1,1,1,1,1,1])

        #U_{k,μ} X_{k+μ}) ⊗ Y_k
        mul!(temp0_g,temp0_f,Y')
        add_U!(UdSfdU[μ],coeff*0.5,temp0_g)
        #println(temp2_g[1,1,1,1,1,1])
        #mul!(dSfdU[μ],U[μ]',temp2_g) #additional term

        #!  Construct P2*U_adj(x,mu)
        # Y_{k+μ}^dag U_{k,μ}^dag
        Yplus = shift_fermion(Y,μ)
        mul!(temp0_f,Yplus',Us')

        #X_k ⊗ Y_{k+μ}^dag U_{k,μ}^dag
        mul!(temp0_g,X,temp0_f)
        add_U!(UdSfdU[μ],coeff*0.5,temp0_g)

        #mul!(temp3_g,U[μ]',temp2_g)
        #add_U!(dSfdU[μ],temp3_g)
    end

end


text:
function residuals(A, b, shifts, x)
  nshifts = size(shifts, 1);
  r = [ (b - A * x[:,i] - shifts[i] * x[:,i]) for i = 1 : nshifts ];
  return r;
end

function test_cg_lanczos()
  cg_tol = 1.0e-6;

  # Cubic splines matrix.
  n = 10;
  A = spdiagm(-1 => ones(n-1), 0 => 4*ones(n), 1 => ones(n-1))
  b = A * [1:n;];
  b_norm = norm(b);

  (x, stats) = cg_lanczos(A, b, itmax=n);
  resid = norm(b - A * x) / b_norm;
  @printf("CG_Lanczos: Relative residual: %8.1e\n", resid);
  @test(resid <= cg_tol);
  @test(stats.solved);

  # Test negative curvature detection.
  A[n-1,n-1] = -4.0
  (x, stats) = cg_lanczos(A, b, check_curvature=true)
  @test(stats.status == "negative curvature")
  A[n-1,n-1] = 4.0

  shifts=[1:6;];

  (x, stats) = cg_lanczos_shift_seq(A, b, shifts, itmax=n);
  r = residuals(A, b, shifts, x);
  resids = map(norm, r) / b_norm;
  @printf("CG_Lanczos: Relative residuals with shifts:");
  for resid in resids
    @printf(" %8.1e", resid);
  end
  @printf("\n");
  @test(all(resids .<= cg_tol));
  @test(stats.solved);

  # Test negative curvature detection.
  shifts = [-4; -3; 2]
  (x, stats) = cg_lanczos_shift_seq(A, b, shifts, check_curvature=true, itmax=n);
  @test(stats.flagged == [true, true, false])

  # Code coverage.
  (x, stats) = cg_lanczos(Matrix(A), b);
  (x, stats) = cg_lanczos_shift_seq(Matrix(A), b, [1:6;], verbose=true);
  show(stats);

  # Test b == 0
  (x, stats) = cg_lanczos(A, zeros(size(A,1)))
  @test x == zeros(size(A,1))
  @test stats.status == "x = 0 is a zero-residual solution"
  (x, stats) = cg_lanczos_shift_seq(A, zeros(size(A,1)), [1:6;])
  @test x == zeros(size(A,1), 6)
  @test stats.status == "x = 0 is a zero-residual solution"

  # Test integer values
  A = ceil.(Int, A) + I
  b = A * ones(Int, size(A,1))
  (x, stats) = cg_lanczos(A, b)
  @test stats.solved
end

test_cg_lanczos()


text:
module Escher

if VERSION < v"0.4.0-dev"
    using Docile
end

if VERSION < v"0.4.0-dev"
    using Markdown
else
    using Base.Markdown
end

export @md_str, @md_mstr
#@docstrings

using Compat
using Patchwork
using Reactive
using JSON

# Polymer Setup
custom_elements() =
    readstring(Pkg.dir("Escher", "assets", "basics/basics.html"))

include("basics/macros.jl")
include("basics/tile.jl")
include("basics/util.jl")
include("basics/length.jl")
include("basics/signal.jl")

include("basics/layout.jl")
include("basics/content.jl")
include("basics/typography.jl")
include("basics/embellishment.jl")
include("basics/behavior.jl")
include("basics/window.jl")
include("basics/component.jl")

include("library/markdown.jl")
include("library/tex.jl")
include("library/widgets.jl")
include("library/layout2.jl")
include("library/slideshow.jl")
include("library/codemirror.jl")

include("deprecate.jl")

# We need to set up convert methods for other packages (e.g. Gadfly, SymPy, Images)
# This will work on 0.3 automatically, when using 0.4 with precompilation,
# an explicit call to conditional_setup is needed after loading the other packages.

if VERSION >= v"0.4.0-"
    external_setup() = include(joinpath(dirname(@__FILE__), "basics", "lazyload.jl"))
else
    external_setup() = nothing
end

serve(port=5555) =
    include(joinpath(dirname(@__FILE__), "cli", "serve.jl"))(port)

# 3rd party package interop
include("basics/lazyload.jl")

end


text:
# ---
# title: 399. Evaluate Division
# id: problem399
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Union Find, Graph
# link: <https://leetcode.com/problems/evaluate-division/description/>
# hidden: true
# ---
#
# You are given an array of variable pairs `equations` and an array of real
# numbers `values`, where `equations[i] = [Ai, Bi]` and `values[i]` represent
# the equation `Ai / Bi = values[i]`. Each `Ai` or `Bi` is a string that
# represents a single variable.
#
# You are also given some `queries`, where `queries[j] = [Cj, Dj]` represents
# the `jth` query where you must find the answer for `Cj / Dj = ?`.
#
# Return _the answers to all queries_. If a single answer cannot be determined,
# return `-1.0`.
#
# **Note:** The input is always valid. You may assume that evaluating the
# queries will not result in division by zero and that there is no
# contradiction.
#
#
#
# **Example 1:**
#
#
#
#     Input: equations = [["a","b"],["b","c"]], values = [2.0,3.0], queries = [["a","c"],["b","a"],["a","e"],["a","a"],["x","x"]]
#     Output: [6.00000,0.50000,-1.00000,1.00000,-1.00000]
#     Explanation:
#     Given: _a / b = 2.0_ , _b / c = 3.0_
#     queries are: _a / c = ?_ , _b / a = ?_ , _a / e = ?_ , _a / a = ?_ , _x / x = ?_
#     return: [6.0, 0.5, -1.0, 1.0, -1.0 ]
#
#
# **Example 2:**
#
#
#
#     Input: equations = [["a","b"],["b","c"],["bc","cd"]], values = [1.5,2.5,5.0], queries = [["a","c"],["c","b"],["bc","cd"],["cd","bc"]]
#     Output: [3.75000,0.40000,5.00000,0.20000]
#
#
# **Example 3:**
#
#
#
#     Input: equations = [["a","b"]], values = [0.5], queries = [["a","b"],["b","a"],["a","c"],["x","y"]]
#     Output: [0.50000,2.00000,-1.00000,-1.00000]
#
#
#
#
# **Constraints:**
#
#   * `1 <= equations.length <= 20`
#   * `equations[i].length == 2`
#   * `1 <= Ai.length, Bi.length <= 5`
#   * `values.length == equations.length`
#   * `0.0 < values[i] <= 20.0`
#   * `1 <= queries.length <= 20`
#   * `queries[i].length == 2`
#   * `1 <= Cj.length, Dj.length <= 5`
#   * `Ai, Bi, Cj, Dj` consist of lower case English letters and digits.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
module FESpacesInterfacesTests

using FillArrays
using Test
using Gridap.Arrays
using Gridap.TensorValues
using Gridap.ReferenceFEs
using Gridap.Geometry
using Gridap.FESpaces
using Gridap.CellData
using LinearAlgebra

order = 1
domain =(0,1,0,1)
partition = (3,3)
model = CartesianDiscreteModel(domain,partition)
grid_topology = get_grid_topology(model)
polytopes = get_polytopes(grid_topology)
reffes = [LagrangianRefFE(Float64,p,order) for p in polytopes]
dirichlet_tags = ["tag_1","tag_6"]
V = GradConformingFESpace(reffes,model,dirichlet_tags)
test_fe_space(V)

zh = zero(V)
@test is_a_fe_function(zh)
@test has_constraints(V) == false

cellids = [1,3,5,2]
@test get_cell_dofs(V,cellids) == [[-1, 1, 4, 5], [2, 3, 6, 7], [5, 6, 9, 10], [1, 2, 5, 6]]

cellidsL = cellids
cellidsR = [2,4,3,1]
cellidsS = SkeletonPair(cellidsL,cellidsR)
@test isa(get_cell_dofs(V,cellidsS)[1],BlockArrayCoo)
@test get_cell_dofs(V,cellidsS)[1] == [-1, 1, 4, 5, 1, 2, 5, 6]

@test get_cell_axes(V) == get_cell_axes_with_constraints(V)

cell_constr = get_cell_constraints(V)
@test cell_constr == [Matrix(I,4,4) for cell in 1:num_cells(model)]

cell_constr = get_cell_constraints(V,cellids)
@test cell_constr == [Matrix(I,4,4) for cell in 1:length(cellids)]

cell_constr = get_cell_constraints(V,cellidsS)
@test isa(cell_constr[1],BlockArrayCoo)

@test get_cell_isconstrained(V,cellids) == Fill(false,length(cellids))

@test get_cell_isconstrained(V,cellidsS) == Fill(false,length(cellids))

cellmat = [rand(4,4) for cell in 1:num_cells(model)]
cellvec = [rand(4) for cell in 1:num_cells(model)]
cellmatvec = pair_arrays(cellmat,cellvec)
cellids = identity_vector(num_cells(model))
matdata = (cellmat,cellids,cellids)
vecdata = (cellvec,cellids)
matvecdata = (cellmatvec,cellids,cellids)

@test cellmat === attach_constraints_rows(V,cellmat,cellids)
@test cellmat === attach_constraints_cols(V,cellmat,cellids)
test_fe_space(V,matvecdata,matdata,vecdata)

end # module


text:
mutable struct Points
	position::Array{Float64,2}
	velocity::Array{Float64,2}
	acceleration::Array{Float64,2}
	mass::Array{Float64,1}
	size::Int
end

function acceleration_update(points::Points, force_field::Array{Float64,3})
	for i = 1:points.size
		points.acceleration[:,i] = force_field[:,ceil(Int,points.position[1,i]/dx),ceil(Int,points.position[2,i]/dx)] / points.mass[i]
	end
end

function velocity_update(points::Points)
	for i = 1:points.size
		points.velocity[:,i] += points.acceleration[:,i] * dt #(accelration+old_acc)/2
	end
end

function position_update(points::Points)
	for i = 1:points.size
		points.position[:,i] += points.velocity[:,i] * dt  #+ acceleration*dt*dt/mass
	end
end

function assign_mass(density::Array{Float64,2}, points::Points, dx)
	for i = 1:points.size
		density[ceil(Int, points.position[1,i]/dx),ceil(Int, points.position[2,i]/dx)] = points.mass[i]
	end
end

function gauss_seidel(potential::Array{Float64,2}, density::Array{Float64,2})
	for i = 2:grid-1
		for j = 2:grid-1
			potential[i,j] = -1/4 * (dx^2*4*pi*GravitationalConstant*density[i,j] + (- potential[i-1,j] - potential[i+1,j] - potential[i,j-1] - potential[i,j+1]))
		end
	end
end

function residuum(potentail::Array{Float64,2}, density::Array{Float64,2})
	r = 0.
	for i = 2:grid-1
		for j = 2:grid-1
			r += (4*pi*GravitationalConstant*density[i,j] - (potentail[i+1,j] + potentail[i-1,j] + potentail[i,j+1] + potentail[i,j-1] - 4potentail[i,j]) * 1/dx^2)^2
		end
	end
	return sqrt(r/(grid-2)^2)
end

using WriteVTK

N = 20
grid = 200
ϵ = 0.0001 #residual
ω = 1.8 #SOR
γ = 0.9 #upwind
itermax = 2000
dims = 2
GravitationalConstant = 6.673e-8
cells = WriteVTK.MeshCell[MeshCell(VTKCellTypes.VTK_VERTEX, [i]) for i = 1:N]

points = Points(rand(dims,N), zeros(dims,N), zeros(dims,N), rand(N)*100, N)
density = zeros(Float64, grid, grid)
potential = zeros(Float64, grid, grid)
force_field = zeros(Float64, 2, grid, grid)

dx = 1/grid
dt = 1

for i = 1:60
	outfile = vtk_grid(string("my_vtr_file",i), [i*dx for i=0:grid],[i*dx for i=0:grid]) do vtk
		vtk_cell_data(vtk, density, "Density")
		vtk_cell_data(vtk, potential, "Potential")
		vtk_cell_data(vtk, vcat(force_field, zeros(1,grid,grid)), "Force")
	end

	outfile = vtk_grid(string("my_vtu_file",i), points.position, cells) do vtk
		vtk_cell_data(vtk, points.mass, "Mass")
		vtk_cell_data(vtk, vcat(points.velocity, zeros(1,N)), "Velocity")

		vtk_point_data(vtk, points.mass, "Pointmass")
		vtk_point_data(vtk, vcat(points.velocity, zeros(1,N)), "Pointvelocity")
	end

	density = zeros(Float64, grid, grid)
	assign_mass(density, points, dx)
	for n=1:itermax
		gauss_seidel(potential, density)
		r = residuum(potential, density)
		println(r)
		if r < ϵ
			break
		end
	end

	for i=2:grid-1
		for j=2:grid-1
			force_field[1,i,j] = -(potential[i+1,j] - potential[i-1,j])/dx
			force_field[2,i,j] = -(potential[i,j+1] - potential[i,j-1])/dx
		end
	end

	acceleration_update(points,force_field)
	velocity_update(points)
	position_update(points)
end


text:
struct CompilerMetadata{StaticNDRange, CheckBounds, I, NDRange, Iterspace}
    groupindex::I
    ndrange::NDRange
    iterspace::Iterspace

    # CPU variant
    function CompilerMetadata{NDRange, CB}(idx, ndrange, iterspace) where {NDRange, CB}
        if ndrange !== nothing
            ndrange = CartesianIndices(ndrange)
        end
        return new{NDRange, CB, typeof(idx), typeof(ndrange), typeof(iterspace)}(idx, ndrange, iterspace)
    end

    # GPU variante: index is given implicit
    function CompilerMetadata{NDRange, CB}(ndrange, iterspace) where {NDRange, CB}
        if ndrange !== nothing
            ndrange = CartesianIndices(ndrange)
        end
        return new{NDRange, CB, Nothing, typeof(ndrange), typeof(iterspace)}(nothing, ndrange, iterspace)
    end
end

@inline __iterspace(cm::CompilerMetadata)  = cm.iterspace
@inline __groupindex(cm::CompilerMetadata) = cm.groupindex
@inline __groupsize(cm::CompilerMetadata) = size(workitems(__iterspace(cm)))
@inline __dynamic_checkbounds(::CompilerMetadata{NDRange, CB}) where {NDRange, CB} = CB <: DynamicCheck
@inline __ndrange(::CompilerMetadata{NDRange}) where {NDRange<:StaticSize}  = CartesianIndices(get(NDRange))
@inline __ndrange(cm::CompilerMetadata{NDRange}) where {NDRange<:DynamicSize} = cm.ndrange
@inline groupsize(ctx) = __groupsize(ctx)
@inline __workitems_iterspace(ctx) = workitems(__iterspace(ctx))

include("compiler/contract.jl")


text:
using AWS
using AWSTools.CloudFormation: stack_output
using AWSBatch
using Dates
using HTTP: HTTP
using Memento
using Memento.TestUtils: @test_log, @test_nolog
using Mocking
using Test

using AWS.AWSExceptions: AWSException

Mocking.activate()

# Controls the running of various tests: "local", "batch"
const TESTS = strip.(split(get(ENV, "TESTS", "local"), r"\s*,\s*"))

# Run the tests on a stack created with the "test/batch.yml" CloudFormation template
# found in AWSClusterMangers.jl
const AWS_STACKNAME = get(ENV, "AWS_STACKNAME", "")
const STACK = !isempty(AWS_STACKNAME) ? stack_output(AWS_STACKNAME) : Dict()
const JOB_TIMEOUT = 900
const LOG_TIMEOUT = 30

const JULIA_BAKED_IMAGE = let
    output = read(`git ls-remote --tags https://github.com/JuliaLang/julia`, String)
    tags = split(replace(output, r".*\/" => ""))
    versions = VersionNumber.(filter(v -> !endswith(v, "^{}"), tags))
    latest_version = maximum(versions)
    docker_tag = VERSION > latest_version ? "nightly" : "$VERSION"

    "468665244580.dkr.ecr.us-east-1.amazonaws.com/julia-baked:$docker_tag"
end

Memento.config!("debug"; fmt="[{level} | {name}]: {msg}")
const logger = getlogger(AWSBatch)
setlevel!(logger, "info")

# We've been having issues with the log stream being created but no log events are present.
# - https://gitlab.invenia.ca/invenia/AWSBatch.jl/issues/28
# - https://gitlab.invenia.ca/invenia/AWSBatch.jl/issues/30
#
# This function allows for us to wait if logs are not present but avoids blocking when
# logs are ready.
#
# Note: The timeout duration expects the job has reached the SUCCEEDED or FAILED state and
# is not expected to last long enough for a job to complete running.
function wait_for_log_events(job::BatchJob)
    events = nothing

    # Convert to Float64 until this is merged and we no longer use versions of Julia
    # without PR (probably Julia 1.5+): https://github.com/JuliaLang/julia/pull/35103
    timedwait(Float64(LOG_TIMEOUT); pollint=Float64(5)) do
        events = log_events(job)

        # Note: These warnings should assist in determining the special circumstances
        # the log events not being present. Eventually warnings should be removed.
        if events === nothing
            notice(logger, "Log stream for $(job.id) does not exist")
        elseif isempty(events)
            notice(logger, "No log events for $(job.id)")
        end

        # Wait for log stream to exist and contain at least one event
        events !== nothing && !isempty(events)
    end

    return events
end

include("mock.jl")

@testset "AWSBatch.jl" begin
    if "local" in TESTS
        # need to define these to make sure we don't inadvertently try to talk to AWS
        withenv("AWS_ACCESS_KEY_ID" => "", "AWS_SECRET_ACCESS_KEY" => "") do
            include("compute_environment.jl")
            include("job_queue.jl")
            include("log_event.jl")
            include("job_state.jl")
            include("batch_job.jl")
            include("run_batch.jl")
        end
    else
        warn(logger, "Skipping \"local\" tests. Set `ENV[\"TESTS\"] = \"local\"` to run.")
    end

    if "batch" in TESTS && !isempty(AWS_STACKNAME)
        @testset "AWS Batch" begin
            info(logger, "Running AWS Batch tests")

            @testset "Job Submission" begin
                definition = "aws-batch-test"

                # Append the job ID to the definition when running on the CI. Doing this
                # will allow this test to successfully reuse the job definition later when
                # concurrent CI jobs are running AWS Batch tests at the same time.
                if haskey(ENV, "CI_JOB_ID")
                    definition *= "-" * ENV["CI_JOB_ID"]
                end

                job = run_batch(;
                    name = "aws-batch-test",
                    definition = definition,
                    queue = STACK["JobQueueArn"],
                    image = JULIA_BAKED_IMAGE,
                    vcpus = 1,
                    memory = 1024,
                    role = STACK["JobRoleArn"],
                    cmd = `julia -e 'println("Hello World!")'`,
                    parameters = Dict{String, String}("region" => "us-east-1"),
                )

                @test wait(job, [AWSBatch.SUCCEEDED]; timeout=JOB_TIMEOUT) == true
                @test status(job) == AWSBatch.SUCCEEDED

                events = wait_for_log_events(job)
                @test length(events) == 1
                @test first(events).message == "Hello World!"

                # Test job details were set correctly
                job_details = describe(job)
                @test job_details["jobName"] == "aws-batch-test"
                @test occursin(STACK["JobQueueArn"], job_details["jobQueue"])
                @test job_details["parameters"] == Dict("region" => "us-east-1")

                # Test job definition and container parameters were set correctly
                job_definition = JobDefinition(job)
                @test isregistered(job_definition) == true

                job_definition_details = first(describe(job_definition)["jobDefinitions"])

                @test job_definition_details["jobDefinitionName"] == definition
                @test job_definition_details["status"] == "ACTIVE"
                @test job_definition_details["type"] == "container"

                container_properties = job_definition_details["containerProperties"]
                @test container_properties["image"] == JULIA_BAKED_IMAGE
                @test container_properties["vcpus"] == 1
                @test container_properties["memory"] == 1024
                @test container_properties["command"] == [
                    "julia",
                    "-e",
                    "println(\"Hello World!\")"
                ]
                @test container_properties["jobRoleArn"] == STACK["JobRoleArn"]

                # Reuse job definition
                job = run_batch(;
                    name = "aws-batch-test",
                    definition = definition,
                    queue = STACK["JobQueueArn"],
                    image = JULIA_BAKED_IMAGE,
                    vcpus = 1,
                    memory = 1024,
                    role = STACK["JobRoleArn"],
                    cmd = `julia -e 'println("Hello World!")'`,
                    parameters = Dict{String, String}("region" => "us-east-1"),
                )

                @test wait(job, [AWSBatch.SUCCEEDED]; timeout=JOB_TIMEOUT) == true
                @test status(job) == AWSBatch.SUCCEEDED

                # Test job definition and container parameters were set correctly
                job_definition_2 = JobDefinition(job)
                @test job_definition_2 == job_definition

                deregister(job_definition)
            end

            @testset "Job registration disallowed" begin
                @test_throws BatchEnvironmentError run_batch(;
                    name = "aws-batch-no-job-registration-test",
                    queue = STACK["JobQueueArn"],
                    image = JULIA_BAKED_IMAGE,
                    role = STACK["JobRoleArn"],
                    cmd = `julia -e 'println("Hello World!")'`,
                    parameters = Dict{String, String}("region" => "us-east-1"),
                    allow_job_registration = false,
                )
            end

            @testset "Job parameters" begin
                # Use parameter substitution placeholders in the command field
                command = Cmd(["julia", "-e", "Ref::juliacmd"])

                # Set a default output string when registering the job definition
                job_definition = register(
                    "aws-batch-parameters-test";
                    image=JULIA_BAKED_IMAGE,
                    role=STACK["JobRoleArn"],
                    vcpus=1,
                    memory=1024,
                    cmd=command,
                    parameters=Dict("juliacmd" => "println(\"Default String\")"),
                )

                # Override the default output string
                job = run_batch(;
                    name = "aws-batch-parameters-test",
                    definition = job_definition,
                    queue = STACK["JobQueueArn"],
                    image = JULIA_BAKED_IMAGE,
                    vcpus = 1,
                    memory = 1024,
                    role = STACK["JobRoleArn"],
                    cmd = command,
                    parameters=Dict("juliacmd" => "println(\"Hello World!\")"),
                )

                @test wait(state -> state < AWSBatch.SUCCEEDED, job; timeout=JOB_TIMEOUT)
                @test status(job) == AWSBatch.SUCCEEDED

                events = wait_for_log_events(job)
                @test length(events) == 1
                @test first(events).message == "Hello World!"

                # Test job details were set correctly
                job_details = describe(job)
                @test job_details["parameters"] == Dict(
                    "juliacmd" => "println(\"Hello World!\")"
                )

                job_definition = JobDefinition(job)
                job_definition_details = first(describe(job_definition)["jobDefinitions"])
                job_definition_details["parameters"] = Dict(
                    "juliacmd" => "println(\"Default String\")"
                )

                container_properties = job_definition_details["containerProperties"]
                @test container_properties["command"] == ["julia", "-e", "Ref::juliacmd"]

                # Deregister job definition
                deregister(job_definition)
            end

            @testset "Array job" begin
                job = run_batch(;
                    name = "aws-batch-array-job-test",
                    definition = "aws-batch-array-job-test",
                    queue = STACK["JobQueueArn"],
                    image = JULIA_BAKED_IMAGE,
                    vcpus = 1,
                    memory = 1024,
                    role = STACK["JobRoleArn"],
                    cmd = `julia -e 'println("Hello World!")'`,
                    num_jobs = 3,
                )

                @test wait(state -> state < AWSBatch.SUCCEEDED, job; timeout=JOB_TIMEOUT)
                @test status(job) == AWSBatch.SUCCEEDED

                # Test array job was submitted properly
                status_summary = Dict(
                    "STARTING" => 0, "FAILED" => 0, "RUNNING" => 0, "SUCCEEDED" => 3,
                    "RUNNABLE" => 0, "SUBMITTED" => 0, "PENDING" => 0,
                )
                job_details = describe(job)
                @test job_details["arrayProperties"]["statusSummary"] == status_summary
                @test job_details["arrayProperties"]["size"] == 3

                # No log stream will exist for the parent job
                events = log_events(job)
                @test events === nothing

                # Test logs for each individual job that is part of the job array
                for i in 0:2
                    job_id = "$(job.id):$i"
                    events = wait_for_log_events(BatchJob(job_id))

                    @test length(events) == 1
                    @test first(events).message == "Hello World!"
                end

                # Deregister the job definition
                job_definition = JobDefinition(job)
                deregister(job_definition)
            end

            @testset "Failed Job" begin
                info(logger, "Testing job failure")

                job = run_batch(;
                    name = "aws-batch-failed-job-test",
                    definition = "aws-batch-failed-job-test",
                    queue = STACK["JobQueueArn"],
                    image = JULIA_BAKED_IMAGE,
                    vcpus = 1,
                    memory = 1024,
                    role = STACK["JobRoleArn"],
                    cmd = `julia -e 'error("Testing job failure")'`,
                )

                job_definition = JobDefinition(job)
                @test isregistered(job_definition) == true

                @test_throws BatchJobError wait(
                    job,
                    [AWSBatch.SUCCEEDED];
                    timeout=JOB_TIMEOUT
                )

                events = wait_for_log_events(job)
                @test first(events).message == "ERROR: Testing job failure"

                deregister(job_definition)
            end
        end
    else
        warn(
            logger,
            "Skipping \"batch\" tests. Set `ENV[\"TESTS\"] = \"batch\"` and " *
            "`ENV[\"AWS_STACKNAME\"]` to run."
        )
    end
end


text:
#=
    "high level" i.e. using Platform
=#
using HetaSimulator, Plots

################################## Model Upload ###########################################
# heta_update_dev()

platform = load_platform("$HetaSimulatorDir/cases/story_2", rm_out=false);
# platform.models
# platform.scenarios

################################## Single Simulation ######################################

model = platform.models[:nameless]

Scenario(model; tspan = (0., 200.)) |> sim |> plot

################################## load and simulate scenarios  #####################

### load scenarios from csv
scn_csv = read_scenarios("./cases/story_2/scenarios.csv")
# scn_xlsx = read_scenarios("./cases/story_2/scenarios.xlsx")
add_scenarios!(platform, scn_csv)

### sim
sim1 = sim(platform, scenarios = [:three]);
sim1 |> plot
sim_all = sim(platform);
sim_all |> plot

### Measurements
# load from csv to model
measurements = read_measurements("./cases/story_2/measurements.csv");
add_measurements!(platform, measurements)

### Fitting

# loss(sim1, measurements.dataone) # why?

fit1 = fit(platform, [:k1=>0.1,:k2=>0.2,:k3=>0.3], scenarios = [:dataone])
fit2 = fit(platform, [:k1=>0.1,:k2=>0.2,:k3=>0.3], scenarios = [:withdata2])
fit3 = fit(platform, [:k1=>0.1,:k2=>0.2,:k3=>0.3], scenarios = [:dataone,:withdata2])
fit_all = fit(platform, [:k1=>0.1,:k2=>0.2,:k3=>0.3])

# save_as_heta("optimal.heta", fit_all, append=false)
# save_as_heta("optimal.heta", [:k1=>0.1,:k2=>0.2,:k3=>0.3])
# save_as_heta("optimal.heta", (k1=0.1,k2=0.2,k3=0.3))
optim(fit_all)

### Monte-Carlo Simulations

mcsim1 = mc(platform, [:k2=>Normal(1e-3,1e-4), :k3=>Normal(1e-4,1e-5)], 1000)
plot(mcsim1)

mcsim2 = mc(platform, [:k1=>0.01, :k2=>Normal(1e-3,1e-4), :k3=>Uniform(1e-4,1e-2)], 1000)
plot(mcsim2)

################################## Monte-Carlo Parallel #####################

using Distributed
addprocs(2)
@everywhere push!(LOAD_PATH, "Y:/")
@everywhere using Pkg
@everywhere Pkg.activate("Y:/HetaSimulator.jl")
@everywhere using HetaSimulator

mcsim0 = mc(platform, [:k2=>Normal(1e-3,1e-4), :k3=>Normal(1e-4,1e-5)], 20)
mcsim1 = mc(platform, [:k2=>Normal(1e-3,1e-4), :k3=>Normal(1e-4,1e-5)], 150, parallel_type=EnsembleDistributed())


text:
# The Computer Language Benchmarks Game
# https://salsa.debian.org/benchmarksgame-team/benchmarksgame/
#
# contributed by David Campbell
# based on the Go version
# modified by Jarret Revels, Alex Arslan, Yichao Yu

using Printf

import Base.isless

function count_data(data::AbstractString, n::Int)
    counts = Dict{SubString{String}, Int}()
    top = length(data) - n + 1
    @inbounds for i = 1:top
        s = SubString(data, i, i+n-1)
        token = Base.ht_keyindex2!(counts, s)
        if token > 0
            counts.vals[token] += 1
        else
            Base._setindex!(counts, 1, s, -token)
        end
    end
    return counts
end

function count_one(data::AbstractString, s::AbstractString)
    d = count_data(data, length(s))
    return haskey(d, s) ? d[s] : 0
end

struct KNuc
    name::SubString{String}
    count::Int
end

# sort down
function isless(x::KNuc, y::KNuc)
    if x.count == y.count
        return x.name > y.name
    end
    x.count > y.count
end

function sorted_array(m::Dict{<:AbstractString, Int})
    kn = Vector{KNuc}(undef, length(m))
    i = 1
    for elem in m
        kn[i] = KNuc(elem...)
        i += 1
    end
    sort(kn)
end

function print_knucs(a::Array{KNuc, 1})
    sum = 0
    for kn in a
        sum += kn.count
    end
    for kn in a
        @printf("%s %.3f\n", kn.name, 100.0kn.count/sum)
    end
    println()
end

function perf_k_nucleotide(io = stdin)
    three = ">THREE "
    while true
        line = readline(io)
        if length(line) >= length(three) && line[1:length(three)] == three
            break
        end
    end
    data = collect(read(io, String))
    # delete the newlines and convert to upper case
    i, j = 1, 1
    while i <= length(data)
        if data[i] != '\n'
            data[j] = uppercase(data[i])
            j += 1
        end
        i += 1
    end
    str = join(data[1:j-1], "")

    arr1 = sorted_array(count_data(str, 1))
    arr2 = sorted_array(count_data(str, 2))

    print_knucs(arr1)
    print_knucs(arr2)

    v = ["GGT", "GGTA", "GGTATT", "GGTATTTTAATT", "GGTATTTTAATTTATAGT"]
    counts = zeros(length(v))
    # Could threads this loop but seems slower?
    for i in 1:length(v)
        counts[i] = count_one(str, v[i])
    end
    for (chain, c) in zip(v, counts)
        @printf("%d\t%s\n", c, chain)
    end
end

perf_k_nucleotide()
#perf_k_nucleotide(open("knucleotide-input.txt"))


text:
function random_walk(X::AbstractMatrix,α,step)
    X̄=X./sum(X,dims=2)
    Xs=Matrix{Float32}[X̄]
    for i=2:step
        push!(Xs,Xs[end]*X̄+α*X̄)
    end
    cat(Xs...,dims=4)
end

max_stars_repo_name: yeruoforever/DrugRepositioning

max_issues_repo_name: yeruoforever/DrugRepositioning

max_forks_repo_name: yeruoforever/DrugRepositioning


text:
# *****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License. You
# may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable
# law or agreed to in writing, software distributed under the License is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the License for the specific language
# governing permissions and limitations under the License.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

module ConfigAST

export defineASTParams, defineAST

using SISLES
using SISLES.GenerativeModel
using AdaptiveStressTesting

const COST_ACCEL_NEAR_RA = -5000.0

function defineASTParams(;
                         max_steps::Int64 = 51,
                         rsg_length::Int64 = 3,
                         init_seed::Int64 = 0,
                         reset_seed::Union{Void,Int64} = nothing)
  p = ASTParams()
  p.max_steps = max_steps
  p.rsg_length = rsg_length
  p.init_seed = init_seed
  p.reset_seed = reset_seed

  return p
end

function get_reward_custom(prob::Float64, event::Bool, terminal::Bool, dist::Float64,
                            ast::AdaptiveStressTest, sim::AbstractGenerativeModel)
  r = log(prob)

  #penalize for maneuvering too close to RA
  pr = sim.pr
  if any(pr[i].accel_near_RA for i=1:length(pr))
    r += COST_ACCEL_NEAR_RA
  end

  if event
    r += 0.0
  elseif terminal #incur distance cost only if !event && terminal
    r += -dist
  end
  return r
end

function defineAST(sim, p::ASTParams)
  return AdaptiveStressTest(p, sim, GenerativeModel.initialize, GenerativeModel.update,
                 GenerativeModel.isterminal, get_reward_custom)
end

end #module

max_stars_repo_name: sisl/RLESCAS

max_issues_repo_name: sisl/RLESCAS

max_forks_repo_name: sisl/RLESCAS


text:
using Test
using PlanktonIndividuals

@testset "PlanktonIndividuals" begin
    @testset "Unit tests" begin
        include("grid_test.jl")
        include("field_test.jl")
        include("parameter_test.jl")
        include("output_test.jl")
    end
    @testset "Example tests" begin
        include("example_3D_test.jl")
        include("example_2D_test.jl")
        include("example_1D_test.jl")
        include("example_0D_test.jl")
    end
end


text:
"""
    get_robust_cov_stats()

    Return all robust covariance estimators.
"""
get_all_robust_cov_stats() = Set([:white, :nw, :hc0, :hc1, :hc2, :hc3])
get_needed_robust_cov_stats(s::String) = return get_needed_robust_cov_stats([s])
get_needed_robust_cov_stats(s::Symbol) = return get_needed_robust_cov_stats([s])
get_needed_robust_cov_stats(s::Vector{String}) = return get_needed_robust_cov_stats(Symbol.(lowercase.(s)))
get_needed_robust_cov_stats(::Vector{Any}) = return get_needed_robust_cov_stats([:none])
get_needed_robust_cov_stats(::Set{Any}) = return get_needed_robust_cov_stats(Set([:none]))
get_needed_robust_cov_stats(s::Set{Symbol}) = return get_needed_robust_cov_stats(collect(s))

function get_needed_robust_cov_stats(s::Vector{Symbol})

    needed_white = Vector{Symbol}()
    needed_hac = Vector{Symbol}()

    length(s) == 0 && return (needed_white, needed_hac)
    :none in s && return (needed_white, needed_hac)
    if :all in s
        s = collect(get_all_robust_cov_stats())
    end

    :white in s && push!(needed_white, :white)
    :hc0 in s && push!(needed_white, :hc0)
    :hc1 in s && push!(needed_white, :hc1)
    :hc2 in s && push!(needed_white, :hc2)
    :hc3 in s && push!(needed_white, :hc3)

    :nw in s && push!(needed_hac, :nw)

    return return (needed_white, needed_hac)

end

"""
    get_all_model_stats()

    Returns all statistics availble for the fitted model.
"""
get_all_model_stats() = Set([:coefs, :sse, :mse, :sst, :rmse, :aic, :sigma, :t_statistic, :vif, :r2, :adjr2, :stderror, :t_values, :p_values, :ci])

get_needed_model_stats(req_stats::String) = return get_needed_model_stats([req_stats])
get_needed_model_stats(req_stats::Symbol) = return get_needed_model_stats(Set([req_stats]))
get_needed_model_stats(req_stats::Vector{String}) = return get_needed_model_stats(Symbol.(lowercase.(req_stats)))
get_needed_model_stats(::Vector{Any}) = return get_needed_model_stats([:none])
get_needed_model_stats(::Set{Any}) = return get_needed_model_stats(Set([:none]))
get_needed_model_stats(req_stats::Set{Symbol}) = get_needed_model_stats(collect(req_stats))

"""
    get_needed_model_stats(req_stats::Vector{Symbol})

    return the list of needed statistics given the list of statistics about the model the caller wants.
"""
function get_needed_model_stats(req_stats::Vector{Symbol})
    needed = Set([:coefs, :sse, :mse])
    full = get_all_model_stats()
    unique!(req_stats)

    length(req_stats) == 0 && return needed
    :none in req_stats && return needed
    :all in req_stats && return full

    :sst in req_stats && push!(needed, :sst)
    :rmse in req_stats && push!(needed, :rmse)
    :aic in req_stats && push!(needed, :aic)
    :sigma in req_stats && push!(needed, :sigma)
    :t_statistic in req_stats && push!(needed, :t_statistic)
    :vif in req_stats && push!(needed, :vif)

    if :r2 in req_stats
        push!(needed, :sst)
        push!(needed, :r2)
    end
    if :adjr2 in req_stats
        push!(needed, :sst)
        push!(needed, :r2)
        push!(needed, :adjr2)
    end
    if :stderror in req_stats
        push!(needed, :sigma)
        push!(needed, :stderror)
    end
    if :t_values in req_stats
        push!(needed, :sigma)
        push!(needed, :stderror)
        push!(needed, :t_values)
    end
    if :p_values in req_stats
        push!(needed, :sigma)
        push!(needed, :stderror)
        push!(needed, :t_values)
        push!(needed, :p_values)
    end
    if :ci in req_stats
        push!(needed, :sigma)
        push!(needed, :stderror)
        push!(needed, :t_statistic)
        push!(needed, :ci)
    end

    return needed
end

"""
    get_all_prediction_stats()

    get all the available statistics about the values predicted by a fitted model
"""
get_all_prediction_stats() = Set([:predicted, :residuals, :leverage, :stdp, :stdi, :stdr, :student, :rstudent, :lcli, :ucli, :lclp, :uclp, :press, :cooksd])

get_prediction_stats(req_stats::String) = return get_prediction_stats([req_stats])
get_prediction_stats(req_stats::Vector{String}) = return get_prediction_stats(Symbol.(lowercase.(req_stats)))
get_prediction_stats(::Vector{Any}) = return get_prediction_stats([:none])
get_prediction_stats(::Set{Any}) = return get_prediction_stats(Set([:none]))
get_prediction_stats(req_stats::Set{Symbol}) = return get_prediction_stats(collect(req_stats))

"""
    function get_prediction_stats(req_stats::Vector{Symbol})

    return the list of needed statistics and the statistics that need to be presentd given the list of statistics about the predictions the caller wants.
"""
function get_prediction_stats(req_stats::Vector{Symbol})
    needed = Set([:predicted])
    full = get_all_prediction_stats()
    present = Set([:predicted])
    unique!(req_stats)

    length(req_stats) == 0 && return needed, present
    :none in req_stats && return needed, present
    :all in req_stats && return full, full

    :leverage in req_stats && push!(present, :leverage)
    :residuals in req_stats && push!(present, :residuals)
    if :stdp in req_stats
        push!(needed, :leverage)
        push!(present, :stdp)
    end
    if :stdi in req_stats
        push!(needed, :leverage)
        push!(present, :stdi)
    end
    if :stdr in req_stats
        push!(needed, :leverage)
        push!(present, :stdr)
    end
    if :student in req_stats
        push!(needed, :leverage)
        push!(needed, :residuals)
        push!(needed, :stdr)
        push!(present, :student)
    end
    if :rstudent in req_stats
        push!(needed, :leverage)
        push!(needed, :residuals)
        push!(needed, :stdr)
        push!(needed, :student)
        push!(present, :rstudent)
    end
    if :lcli in req_stats
        push!(needed, :leverage)
        push!(needed, :stdi)
        push!(present, :lcli)
    end
    if :ucli in req_stats
        push!(needed, :leverage)
        push!(needed, :stdi)
        push!(present, :ucli)
    end
    if :lclp in req_stats
        push!(needed, :leverage)
        push!(needed, :stdp)
        push!(present, :lclp)
    end
    if :uclp in req_stats
        push!(needed, :leverage)
        push!(needed, :stdp)
        push!(present, :uclp)
    end
    if :press in req_stats
        push!(needed, :residuals)
        push!(needed, :leverage)
        push!(present, :press)
    end
    if :cooksd in req_stats
        push!(needed, :leverage)
        push!(needed, :stdp)
        push!(needed, :stdr)
        push!(needed, :residuals)
        push!(needed, :student)
        push!(present, :cooksd)
    end

    union!(needed, present)

    return needed, present

end

"""
    function encapsulate_string(s)

    (internal) Only used to encapsulate a string into an array.

    used exclusively to handle the function ```StatsBase.coefnames``` which sometime return an array or when there is only one element the element alone.
"""
function encapsulate_string(s::String)
    return [s]
end

"""
    function encapsulate_string(v)

    (internal) Only used to encapsulate a string into an array.

    used exclusively to handle the function ```StatsBase.coefnames``` which sometime return an array or when there is only one element the element alone.
"""
function encapsulate_string(v::Vector{String})
    return v
end

import Printf
"""
    macro gprintf(fmt::String)

    (internal) used to format with %g

    Taken from message published by user o314 at https://discourse.julialang.org/t/printf-with-variable-format-string/3805/6
"""
macro gprintf(fmt::String)
    :((io::IO, arg) -> Printf.@printf(io, $fmt, arg))
end

"""
    function fmt_pad(s::String, value, pad=0)
    (internal) helper to format and pad string for results display
"""
function fmt_pad(s::String, value, pad=0)
    fmt = @gprintf("%g")
    return rpad(s * sprint(fmt, value), pad)
end

using NamedArrays
"""
    function my_namedarray_print(io::IO, n::NamedArray)

    (internal) Print the NamedArray with the type annotation (on the first line).
"""
function my_namedarray_print(io::IO, n::NamedArray)
    tmpio = IOBuffer()
    show(tmpio, n)
    println(io, split(String(take!(tmpio)), "\n", limit=2)[2])
end

"""
    function helper_print_table(io, title, stats::Vector, stats_name::Vector, updformula)

    (Internal) Convenience function to display a table of statistics to the user.
"""
function helper_print_table(io::IO, title, stats::Vector, stats_name::Vector, updformula)
    println(io, "\n$title")
    todelete = [i for (i, v) in enumerate(stats) if isnothing(v)]
    deleteat!(stats, todelete)
    deleteat!(stats_name, todelete)
    m_all_stats = reduce(hcat, stats)
    if m_all_stats isa Vector
        m_all_stats = reshape(m_all_stats, length(m_all_stats), 1)
    end
    na = NamedArray(m_all_stats)
    setnames!(na, encapsulate_string(string.(StatsBase.coefnames(updformula.rhs))), 1)
    setnames!(na, encapsulate_string(string.(stats_name)), 2)
    setdimnames!(na, ("Terms", "Stats"))
    my_namedarray_print(io, na)
end


text:
using Colors
using CoordinateTransformations
using GeometryBasics
using MeshCat
using MeshIO
using Rotations

function default_background!(vis)
    setvisible!(vis["/Background"], true)
    setprop!(vis["/Background"], "top_color", RGBA(1.0, 1.0, 1.0, 1.0))
    setprop!(vis["/Background"], "bottom_color", RGBA(1.0, 1.0, 1.0, 1.0))
    setvisible!(vis["/Axes"], false)
end

function visualize_double_pendulum!(vis, model, x;
      r = 0.025, Δt = 0.1)

	default_background!(vis)

	# kinematics
	function k1com(model::DoublePendulum, q)
		θ1, θ2 = q

		[0.5 * model.l1 * sin(θ1);
		 -0.5 * model.l1 * cos(θ1)]
	end

	function k1ee(model::DoublePendulum, q)
		θ1, θ2 = q

		[model.l1 * sin(θ1);
		 -model.l1 * cos(θ1)]
	end

	function k2com(model::DoublePendulum, q)
		θ1, θ2 = q

		[model.l1 * sin(θ1) + 0.5 * model.l2 * sin(θ1 + θ2);
		 -model.l1 * cos(θ1) - 0.5 * model.l2 * cos(θ1 + θ2)]
	end

	function k2ee(model::DoublePendulum, q)
		θ1, θ2 = q

		[model.l1 * sin(θ1) + model.l2 * sin(θ1 + θ2);
		 -model.l1 * cos(θ1) - model.l2 * cos(θ1 + θ2)]
	end

	link1 = Cylinder(Point3f0(0.0,0.0,0.0), Point3f0(0.0, 0.0, model.l1),
		convert(Float32, 0.015))
	setobject!(vis["link1"], link1,
		MeshPhongMaterial(color = RGBA(0.0, 0.0, 0.0, 1.0)))

	link2 = Cylinder(Point3f0(0.0,0.0,0.0), Point3f0(0.0, 0.0, model.l1),
		convert(Float32, 0.015))
	setobject!(vis["link2"], link2,
		MeshPhongMaterial(color = RGBA(0.0, 0.0, 0.0, 1.0)))

	T = length(x)

	anim = MeshCat.Animation(convert(Int, floor(1.0 / Δt)))

	for t = 1:T
		q = x[t][1:2]

		MeshCat.atframe(anim, t) do

			_k1ee = k1ee(model, q)
			p1ee = [_k1ee[1], 0.0, _k1ee[2]]

			_k2ee = k2ee(model, q)
			p2ee = [_k2ee[1], 0.0, _k2ee[2]]

			settransform!(vis["link1"], cable_transform([0.0; 0.0; 0.0], p1ee))
			settransform!(vis["link2"], cable_transform(p1ee, p2ee))
		end
	end

	settransform!(vis["/Cameras/default"],
	    compose(Translation(0.0, 0.0, -1.0), LinearMap(RotZ(-pi / 2.0))))

	MeshCat.setanimation!(vis, anim)
end

function cable_transform(y, z)
    v1 = [0.0, 0.0, 1.0]
    v2 = y[1:3,1] - z[1:3,1]
    normalize!(v2)
    ax = cross(v1, v2)
    ang = acos(v1'*v2)
    R = AngleAxis(ang, ax...)

    if any(isnan.(R))
        R = I
    else
        nothing
    end

    compose(Translation(z), LinearMap(R))
end


text:
using PkgSRA
using Plots, Random, Distributions, Interpolations
using StatsPlots, StatsBase
pyplot()
Random.seed!(11)
using BSON: @save
include("../scripts/paper_settings.jl");

################################################################################
#####
##### Lotka-Volterra example: get data
#####
include(EXAMPLE_FOLDERNAME*"example_lotkaVolterra.jl")

# Define the multivariate forcing function
function local_make_data(num_ctr)
    # num_ctr = 7;
    U_starts = rand(1, num_ctr) .* tspan[2]
    U_widths = 0.5;
    amplitude = 2.0
    my_U_func_time2(t) = U_func_time(t, u0,
                            U_widths, U_starts,
                            F_dim=1,
                            amplitude=amplitude)
    # Get data
    sol = solve_lv_system(U_func_time=my_U_func_time2)
    dat = Array(sol)

    # Also get truth
    true_grad = core_dyn_true(dat)
    U_true = zeros(size(dat))
    for (i, t) in enumerate(ts)
        U_true[:,i] = my_U_func_time2(t)
    end

    ## Also get baseline true/ideal cases
    # Uncontrolled
    dat_raw = Array(solve_lv_system())
    numerical_grad_raw = numerical_derivative(dat_raw, ts)
    true_grad_raw = zeros(size(dat))
    #
    # # Intialize truth object
    this_truth = sra_truth_object(dat_raw, true_grad, U_true, core_dyn_true)

    # End
    return this_truth, dat
end

#####
##### Loop over noise and number of spikes
#####
noise_vals = 0.0:0.05:0.5
num_models = 50
control_signal_vals = 0:4:20
# noise_vals = 0.0:0.05:0.2
# num_models = 10
# control_signal_vals = [0, 5, 10]

## Output values
this_size = (length(noise_vals), length(control_signal_vals), num_models)
all_err = zeros(this_size)
# all_err_deriv = zeros(this_size)
all_err_deriv_subsample = zeros(this_size)
all_naive_err = zeros(this_size)
all_naive_err_deriv = zeros(this_size)

# all_i = zeros(this_size)

# Make model template
function local_param_reset!(this_model)
    variable_names = ["x", "y"]
    optimizer = slstNumber((1, 1)); # Number of terms-based sparse regression
    library = Dict("cross_terms"=>[2], "constant"=>nothing)
    model_template = sindyModel(convert_string2function(library),
        variable_names, optimizer)
    this_model.parameters = get_sra_defaults(false)
    prams = this_model.parameters
    prams.model_template = model_template
    prams.sindy_terms_list = Iterators.product(1:3, 1:3)
    prams.initial_subsampling = true
end

for (i_noise,σ) in enumerate(noise_vals)
    for (i_ctr, num_ctr) in enumerate(control_signal_vals)
    this_truth, dat = local_make_data(num_ctr)
        for i_model in 1:num_models
            # Add noise, THEN take derivatives
            noisy_dat = dat .+ σ.*randn(size(dat))
            noisy_grad = numerical_derivative(noisy_dat, ts)
            # Initialize
            global this_model =
                sra_stateful_object(ts, tspan, dat, u0, noisy_grad)
            # Reset parameters because we are NOT using sindyc
            local_param_reset!(this_model)

            fit_first_model(this_model, 30);
            # print_true_equations(this_truth)
            # print_current_equations(this_model)
            all_naive_err[i_noise, i_ctr, i_model] =
                    calc_coefficient_error(this_model, this_truth)
            all_naive_err_deriv[i_noise, i_ctr, i_model] = rss_sindy_derivs(
                    this_model.sindy_model, noisy_dat, this_truth.true_grad) ./
                    length(dat)

            #################################################################
            ### Iterate
            is_improved, max_iter = true, 10
            while is_improved
                calculate_subsampled_ind(this_model);
                all_models, is_improved = fit_model(this_model);
                # is_improved && print_current_equations(this_model)
                this_model.i > max_iter && break
            end
            println("Finished with $(this_model.i) iterations")

            all_err[i_noise, i_ctr, i_model] =
                calc_coefficient_error(this_model, this_truth)
            # all_err_deriv[i_noise, i_ctr, i_model] = rss_sindy_derivs(
            #         this_model.sindy_model, dat, numerical_grad)
            ind = this_model.subsample_ind
            all_err_deriv_subsample[i_noise, i_ctr, i_model] = rss_sindy_derivs(
                    this_model.sindy_model, noisy_dat[:,ind], this_truth.true_grad[:,ind]) ./
                    (length(ind)*size(dat,1))
            # all_i[i_noise, i_ctr, i_model] = this_model.i
            # println("Final error in coefficients = $(all_err[i_noise, i_model])")
            # println("")
        end
        println("Finished number of control signals $num_ctr")
    end
    println("Finished noise level $σ")
    println("=====================================")
end

coef_norm = sum(core_dyn_true.A.^2)

# heatmap_coef_naive = mean(all_naive_err./coef_norm, dims=3)[:,:,1]
# heatmap_coef_final = mean(all_err./coef_norm, dims=3)[:,:,1]
# heatmap(control_signal_vals, noise_vals, heatmap_coef_final)
#     xlabel!("Controllers")
#     ylabel!("Noise")
#     title!("Final Model Error")
# heatmap(control_signal_vals, noise_vals, heatmap_coef_naive .- heatmap_coef_final)
#     xlabel!("Controllers")
#     ylabel!("Noise")
#     title!("Improvement in error")

# vec_naive, std_naive = mean_and_std(all_naive_err./coef_norm, 2)
# vec_err, std_err = mean_and_std(all_err./coef_norm, 2)
# # Error in derivatives: only look at where there's no control
# vec_deriv, std_deriv = mean_and_std(all_err_deriv_subsample, 2)
# vec_naive_deriv, std_naive_deriv = mean_and_std(all_naive_err_deriv, 2)

#####
##### Save Lotka-Volterra data
#####
this_dat_name = "dat_library_with_noise_lv_"

# save_for_plotting(noise_vals,
#                 vec_naive, std_naive, vec_err, std_err,
#                 vec_deriv, std_deriv, vec_naive_deriv, std_naive_deriv,
#                 this_dat_name)

save_for_plotting(noise_vals, control_signal_vals, coef_norm,
                all_naive_err, all_err,
                all_err_deriv_subsample, all_naive_err_deriv,
                this_dat_name)


text:
# __precompile__()

module MotionPlanning

using Reexport
@reexport using StaticArrays
@reexport using Distances
# import PyPlot; const plt = PyPlot
using Distributions; import Base.eltype; eltype{T}(::Dirichlet{T}) = T    # TODO: submit Distributions PR
using SCS
using Iterators
using NearestNeighbors
using DataStructures

### Includes
include("primitivetypes.jl")
include("utilities/utils.jl")
include("collisioncheckers.jl")
include("goals.jl")
include("nearneighbors.jl")
include("statespaces.jl")
include("problems.jl")
include("sampling.jl")
include("planners.jl")
include("postprocessors.jl")
# include("plotting.jl")

# TODO: test suite!
# TODO: T<:AbstractFloat - either ensure type stability (at least for states), or just go with Float64 everything
# TODO: change up lazy vcat [.,.] syntax in anticipation of v0.4; also consider push!/append!
# TODO: half-baked controls caching - do it or kill it
# TODO: "prototype" (.h-ish) files to specify required methods; julia type system makes OOP a pain
# TODO: FMT re-solving is wonky with different sample sizes
# TODO: Think about parameterizing MPProblem on the StateSpace type (for ease of postprocessing)
# TODO: Requires.jl for conditional dependency loading

end # module


text:
"""
    GlobalPencilArray{T,N} <: AbstractArray{T,N}

Alias for an `OffsetArray` wrapping a [`PencilArray`](@ref).

Unlike `PencilArray`s, `GlobalPencilArray`s take *global* indices, which
in general don't start at 1 for a given MPI process.

The [`global_view`](@ref) function should be used to create a
`GlobalPencilArray` from a `PencilArray`.
"""
const GlobalPencilArray{T,N} = OffsetArray{T,N,A} where {A <: PencilArray}

"""
    global_view(x::PencilArray)

Create an [`OffsetArray`](https://github.com/JuliaArrays/OffsetArrays.jl) of a
`PencilArray` that takes global indices in logical order.
"""
function global_view(x::PencilArray)
    r = range_local(x, LogicalOrder())
    offsets = first.(r) .- 1
    xo = OffsetArray(x, offsets)
    @assert parent(xo) === x  # OffsetArray shouldn't create a copy...
    xo :: GlobalPencilArray
end

Base.similar(x::GlobalPencilArray, ::Type{S}) where {S} =
    global_view(similar(parent(x), S))

# Account for index permutation in global views.
@inline Base._sub2ind(x::GlobalPencilArray, I...) =
    Base._sub2ind(parent(x), (I .- x.offsets)...)


text:
using SurrogateCube
# Define some Parameters
Ncomp = 3 #We want 3 independent components
Nvar  = 5 # we want 10 measured variables
noise = WhiteNoise(0.2) #Each component has a noise with sigma=0.5
means = [ConstantBaseline(0.0),ConstantBaseline(0.0),SineBaseline(5,2.0)] # One of the Conmponents has a seasonal cycle
dists = [EmptyEvent(),EmptyEvent(),CubeEvent(0.2)] # Only one of the components is disturbed with a size of 20% in each dimension
varNoise = WhiteNoise(0.1) # each variable generated has some independent white noise
kb=-1  # We increase the amplitude of the MAC by a factor of 2
kn=0.0  # No disturbance in the variance
ks=0.0  # No disturbance in the mean
#Set dimensions
Nlon = 10
Nlat = 10
Ntime =100
ds = genDataCube(means,noise,dists,varNoise,Ncomp,Nvar,Ntime,Nlat,Nlon,kb,kn,ks);
#Save the datacube to a file
save_cube(ds)
#Plot something
using PyPlot
pygui(true)
f=figure()
#plot(int(d[1,5,5,:][:]))
for i=1:Nvar
plot(ds.variables[:,5,5,i][:])
end

max_stars_repo_name: CAB-LAB/DataFarm

max_issues_repo_name: esa-esdl/DataFarm

max_forks_repo_name: esa-esdl/DataFarm


text:
function refine(P::LazySet, ζ::Number)
    constraints = Vector{HalfSpace{Float64, Vector{Float64}}}()
    for c in constraints_list(P)
        c_tightened = HalfSpace(c.a, c.b - (ζ * norm(c.a, 2)))
        push!(constraints, c_tightened)
    end
    return HPolytope(constraints)
end

function refine!(vertices_sat, P::LazySet, vertices_unsat, ζ::Number,
                 dyn1⁻, dyn2⁻, x1, time, ε; log::Bool=false)
    if !isempty(vertices_unsat)
        while true
            P = refine(P, ζ)
            vertices = vertices_list(P)
            vertices_sat2, vertices_unsat2 =
                check_sufficient_points(dyn1⁻, dyn2⁻, x1, vertices, time, ε; log=log)
            append!(vertices_sat, vertices_sat2)
            if isempty(vertices_unsat2)
                break
            end
            if log
                println("refinement continues")
            end
        end
    end
    return vertices_sat
end

max_stars_repo_name: HySynth/HySynth

max_issues_repo_name: HySynth/HySynth

max_forks_repo_name: HySynth/HySynth


text:
showtail(b::Bool) = b ? "…" : ""
ngaps(t::Array{Int64,2}) = max(0, size(t,1)-2 + (t[end,2] == 0 ? 0 : 1))

function str_trunc(str::String, W::Int64)
  i = 0
  j = firstindex(str)
  for c in eachindex(str)
   i += 1
   if i > W-2
     return str[1:j] * "… "
   end
   j = c
 end
 return(rpad(str, W))
end

function show_str(io::IO, S::Array{String,1}, w::Int, N::Int64)
  for i = 1:N
    print(io, str_trunc(S[i], w))
  end
  println(io, showtail(N<length(S)))
  return nothing
end

function show_t(io::IO, T::Array{Array{Int64,2},1}, w::Int, N::Int64, fs::Array{Float64,1})
  for i = 1:N
    if isempty(T[i])
      s = ""
    elseif fs[i] == 0.0
      s = string(timestamp(T[i][1,2]*μs), " (", size(T[i], 1), " vals)")
    else
      s = string(timestamp(T[i][1,2]*μs), " (", ngaps(T[i]), " gaps)")
    end
    print(io, str_trunc(s, w))
  end
  println(io, showtail(N<length(T)))
  return
end

function mkxstr(N::Int64, X::Union{ Array{AbstractArray{Float64,1},1},
                                    Array{AbstractArray{Float32,1},1},
                                    Array{Union{AbstractArray{Float64,1},
                                          AbstractArray{Float32,1}},1} })

  # Fill matrix of X values
  vx = 5
  X_str = Array{String,2}(undef, vx, N)
  fill!(X_str, "")
  for j = 1:N
    x = getindex(X, j)
    nx = length(x)
    if nx == 0
      X_str[1,j] = "(empty)"
      continue
    else
      if nx < vx
        for i = 1:nx
          X_str[i,j] = @sprintf("%+10.3e", x[i])
        end
      else
        for i = 1:vx-3
          X_str[i,j]    = @sprintf("%+10.3e", x[i])
        end
        X_str[vx-2,j]   = "    ..."
        X_str[vx-1,j]   = @sprintf("%+10.3e", x[nx])
      end
      X_str[vx,j] = string("(nx = ", nx, ")")
    end
  end
  return X_str
end

# Fill matrix of X value strings
function mkxstr(X::FloatArray)
  vx = 5
  X_str = Array{String,2}(undef, vx, 1)
  fill!(X_str, "")
  nx = lastindex(X)
  if nx == 0
    X_str[1,1]      = "(empty)"
  elseif nx < vx
    for i = 1:nx
      X_str[i,1]    = @sprintf("%+10.3e", X[i])
    end
  else
    nx_str          = string(nx)
    for i = 1:vx-3
      X_str[i,1]    = @sprintf("%+10.3e", X[i])
    end
    X_str[vx-2,1]   = "    ..."
    X_str[vx-1,1]   = @sprintf("%+10.3e", X[nx])
    X_str[vx,1]     = string("(nx = ", nx_str, ")")
  end
  return X_str
end

function show_x(io::IO, X_str::Array{String,2}, w::Int64, b::Bool)
  (vx, N) = size(X_str)

  # Display
  for i = 1:vx
    if i > 1
      print(io, " "^show_os)
    end

    for j = 1:N
      x_str = X_str[i,j]
      L = length(x_str)
      print(io, x_str)
      if (x_str == "(empty)" || x_str == "") && N == 1
        return nothing
      end
      print(io, " "^(w-L))
    end
    print(io, showtail(b))
    if i < vx
      print(io, "\n")
    end
  end
  return nothing
end

function show_conn(io::IO, C::Array{TCPSocket,1})
  print(io, string(sum([isopen(i) for i in C]), " open, ", length(C), " total"))
  if !isempty(C)
    for (m,c) in enumerate(C)
      print(io, "\n")
      if isopen(c)
        (url,port) = getpeername(c)
        print(io, " "^show_os, "(", m, ") ", url, ":", Int(port))
      else
        print(io, " "^show_os, "(", m, ") (closed)")
      end
    end
  end
  return nothing
end

summary(S::GphysData) = string(typeof(S), " with ", S.n, " channel",  S.n == 1 ? "" : "s")
summary(S::GphysChannel) = string(typeof(S), " with ",
  length(S.x), " sample", (length(S.x) == 1 ? "" : "s"), ", gaps: ", ngaps(S.t))

# GphysData
function show(io::IO, S::T) where {T<:GphysData}
  W = max(80, displaysize(io)[2]) - show_os
  w = min(W, 35)
  nc = getfield(S, :n)
  N = min(nc, div(W-1, w))
  M = min(N+1, nc)
  println(io, T, " with ", nc, " channels (", N, " shown)")
  F = fieldnames(T)
  for f in F
    if (f in unindexed_fields) == false
      targ = getfield(S, f)
      t = typeof(targ)
      fstr = uppercase(String(f))
      print(io, lpad(fstr, show_os-2), ": ")
      if t == Array{String,1}
        show_str(io, targ, w, N)
      elseif f == :notes || f == :misc
        show_str(io, String[string(length(getindex(targ, i)), " entries") for i = 1:M], w, N)
      elseif f == :pha
        show_str(io, String[string(length(getindex(targ, i)), " phases") for i = 1:M], w, N)
      elseif f == :t
        show_t(io, targ, w, N, S.fs)
      elseif f == :x
        x_str = mkxstr(N, getfield(S, :x))
        show_x(io, x_str, w, N<nc)
      else
        show_str(io, String[repr("text/plain", targ[i], context=:compact=>true) for i = 1:M], w, N)
      end
    elseif f == :c
      print(io, "\n", lpad("C", show_os-2), ": ")
      show_conn(io, S.c)
    end
  end
  return nothing
end
show(S::SeisData) = show(stdout, S)

# GphysChannel
function show(io::IO, C::T) where T<:GphysChannel
  W = max(80,displaysize(io)[2]-2)-show_os
  w = min(W, 36)
  nx = length(C.x)
  F = fieldnames(T)

  println(io, T, " with ", nx, " samples")
  for f in F
    targ = getfield(C, f)
    t = typeof(targ)
    fstr = uppercase(String(f))
    print(io, lpad(fstr, show_os-2), ": ")
    if t == String
      println(io, targ)
    elseif (t <: AbstractFloat || t <: InstrumentPosition || t<: InstrumentResponse)
      println(io, repr("text/plain", targ, context=:compact=>true))
    elseif f == :notes
      println(io, string(length(targ), " entries"))
    elseif f == :misc
      println(io, string(length(targ), " entries"))
    elseif f == :pha
      println(io, string(length(targ), " phases"))
    elseif f == :t
      show_t(io, [targ], w, 1, [C.fs])
    else
      x_str = mkxstr(getfield(C, :x))
      show_x(io, x_str, w, false)
    end
  end
  return nothing
end
show(C::GphysChannel) = show(stdout, C)