text:
abstract type GeometricEquation{invType,parType,perType} end

abstract type DifferentialEquation{invType,parType,perType} <: GeometricEquation{invType,parType,perType} end
abstract type DifferentialAlgebraicEquation{invType,parType,perType,secType} <: GeometricEquation{invType,parType,perType} end
abstract type StochasticDifferentialEquation{invType,parType,perType} <: GeometricEquation{invType,parType,perType} end

abstract type AbstractEquationODE{invType,parType,perType} <: GeometricEquation{invType,parType,perType} end
abstract type AbstractEquationDAE{invType,parType,perType,secType} <: DifferentialAlgebraicEquation{invType,parType,perType,secType} end
abstract type AbstractEquationSDE{invType,parType,perType} <: StochasticDifferentialEquation{invType,parType,perType} end
abstract type AbstractEquationPODE{invType,parType,perType} <: GeometricEquation{invType,parType,perType} end
abstract type AbstractEquationPDAE{invType,parType,perType,secType} <: DifferentialAlgebraicEquation{invType,parType,perType,secType} end
abstract type AbstractEquationPSDE{invType,parType,perType} <: StochasticDifferentialEquation{invType,parType,perType} end

const GEinvType{invType,parType,perType} = GeometricEquation{invType,parType,perType} # type alias for dispatch on invariants type parameter
const GEparType{parType,invType,perType} = GeometricEquation{invType,parType,perType} # type alias for dispatch on parameters type parameter
const GEperType{perType,invType,parType} = GeometricEquation{invType,parType,perType} # type alias for dispatch on periodicity type parameter

hasinvariants(::GEinvType{<:NullInvariants}) = false
hasinvariants(::GEinvType{<:Nothing}) = false
hasinvariants(::GEinvType{<:NamedTuple}) = true

hasparameters(::GEparType{<:NullParameters}) = false
hasparameters(::GEparType{<:Nothing}) = false
hasparameters(::GEparType{<:NamedTuple}) = true

hasperiodicity(::GEperType{<:NullPeriodicity}) = false
hasperiodicity(::GEperType{<:Nothing}) = false
hasperiodicity(::GEperType{<:AbstractArray}) = true

const DAEsecType{secType,invType,parType,perType} = DifferentialAlgebraicEquation{invType,parType,perType,secType} # type alias for dispatch on secondary constraint type parameter

hassecondary(::DAEsecType{<:Nothing}) = false
hassecondary(::DAEsecType{<:Callable}) = true

_functions(equ::GeometricEquation) = error("_functions(::GeometricEquation) not implemented for ", typeof(equ), ".")
_solutions(equ::GeometricEquation) = error("_solutions(::GeometricEquation) not implemented for ", typeof(equ), ".")
_functions(equ::GeometricEquation, ::OptionalParameters) = error("_functions(::GeometricEquation, ::OptionalParameters) not implemented for ", typeof(equ), ".")
_solutions(equ::GeometricEquation, ::OptionalParameters) = error("_solutions(::GeometricEquation, ::OptionalParameters) not implemented for ", typeof(equ), ".")

function functions(equ::GeometricEquation)
    if hasvectorfield(equ)
        return _functions(equ)
    else
        return NamedTuple()
    end
end

function functions(equ::GeometricEquation, params::OptionalParameters)
    @assert check_parameters(equ, params)
    if hasvectorfield(equ)
        return _functions(equ, params)
    else
        return NamedTuple()
    end
end

function solutions(equ::GeometricEquation)
    if hassolution(equ)
        return _solutions(equ)
    else
        return NamedTuple()
    end
end

function solutions(equ::GeometricEquation, params::OptionalParameters)
    @assert check_parameters(equ, params)
    if hassolution(equ)
        return _solutions(equ, params)
    else
        return NamedTuple()
    end
end

GeometricBase.invariants(equ::GeometricEquation) = error("invariants(::GeometricEquation) not implemented for ", typeof(equ), ".")
GeometricBase.parameters(equ::GeometricEquation) = error("parameters(::GeometricEquation) not implemented for ", typeof(equ), ".")
GeometricBase.periodicity(equ::GeometricEquation) = error("periodicity(::GeometricEquation) not implemented for ", typeof(equ), ".")

hassolution(::GeometricEquation) = false
hasvectorfield(::GeometricEquation) = false
hasprimary(::GeometricEquation) = false
hassecondary(::GeometricEquation) = false

hasinvariants(::GeometricEquation) = false
hasparameters(::GeometricEquation) = false
hasperiodicity(::GeometricEquation) = false

hashamiltonian(::GeometricEquation) = false
haslagrangian(::GeometricEquation) = false

datatype(equ::GeometricEquation, ics::NamedTuple) = error("datatype(::GeometricEquation, ::NamedTuple) not implemented for ", typeof(equ), ".")
arrtype(equ::GeometricEquation, ics::NamedTuple) = error("arrtype(::GeometricEquation, ::NamedTuple) not implemented for ", typeof(equ), ".")

check_initial_conditions(equ::GeometricEquation, ics::NamedTuple) = error("check_initial_conditions(::GeometricEquation, ::NamedTuple) not implemented for ", typeof(equ), ".")
check_methods(equ::GeometricEquation, tspan, ics, params) = error("check_methods(::GeometricEquation, ::Tuple, ::NamedTuple, ::OptionalParameters) not implemented for ", typeof(equ), ".")

function check_parameters(equ::GeometricEquation, params::NamedTuple)
    typeof(parameters(equ)) <: NamedTuple || return false
    keys(parameters(equ)) == keys(params) || return false
    for key in keys(params)
        haskey(parameters(equ), key) || return false
        typeof(params[key]) <: parameters(equ)[key] || return false
    end
    return true
end

function check_parameters(equ::GeometricEquation, ::NullParameters)
    if parameters(equ) == NullParameters()
        return true
    else
        return false
    end
end

function check_parameters(::GeometricEquation, ::Any)
    error("Parameters are neither a NamedTuple nor NullParameters.")
end


text:
using Test, Images, Colors
using Pyramids

function end_to_end(test_im, T)
    if typeof(T) <: ComplexSteerablePyramid
        scale = 0.5^(1/4);
        pyramid = ImagePyramid(test_im, T, scale=scale, num_orientations=8, max_levels=23, min_size=15, twidth=1.0)
    else
        pyramid = ImagePyramid(test_im, T, max_levels=23, min_size=15)
    end

    test_im_recon = toimage(pyramid)

    return all((test_im_recon .- test_im) .< 0.0002)
end

rand_im = rand(128, 128)
@testset begin
    println("Running end-to-end image comparison test for Complex Steerable Pyramid.")
    @test end_to_end(rand_im, ComplexSteerablePyramid())

    println("Running end-to-end image comparison test for Gaussian Pyramid.")
    @test end_to_end(rand_im, GaussianPyramid())

    println("Running end-to-end image comparison test for Laplacian Pyramid.")
    @test end_to_end(rand_im, LaplacianPyramid())
end

# TODO:
#    Test pyramid copying
#       Test that pyramid levels don't unexpectedly change
#    Test subband editing
#    Test another scale factor

ext: jl


text: sum(map(int,split(string(big(2)^1000),"")))

max_stars_repo_name: cmey/euler

max_issues_repo_name: cmey/euler

max_forks_repo_name: cmey/euler


text:
module PSpaces

import GridInterpolations: RectangleGrid, ind2x, interpolants

export
  PSpace,
  RangeSpace,
  ValuesSpace,
  ind2x,
  ind2x!,
  x2ind

abstract Space

type RangeSpace <: Space

  values::Range

  RangeSpace(range::Range) = new(range)
  RangeSpace{T<:Real}(minval::T, maxval::T) = new(minval:maxval)
  RangeSpace{T<:Real}(minval::T, step::T, maxval::T) = new(minval:step:maxval)

end

type ValuesSpace <: Space

  values::Vector

  function ValuesSpace(vector::Vector)
    if isunique(vector)
      new(vector)
    else
      error("Input values need to be unique.")
    end
  end

end

type PSpace

  cutpoints::Vector{Vector{Any}}
  grid::RectangleGrid

  function PSpace(spaces::Space...)
    cutpoints = [collect(spaces[i].values) for i in 1:length(spaces)]
    spaces = [collect(1.0:length(spaces[i].values)) for i in 1:length(spaces)]
    grid = RectangleGrid(spaces...)
    new(cutpoints, grid)
  end

end

function isunique(vector::Vector)
  for i in 1:length(vector)
    for j in i+1:length(vector)
      if vector[i] == vector[j]
        return false
      end
    end
  end
  return true
end

function ind2x(pspace::PSpace, ind::Int)
  x = Array(Any, length(pspace.cutpoints))
  ind2x!(x, pspace, ind)
  return x
end

function ind2x!(x::Vector, pspace::PSpace, ind::Int)
  gridx = ind2x(pspace.grid, ind)
  grid2x!(x, pspace, gridx)
end

function grid2x!(x::Vector, pspace::PSpace, gridx::Vector{Float64})
  intgridx = map(Int64, gridx)
  for i in 1:length(x)
    x[i] = pspace.cutpoints[i][intgridx[i]]
  end
end

# warning: ensure that |x| is a valid vector, otherwise you'll get junk
function x2ind(pspace::PSpace, x::Vector)
  gridx = [Float64(findfirst(pspace.cutpoints[i], x[i])) for i in 1:length(x)]
  # @assert !in(0, gridx)  # gridx[i] = 0.0 if x[i] isn't a valid variable value
  indices, weights = interpolants(pspace.grid, gridx)
  indices[indmax(weights)]
end

end


text:
using ProjAssistant
@quickactivate

# ------------------------------------------------------------------
@time begin
    using MAT

    import Chemostat
    const ChLP = Chemostat.LP
    const ChSS = Chemostat.SteadyState
    const ChSU = Chemostat.SimulationUtils
    const ChU = Chemostat.Utils

    import Chemostat_Folsom2014
    const ChF = Chemostat_Folsom2014

    const iJR = ChF.iJR904
    const Fd = ChF.FolsomData # experimental data

    using SimTools
    const SimT = SimTools

    using ProgressMeter
    using Plots
    import GR
    !isinteractive() && GR.inline("png")
    import SparseArrays
end

## ----------------------------------------------------------------------------
MODELS_FILE = procdir(iJR, "base_models.bson")
BASE_MODELS = ldat(MODELS_FILE);

## ----------------------------------------------------------------------------
# Biomass medium sensibility
let
    model = BASE_MODELS["max_model"]
    obj_ider = iJR.BIOMASS_IDER
    xi = Fd.val("xi") |> minimum
    intake_info = iJR.load_base_intake_info()
    results = Dict()
    factors = 0.0:0.01:1.0
    prog = Progress(length(factors) * length(intake_info))
    for (exch, dat) in intake_info
        c0 = dat["c"]
        res = get!(results, exch, [])
        for f in factors
            dat["c"] = f * c0
            ChSS.apply_bound!(model, xi, intake_info; emptyfirst = true)
            growth = try
                    fbaout = ChLP.fba(model, obj_ider, iJR.COST_IDER)
                    ChU.av(model, fbaout, obj_ider)
                catch err; 0.0 end
            push!(res, growth)

            next!(prog; showvalues = [
                    (:exch, exch),
                    (:f, f),
                    (:c0, c0),
                    (:c, dat["c"]),
                    (:growth, growth)
                ]
            )
        end
        dat["c"] = c0
    end
    finish!(prog)

    p = plot(
            title = "Biomass medium sensivility",
            xlabel = "fraction of initial conc",
            ylabel = "growth"
        )
    sresults = sort(collect(results); by = (x) -> sum(x[2]))
    lcount = 4
    for (exch, res) in sresults
        lb_ = lcount > 0 ? exch : ""
        plot!(p, factors, res; label = lb_, lw = 3)
        lcount -= 1
    end
    sfig(p,
        @fileid, "medium_sesitivity_study", ".png"
    )
end

## ----------------------------------------------------------------------------
# Checking fba_obj_val < exp_obj_val
let
    to_map = Fd.val("D") |> enumerate
    for (exp, D) in to_map

        model = BASE_MODELS["fva_models"][exp]

        fbaout = ChLP.fba(model, iJR.BIOMASS_IDER, iJR.COST_IDER);
        fba_obj_val = ChU.av(model, fbaout, iJR.BIOMASS_IDER)
        fba_obj_val = ChU.av(model, fbaout, iJR.BIOMASS_IDER)
        fba_ex_glc_val = ChU.av(model, fbaout, iJR.EX_GLC_IDER)
        fba_ex_glc_b = ChU.bounds(model, iJR.EX_GLC_IDER)
        exp_obj_val = Fd.val("D", exp)

        ChU.tagprintln_inmw("FBA SOLUTION",
            "\nobj_ider:                ", iJR.BIOMASS_IDER,
            "\nfba fba_ex_glc_val:      ", fba_ex_glc_val,
            "\nfba fba_ex_glc_b:        ", fba_ex_glc_b,
            "\nfba obj_val:             ", fba_obj_val,
            "\nexp obj_val:             ", exp_obj_val,
            "\ncost_ider:               ", iJR.COST_IDER,
            "\nfba cost_val:            ", ChU.av(model, fbaout, iJR.COST_IDER),
            "\n\n"
        )
        (fba_obj_val < exp_obj_val) && @warn "fba objval < exp objval" fba_obj_val exp_obj_val
    end
end

## ----------------------------------------------------------------------------
# Find cGLC that fit experimental growth
let
    to_map = Fd.val("D") |> enumerate
    for (Di, D) in to_map
        Kd_cGLC = Fd.val(:cGLC, Di)
        Kd_growth = Fd.val(:D, Di)
        xi = Fd.val(:xi, Di)
        model = BASE_MODELS["max_model"]
        intake_info = iJR.load_base_intake_info()
        for (exch, info) in intake_info
            info["c"] = iJR.MAX_CONC # open medium
        end

        function work_fun(cGLC)
            ## Open intakes except Glucose
            intake_info[iJR.EX_GLC_IDER]["c"] = first(cGLC)

            # impose constraint
            ChSS.apply_bound!(model, xi, intake_info; emptyfirst = true)

            ## fba
            fbaout = ChLP.fba(model, iJR.BIOMASS_IDER, iJR.COST_IDER)
            fba_growth = ChU.av(model, fbaout, iJR.BIOMASS_IDER)
            return [fba_growth]
        end

        cGLC = SimT.grad_desc(work_fun; x0 = [Kd_cGLC], x1 = [Kd_cGLC * 0.9], th = 1e-5,
            C = [Kd_cGLC * 0.1], target = [Kd_growth], maxiters = 500) |> first
        @info "Results" Di Kd_cGLC cGLC
    end

end

## ----------------------------------------------------------------------------
# Testing scaled model
let
    model = BASE_MODELS["max_model"]
    fbaout = ChLP.fba(model, iJR.BIOMASS_IDER, iJR.COST_IDER);
    ChU.tagprintln_inmw("FBA SOLUTION",
        "\nobj_ider:         ", iJR.BIOMASS_IDER,
        "\nsize:             ", size(model),
        "\nfba obj_val:      ", ChU.av(model, fbaout, iJR.BIOMASS_IDER),
        "\ncost_ider:        ", iJR.COST_IDER,
        "\nfba cost_val:     ", ChU.av(model, fbaout, iJR.COST_IDER),
        "\n\n"
    )
    model = ChU.well_scaled_model(model, 100.0; verbose = false)
    fbaout = ChLP.fba(model, iJR.BIOMASS_IDER, iJR.COST_IDER);
    ChU.tagprintln_inmw("FBA SOLUTION",
        "\nobj_ider:         ", iJR.BIOMASS_IDER,
        "\nsize:             ", size(model),
        "\nfba obj_val:      ", ChU.av(model, fbaout, iJR.BIOMASS_IDER),
        "\ncost_ider:        ", iJR.COST_IDER,
        "\nfba cost_val:     ", ChU.av(model, fbaout, iJR.COST_IDER),
        "\n\n"
    )
end


text:
# ---
# title: 1. Two Sum
# id: problem1
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Array, Hash Table
# link: <https://leetcode.com/problems/two-sum/description/>
# hidden: true
# ---
#
# Given an array of integers `nums` and an integer `target`, return _indices of
# the two numbers such that they add up to`target`_.
#
# You may assume that each input would have **_exactly_ one solution**, and you
# may not use the _same_ element twice.
#
# You can return the answer in any order.
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [2,7,11,15], target = 9
#     Output: [0,1]
#     Output: Because nums[0] + nums[1] == 9, we return [0, 1].
#
#
# **Example 2:**
#
#
#
#     Input: nums = [3,2,4], target = 6
#     Output: [1,2]
#
#
# **Example 3:**
#
#
#
#     Input: nums = [3,3], target = 6
#     Output: [0,1]
#
#
#
#
# **Constraints:**
#
#   * `2 <= nums.length <= 103`
#   * `-109 <= nums[i] <= 109`
#   * `-109 <= target <= 109`
#   * **Only one valid answer exists.**
#
#
## @lc code=start
using LeetCode

function two_sum(nums::Vector{Int}, target::Int)::Union{Nothing,Tuple{Int,Int}}
    seen = Dict{Int,Int}()
    for (i, n) in enumerate(nums)
        m = target - n
        if haskey(seen, m)
            return seen[m], i
        else
            seen[n] = i
        end
    end
end
## @lc code=end


text:
using Pkg

cd("Julia")
Pkg.activate(".")
Pkg.instantiate()

using Revise

include("01_wrangling.jl")
include("../src/tweet_helpers.jl")
include("../src/tweet_graphing.jl")

using Chain
using CairoMakie
using AlgebraOfGraphics
using DataFrames
using DataFramesMeta
using Dates
using Graphs
using SparseArrays
using DataStructures
using MetaGraphs
using GraphPlot

set_aog_theme!()

function _plot_stats(df, date_func, dir)
    save(dir * "/graph1.png", wrap_in_makie(generate_unique_users(df, date_func), ["Day of year", "N unique users", "Daily unique users"]), px_per_unit = 3)
    save(dir * "/graph2.png", graph2(df, date_func), px_per_unit = 3)
    save(dir * "/graph3.png", wrap_in_makie(generate_unique_user_rate(df, dayofyear), ["Day of year", "N unique users per min", "Daily unique user rate"], true), px_per_unit = 3)
    save(dir * "/graph4.png", graph4(df, date_func), px_per_unit = 3)
    save(dir * "/graph5.png", graph5(df, date_func), px_per_unit = 3)
end

function gen_stat_plots(date_func, out_dir)
    dfs = read_files("Data")
    scandals = read_names("Data")

    for s in scandals
        if !isdir(out_dir * "/" * s)
            mkdir(out_dir * "/" * s)
        end
    end

    for i = 1:length(dfs)
        _plot_stats(dfs[i], date_func, out_dir * "/" * scandals[i])
    end
end

function _plot_net_met(df, met_funcs, per, cum, dir)
    g = generate_reply_quote_graph(df)
    mg = add_meta_reply_quote_graph(df, g)
    for met_func in met_funcs
        save(dir * "/" * string(Symbol(met_func)) * ".png", plot_graph_metric(mg, met_func, per, Date.(ZonedDateTime.(df.created_at))[1], Date.(ZonedDateTime.(df.created_at))[size(df,1)], cum), px_per_unit = 3)
    end
end

function gen_net_plots(met_funcs, per, cum, out_dir)
    dfs = read_files("Data")
    scandals = read_names("Data")

    for s in scandals
        if !isdir(out_dir * "/" * s)
            mkdir(out_dir * "/" * s)
        end
    end

    for i = 1:length(dfs)
        _plot_net_met(dfs[i], met_funcs, per, cum, out_dir * "/" * scandals[i])
    end
end

function graph_metric(x)
    return mean(degree(x))
end

functions = [graph_metric]
gen_net_plots(functions, Day(1), false, "Plots")

gen_stat_plots(dayofyear, "Plots")


text:
const PROJ_IDER = "EColiCore"
const BIOMASS_IDER = "BIOMASS_Ecoli_core_w_GAM"
const ATPM_IDER = "ATPM"
const MAX_ABS_BOUND = 100


text:
using CSV
using DataFrames
using Tables

function test_corkendall_fromfile()
    corkendall_fromfile("c:/temp/example.csv", "c:/temp/examplecor.csv", true, true, true, true)
end

function the_mother_of_all_tests()
    inputfile = "c:/temp/1044rows11214cols.csv"
    outputfile = "c:/temp/didthatwork.csv"
    @time corkendall_fromfile(inputfile, outputfile, true, true, true, true)
end

function corkendall_fromfile(inputfile::String, outputfile::String, inputhasheaderrow::Bool, inputhasheadercol::Bool,
                            outputhasheaderrow::Bool, outputhasheadercol::Bool)

    header = inputhasheaderrow ? 1 : 0
    drop = inputhasheadercol ? [1] : [0]

    filedata = CSV.File(inputfile;header=header,drop=drop)

    println("file '$inputfile' read into memory")
    names = CSV.getnames(filedata)

    # convert to a matrix, as required by corkendall
    data = Tables.matrix(filedata)

    println("converted to matrix")
    filedata = nothing
    GC.gc()
    println("freed memory first time")

    # @time res = corkendallthreads_v8(data)
    @time res = corkendall(data, skipmissing=:pairwise)

    datatowrite = DataFrame(res, names)
    res = nothing
    GC.gc()
    println("freed memory second time")

    if outputhasheadercol
        insertcols!(datatowrite, 1, Symbol("") => String.(names))
    end

    CSV.write(outputfile, datatowrite, writeheader=outputhasheaderrow)

end


text:
export samplecdf,normalizedcumsum!

##bisection inverse

bisectioninv{S,T}(f::Fun{S,T},x::Real;opts...) = first(bisectioninv(f,[x];opts...))

function bisectioninv{S,T}(f::Fun{S,T},x::Float64;numits::Int=47)
    d=domain(f)
    a = first(d);b = last(d)

    for k=1:numits  #TODO: decide 47
        m=.5*(a+b)
        val = f(m)

            (val<= x) ? (a = m) : (b = m)
    end
    .5*(a+b)
end

bisectioninv{S,T}(f::Fun{S,T},x::Vector;opts...)=Float64[bisectioninv(f,xx;opts...) for xx in x]

## Clenshaw bisection

function chebbisectioninv(c::Vector{Float64},x::Float64;numits::Int=47)
    a = -1.;b = 1.

    C=Chebyshev()
    for k=1:numits  #TODO: decide 47
        m=.5*(a+b)
        val = clenshaw(C,c,m)

            (val<= x) ? (a = m) : (b = m)
    end
    .5*(a+b)
end

chebbisectioninv(c::Vector{Float64},xl::Vector{Float64})=(n=length(xl);chebbisectioninv(c,xl,ClenshawPlan(Float64,Chebyshev(),length(c),n)))
function chebbisectioninv{D<:Domain}(c::Vector{Float64},xl::Vector{Float64},plan::ClenshawPlan{Chebyshev{D},Float64})
    n = length(xl)
    a = -ones(n)
    b = ones(n)

    for k=1:47  #TODO: decide 47
        m=.5*(a+b);
        vals = clenshaw(c,m,plan)

        for j = 1:n
            (vals[j] <= xl[j]) ? (a[j] = m[j]) : (b[j] = m[j])
        end
    end
    m=.5*(a+b)
end

#here, xl is vector w/ length == #cols of c
chebbisectioninv(c::Array{Float64,2},xl::Vector{Float64})=(n=length(xl);chebbisectioninv(c,xl,ClenshawPlan(Float64,Chebyshev(),size(c,1),n)))
function chebbisectioninv{D<:Domain}(c::Array{Float64,2},xl::Vector{Float64},plan::ClenshawPlan{Chebyshev{D},Float64})
    @assert size(c)[2] == length(xl)

    n = length(xl)
    a = -ones(n)
    b = ones(n)

    for k=1:47  #TODO: decide 47
        m=.5*(a+b);
        vals = clenshaw(c,m,plan)

        for j = 1:n
            (vals[j] <= xl[j]) ? (a[j] = m[j]) : (b[j] = m[j])
        end
    end
    m=.5*(a+b)
end

for TYP in (:Vector,:Float64)
    @eval begin
        bisectioninv{SP<:Chebyshev}(cf::Fun{SP,Float64},x::$TYP;opts...)=fromcanonical(cf,chebbisectioninv(coefficients(cf),x;opts...))
#        bisectioninv{SP<:LineSpace}(cf::Fun{SP,Float64},x::$TYP;opts...)=fromcanonical(cf,chebbisectioninv(coefficients(cf),x;opts...))
    end
end

##normalized cumsum

function normalizedcumsum(f::Fun)
    cf = cumsum(f)
    cf = cf/last(cf)

    cf
end

function subtract_zeroatleft!(f::Vector{Float64})
    for k=2:length(f)
        @inbounds f[1] += (-1.)^k.*f[k]
    end

    f
end

function subtract_zeroatleft!(f::Array{Float64,2})
    for k=2:size(f)[1],j=1:size(f)[2]
        @inbounds f[1,j] += (-1.)^k.*f[k,j]
    end

    f
end

function multiply_oneatright!(f::Vector{Float64})
    val=0.
    for k=1:length(f)
        val+=f[k]
    end

    val=1./val

    for k=1:length(f)
        @inbounds f[k] *= val
    end

    f
end

function multiply_oneatright!(f::Array{Float64,2})

    for j=1:size(f)[2]
        val=0.
        for k=1:size(f)[1]
            val+=f[k,j]
        end

        val=1./val

        for k=1:size(f)[1]
            @inbounds f[k,j] *= val
        end
    end

    f
end

function chebnormalizedcumsum!(f)
    ultraconversion!(f)
    ultraint!(f)
    subtract_zeroatleft!(f)
    multiply_oneatright!(f)
end

# For RandomMatrices compatibility
normalizedcumsum!(f::Vector{Float64})=chebnormalizedcumsum!(f)

## Sampling

sample(f::Fun,n::Integer)=samplecdf(normalizedcumsum(f),n)

samplecdf(cf::Fun,n::Integer)=bisectioninv(cf,rand(n))

sample(f::Fun)=sample(f,1)[1]
samplecdf(f::Fun)=samplecdf(f,1)[1]
samplecdf(v::Vector)=chebbisectioninv(v,rand())

##2D sample

sample{TS<:AbstractProductSpace}(f::Fun{TS},k::Integer)=sample(ProductFun(f),k)

function sample(f::LowRankFun,n::Integer)
    rx=sample(sum(f,2),n)
    fA=evaluate(f,rx,:)
    ry=map(sample,fA)
    [rx ry]
end

function sample{C<:Chebyshev}(f::LowRankFun{C,C,TensorSpace{Tuple{C,C},RealBasis,2},Float64},n::Integer)
    ry=sample(sum(f,1),n)
    fA=evaluate(f.A,ry)
    CB=coefficients(f.B)
    AB=CB*fA
    chebnormalizedcumsum!(AB)
    rx=chebbisectioninv(AB,rand(n))
  [fromcanonical(domain(f,1),rx) ry]
end

sample(f::ProductFun,n)=sample(LowRankFun(f),n)
sample(f::MultivariateFun)=sample(f,1)[1,:]

## Special spaces

# Rays may be schwartz at right endpoint so we project
function sample{SS,DD<:Ray}(f::Fun{JacobiWeight{SS,DD},Float64},n::Integer)
    if space(f).β == 0
        samplecdf(normalizedcumsum(f),n)
    else
        sample(Fun(x->f(x),JacobiWeight(space(f).α,1,domain(f))),n)
    end
end

# Line/Ray have unbounded endpoints so we map
# for TYP in (:Line,:Ray)
#     @eval bisectioninv{SS,DD<:$TYP,TT}(f::Fun{MappedSpace{SS,DD,TT},Float64},x::Vector)=fromcanonical(f,
#                                                                                bisectioninv(Fun(f.coefficients,space(f).space),x))
# end

# function sample{SS}(f::LowRankFun{LineSpace{SS},LineSpace{SS},TensorSpace{Tuple{LineSpace{SS},LineSpace{SS}},RealBasis,2},Float64},n::Integer)
#     cf=normalizedcumsum(sum(f,1))
#     CB=coefficients(map(cumsum,f.B))
#
#     ry=samplecdf(cf,n)
#     fA=evaluate(f.A,ry)
#     CBfA=CB*fA  #cumsums at points
#     multiply_oneatright!(CBfA)
#
#     rx=fromcanonical(first(f.B),chebbisectioninv(CBfA,rand(n)))
#
#     [rx ry]
# end


text:
# Concrete interpreter semantics

"Concrete PDDL interpreter."
struct ConcreteInterpreter <: Interpreter end

include("interface.jl")
include("satisfy.jl")
include("initstate.jl")
include("goalstate.jl")
include("update.jl")


text:
force_scale = 2.0 #50
tid = 200
fiber_size = 1
porder = 2
if length(ARGS) == 4
   global tid = parse(Int64, ARGS[1])
   global force_scale = parse(Float64, ARGS[2])
   global fiber_size = parse(Int64, ARGS[3])
   global porder = parse(Int64, ARGS[4])
end
printstyled("force_scale=$force_scale, tid=$tid, fiber_size=$fiber_size, porder=$porder\n", color=:green)
include("CommonFuncs.jl")
if Sys.MACHINE=="x86_64-apple-darwin18.6.0"
    matplotlib.use("macosx")
end
np = pyimport("numpy")

testtype = "PlaneStressIncompressibleRivlinSaunders"
#https://classes.engineering.wustl.edu/2009/spring/mase5513/abaqus/docs/v6.6/books/bmk/default.htm?startat=ch01s01ach08.html
prop = Dict("name"=> testtype, "rho"=> 0.800, "C1"=>0.1863e6, "C2"=>0.00979e6)
#prop = Dict("name"=> testtype, "rho"=> 8.0, "C1"=>0.1863e1, "C2"=>-0.00979e1)  #g/dm^2
# testtype = "PlaneStress"
# prop = Dict("name"=> testtype, "rho"=> 8000.0, "E"=> 200e+9, "nu"=> 0.45)

T = 0.2  #s
NT = 200

nxc, nyc = 10,5
nx, ny =  nxc*fiber_size, nyc*fiber_size
Lx, Ly = 0.1, 0.05 #m
#Lx, Ly = 1.0, 0.5 #dm
nodes, EBC, g, gt, FBC, fext, ft, npoints, node_to_point = BoundaryCondition(tid, nx, ny, porder, Lx, Ly; force_scale=force_scale)
elements = []
for j = 1:ny
    for i = 1:nx
        n = (nx*porder+1)*(j-1)*porder + (i-1)porder+1
        #element (i,j)
        if porder == 1
            #   4 ---- 3
            #
            #   1 ---- 2

            elnodes = [n, n + 1, n + 1 + (nx + 1), n + (nx + 1)]
        elseif porder == 2
            #   4 --7-- 3
            #   8   9   6
            #   1 --5-- 2
            elnodes = [n, n + 2, n + 2 + 2*(2*nx+1),  n + 2*(2*nx+1), n+1, n + 2 + (2*nx+1), n + 1 + 2*(2*nx+1), n + (2*nx+1), n+1+(2*nx+1)]
        else
            error("polynomial order error, porder= ", porder)
        end

        coords = nodes[elnodes,:]
        push!(elements,FiniteStrainContinuum(coords,elnodes, prop, 3))
    end
end

ndofs = 2
domain = Domain(nodes, elements, ndofs, EBC, g, FBC, fext)
setGeometryPoints!(domain, npoints, node_to_point)

state = zeros(domain.neqs)
∂u = zeros(domain.neqs)
globdat = GlobalData(state,zeros(domain.neqs), zeros(domain.neqs),∂u, domain.neqs, gt, ft)

assembleMassMatrix!(globdat, domain)
updateStates!(domain, globdat)

Δt = T/NT

ρ_oo = 0.0
αm = (2*ρ_oo - 1)/(ρ_oo + 1)
αf = ρ_oo/(ρ_oo + 1)

for i = 1:NT
    @info i, "/" , NT

    solver = NewmarkSolver(Δt, globdat, domain, αm, αf, 1e-4, 1e-6, 10)
end

# error()
if !isdir("$(@__DIR__)/Data")
    mkdir("$(@__DIR__)/Data")
end

if !isdir("$(@__DIR__)/Debug")
    mkdir("$(@__DIR__)/Debug")
end

if !isdir("$(@__DIR__)/Data/order$porder")
    mkdir("$(@__DIR__)/Data/order$porder")
end
if !isdir("$(@__DIR__)/Debug/order$porder")
    mkdir("$(@__DIR__)/Debug/order$porder")
end

write_data("$(@__DIR__)/Data/order$porder/$(tid)_$(force_scale)_$(fiber_size).dat", domain)
@save "Data/order$porder/domain$(tid)_$(force_scale)_$(fiber_size).jld2" domain

close("all")
visσ(domain)
axis("equal")
savefig("Debug/order$porder/terminal$(tid)_$(force_scale)_$(fiber_size).png")

close("all")
ux = [reshape(domain.history["state"][i][1:(nx*porder+1)*(ny*porder+1)], ny*porder+1, nx*porder+1)[end,end] for i = 1:length(domain.history["state"])]
plot(ux)
savefig("Debug/order$porder/ux$(tid)_$(force_scale)_$(fiber_size).png")

close("all")
uy = [reshape(domain.history["state"][i][(nx*porder+1)*(ny*porder+1)+1:end], ny*porder+1, nx*porder+1)[end,end] for i = 1:length(domain.history["state"])]
plot(uy)
savefig("Debug/order$porder/uy$(tid)_$(force_scale)_$(fiber_size).png")


text:
mutable struct W_ACRModel <: PM.AbstractWRModel PM.@pm_fields end
mutable struct WU_ACRModel <: PM.AbstractWRModel PM.@pm_fields end

function build_acopf_with_free_lines(pm::W_ACRModel)
    build_opf_mod(pm)
    nw = PM.nw_id_default
    vr = PM.var(pm, nw)[:vr] = JuMP.@variable(pm.model, [i in PM.ids(pm, nw, :bus)], base_name="$(nw)_vr", start = PM.comp_start_value(PM.ref(pm, nw, :bus, i), "vr_start", 1.0))
    vi = PM.var(pm, nw)[:vi] = JuMP.@variable(pm.model, [i in PM.ids(pm, nw, :bus)], base_name="$(nw)_vi", start = PM.comp_start_value(PM.ref(pm, nw, :bus, i), "vi_start"))
    for (i, bus) in PM.ref(pm, nw, :bus)
        JuMP.set_lower_bound(vr[i], -bus["vmax"])
        JuMP.set_upper_bound(vr[i],  bus["vmax"])
        JuMP.set_lower_bound(vi[i], -bus["vmax"])
        JuMP.set_upper_bound(vi[i],  bus["vmax"])
    end

    w = PM.var(pm, :w)
    for (i, bus) in PM.ref(pm, nw, :bus)
        JuMP.@constraint(pm.model, w[i] == vr[i]^2 + vi[i]^2)
    end

    wr = PM.var(pm, :wr)
    wi = PM.var(pm, :wi)
    for (i,j) in PM.ids(pm, :buspairs)
        JuMP.@constraint(pm.model, wr[(i,j)] == vr[i] * vr[j] + vi[i] * vi[j])
        JuMP.@constraint(pm.model, wi[(i,j)] == vi[i] * vr[j] - vr[i] * vi[j])
    end
end

function build_acopf_with_free_lines(pm::WU_ACRModel)
    build_opf_mod(pm)
    nw = PM.nw_id_default
    vr = PM.var(pm, nw)[:vr] = JuMP.@variable(pm.model, [i in PM.ids(pm, nw, :bus)], base_name="$(nw)_vr", start = PM.comp_start_value(PM.ref(pm, nw, :bus, i), "vr_start", 1.0))
    vi = PM.var(pm, nw)[:vi] = JuMP.@variable(pm.model, [i in PM.ids(pm, nw, :bus)], base_name="$(nw)_vi", start = PM.comp_start_value(PM.ref(pm, nw, :bus, i), "vi_start"))
    for (i, bus) in PM.ref(pm, nw, :bus)
        JuMP.set_lower_bound(vr[i], -bus["vmax"])
        JuMP.set_upper_bound(vr[i],  bus["vmax"])
        JuMP.set_lower_bound(vi[i], -bus["vmax"])
        JuMP.set_upper_bound(vi[i],  bus["vmax"])
    end

    w = PM.var(pm, :w)
    ur = PM.var(pm, nw)[:ur] = JuMP.@variable(pm.model, [i in PM.ids(pm, nw, :bus)], base_name="$(nw)_ur", start = 1.0)
    ui = PM.var(pm, nw)[:ui] = JuMP.@variable(pm.model, [i in PM.ids(pm, nw, :bus)], base_name="$(nw)_ui", start = 1.0)
    for (i, bus) in PM.ref(pm, nw, :bus)
        JuMP.@constraint(pm.model, w[i] == ur[i] + ui[i])
        JuMP.@constraint(pm.model, ur[i] == vr[i]^2)
        JuMP.@constraint(pm.model, ui[i] == vi[i]^2)
        JuMP.@constraint(pm.model, ur[i] <= +bus["vmax"] * vr[i])
        JuMP.@constraint(pm.model, ur[i] >= -bus["vmax"] * vr[i])
        JuMP.@constraint(pm.model, ui[i] <= +bus["vmax"] * vi[i])
        JuMP.@constraint(pm.model, ui[i] >= -bus["vmax"] * vi[i])

        # SOC relaxation constraints
        # JuMP.@constraint(pm.model, ur[i] >= vr[i]^2)
        # JuMP.@constraint(pm.model, ui[i] >= vi[i]^2)
    end

    urr = PM.var(pm, nw)[:urr] = JuMP.@variable(pm.model, [(i,j) in PM.ids(pm, :buspairs)], base_name="$(nw)_urr", start = 1.0)
    uii = PM.var(pm, nw)[:uii] = JuMP.@variable(pm.model, [(i,j) in PM.ids(pm, :buspairs)], base_name="$(nw)_uii", start = 1.0)
    uir = PM.var(pm, nw)[:uir] = JuMP.@variable(pm.model, [(i,j) in PM.ids(pm, :buspairs)], base_name="$(nw)_uir", start = 1.0)
    uri = PM.var(pm, nw)[:uri] = JuMP.@variable(pm.model, [(i,j) in PM.ids(pm, :buspairs)], base_name="$(nw)_uri", start = 1.0)
    wr = PM.var(pm, :wr)
    wi = PM.var(pm, :wi)
    for (i,j) in PM.ids(pm, :buspairs)
        JuMP.@constraint(pm.model, wr[(i,j)] == urr[(i,j)] + uii[(i,j)])
        JuMP.@constraint(pm.model, wi[(i,j)] == uir[(i,j)] - uri[(i,j)])

        # Gurobi will branch on these guys...
        JuMP.@constraint(pm.model, urr[(i,j)] == vr[i] * vr[j])
        JuMP.@constraint(pm.model, uii[(i,j)] == vi[i] * vi[j])
        JuMP.@constraint(pm.model, uir[(i,j)] == vi[i] * vr[j])
        JuMP.@constraint(pm.model, uri[(i,j)] == vr[i] * vi[j])

        busi = PM.ref(pm, :bus)[i]
        busj = PM.ref(pm, :bus)[j]

        # McCormick relaxation of the bilinear constraints
        # JuMP.@constraint(pm.model, -urr[(i,j)] - busi["vmax"] * vr[j] - busj["vmax"] * vr[i] <= +busi["vmax"] * busj["vmax"])
        # JuMP.@constraint(pm.model, -urr[(i,j)] + busi["vmax"] * vr[j] + busj["vmax"] * vr[i] <= +busi["vmax"] * busj["vmax"])
        # JuMP.@constraint(pm.model, -urr[(i,j)] + busi["vmax"] * vr[j] - busj["vmax"] * vr[i] >= -busi["vmax"] * busj["vmax"])
        # JuMP.@constraint(pm.model, -urr[(i,j)] - busi["vmax"] * vr[j] + busj["vmax"] * vr[i] >= -busi["vmax"] * busj["vmax"])
        # JuMP.@constraint(pm.model, -uii[(i,j)] - busi["vmax"] * vi[j] - busj["vmax"] * vi[i] <= +busi["vmax"] * busj["vmax"])
        # JuMP.@constraint(pm.model, -uii[(i,j)] + busi["vmax"] * vi[j] + busj["vmax"] * vi[i] <= +busi["vmax"] * busj["vmax"])
        # JuMP.@constraint(pm.model, -uii[(i,j)] + busi["vmax"] * vi[j] - busj["vmax"] * vi[i] >= -busi["vmax"] * busj["vmax"])
        # JuMP.@constraint(pm.model, -uii[(i,j)] - busi["vmax"] * vi[j] + busj["vmax"] * vi[i] >= -busi["vmax"] * busj["vmax"])
        # JuMP.@constraint(pm.model, -uir[(i,j)] - busi["vmax"] * vr[j] - busj["vmax"] * vi[i] <= +busi["vmax"] * busj["vmax"])
        # JuMP.@constraint(pm.model, -uir[(i,j)] + busi["vmax"] * vr[j] + busj["vmax"] * vi[i] <= +busi["vmax"] * busj["vmax"])
        # JuMP.@constraint(pm.model, -uir[(i,j)] + busi["vmax"] * vr[j] - busj["vmax"] * vi[i] >= -busi["vmax"] * busj["vmax"])
        # JuMP.@constraint(pm.model, -uir[(i,j)] - busi["vmax"] * vr[j] + busj["vmax"] * vi[i] >= -busi["vmax"] * busj["vmax"])
        # JuMP.@constraint(pm.model, -uri[(i,j)] - busi["vmax"] * vi[j] - busj["vmax"] * vr[i] <= +busi["vmax"] * busj["vmax"])
        # JuMP.@constraint(pm.model, -uri[(i,j)] + busi["vmax"] * vi[j] + busj["vmax"] * vr[i] <= +busi["vmax"] * busj["vmax"])
        # JuMP.@constraint(pm.model, -uri[(i,j)] + busi["vmax"] * vi[j] - busj["vmax"] * vr[i] >= -busi["vmax"] * busj["vmax"])
        # JuMP.@constraint(pm.model, -uri[(i,j)] - busi["vmax"] * vi[j] + busj["vmax"] * vr[i] >= -busi["vmax"] * busj["vmax"])

        # Simple bounds
        JuMP.@constraint(pm.model, urr[(i,j)] >= -busj["vmax"] * vr[i])
        JuMP.@constraint(pm.model, urr[(i,j)] >= -busi["vmax"] * vr[j])
        JuMP.@constraint(pm.model, urr[(i,j)] <= busj["vmax"] * vr[i])
        JuMP.@constraint(pm.model, urr[(i,j)] <= busi["vmax"] * vr[j])
        JuMP.@constraint(pm.model, uii[(i,j)] >= -busj["vmax"] * vi[i])
        JuMP.@constraint(pm.model, uii[(i,j)] >= -busi["vmax"] * vi[j])
        JuMP.@constraint(pm.model, uii[(i,j)] <= busj["vmax"] * vi[i])
        JuMP.@constraint(pm.model, uii[(i,j)] <= busi["vmax"] * vi[j])
        JuMP.@constraint(pm.model, uir[(i,j)] >= -busj["vmax"] * vi[i])
        JuMP.@constraint(pm.model, uir[(i,j)] >= -busi["vmax"] * vr[j])
        JuMP.@constraint(pm.model, uir[(i,j)] <= busj["vmax"] * vi[i])
        JuMP.@constraint(pm.model, uir[(i,j)] <= busi["vmax"] * vr[j])
        JuMP.@constraint(pm.model, uri[(i,j)] >= -busj["vmax"] * vr[i])
        JuMP.@constraint(pm.model, uri[(i,j)] >= -busi["vmax"] * vi[j])
        JuMP.@constraint(pm.model, uri[(i,j)] <= busj["vmax"] * vr[i])
        JuMP.@constraint(pm.model, uri[(i,j)] <= busi["vmax"] * vi[j])

        # SOC relaxation constraints
        # JuMP.@constraint(pm.model, wr[(i,j)]^2 + wi[(i,j)]^2 <= w[i] * w[j])
    end
end


text:
@testset "testing GoF, multinomial distribution" begin

df = CSV.read(joinpath(dirname(Base.find_package("PhyloNetworks")),"..","examples","buckyCF.csv"))
d0 = readTableCF(df)
d = deepcopy(d0)
net3 = readTopology("((((D:0.4,C:0.4):4.8,((A:0.8,B:0.8):2.2)#H1:2.2::0.7):4.0,(#H1:0::0.3,E:3.0):6.2):2.0,O:11.2);");

@testset "using Pearson statistic" begin
# without optimizing branch lengths
netresult1 = quarnetGoFtest!(net3,df,false; quartetstat=:pearson, correction=:none);
@test length(netresult1) == 6
@test netresult1[2] ≈ 8.589058727506838  # z stat
@test netresult1[1] ≈ 4.384234705965304e-18   # p-value
@test df[!,:p_value] ≈ [1.2435633419824544e-11,0.0009280577186649157,0.0009280577186649157,1.2435633419824544e-11,0.007580817260552542,0.9998389684303506,0.012895416058219087,0.9438090044657973,0.9471103615208266,0.9438090044657973,0.012895416058219087,0.9471103615208266,0.9956830628718893,0.24055486863965628,0.007580817260552542]
@test net3.loglik ≈ 105.30058282649648
end

@testset "LRT statistic for 4-taxon set with small expected values" begin
# without optimizing branch lengths
netresult1 = quarnetGoFtest!(net3,d,false; correction=:none); # modified d: fills in CFs expected from net3
@test netresult1[2] ≈ 6.21966321647047 # z stat
@test netresult1[3] ≈ 1.0 # sigma
@test netresult1[1] ≈ 2.491115579898031e-10  # p-value
@test netresult1[4] ≈ [0.0024449826689709165,0.01496306673600063,0.01496306673600063,0.0024449826689709165,0.04086460431063039,0.9998541057240138,0.1901450501005025,0.8909735618259936,0.9058717147295428,0.8909735618259936,0.1901450501005025,0.9058717147295428,0.9913859984840471,0.3656465603640152,0.04086460431063039]
end

@testset "Qlog statistic for 4-taxon set with small expected values" begin
netresult1 = quarnetGoFtest(d.quartet, QuartetNetworkGoodnessFit.multinom_qlog!);
@test netresult1[1] ≈ 6.21966321647047  # z stat
@test netresult1[2] ≈ [6.188115021844278e-7,.0038172375225786225,.0038172375225786225,6.188115021844278e-7,.016885279326279496,.999846547125979,.054779352829667574,.8909734349100712,.9229400733753065,.8909734349100712,.054779352829667574,.9229400733753065,.9913859984246092,.2975827064317875,.016885279326279496]
end

@testset "check for exceptions" begin
@test_throws ErrorException ticr!(net3,d,false; quartetstat = :Qlog)
@test_throws ErrorException ticr!(net3,d,false; test = :bad);
@test_throws ErrorException quarnetGoFtest!(net3,d,false; quartetstat=:maxCF);
@test_throws ErrorException quarnetGoFtest!(net3,d,false; correction=:foo)
d.quartet[1].ngenes = -1
@test_throws ErrorException quarnetGoFtest!(net3,d,false);
end

d = deepcopy(d0)
@testset "LRT stat with dependence correction" begin
# test of expectedCF_ordered
for i in 1:15 d.quartet[i].qnet.expCF = [i+0.1, i+0.2, i+0.3]; end
expCF, taxa = QuartetNetworkGoodnessFit.expectedCF_ordered(d, net3)
taxa == ["A","B","C","D","E","O"]
expCF ≈ [13.1 13.2 13.3; 3.2 3.1 3.3; 2.2 2.1 2.3; 8.2 8.3 8.1; 10.3 10.2 10.1; 15.2 15.1 15.3; 5.2 5.1 5.3; 12.2 12.3 12.1; 9.3 9.2 9.1; 14.3 14.1 14.2; 4.1 4.3 4.2; 7.3 7.2 7.1; 1.1 1.3 1.2; 11.3 11.2 11.1; 6.1 6.2 6.3]
# test of quarnetGoFtest! with correction
netresult1 = quarnetGoFtest!(net3,d,false; seed=4321, verbose=true, nsim=1, keepfiles=true);
@test netresult1[4] ≈ [0.0024449826689709165,0.01496306673600063,0.01496306673600063,0.0024449826689709165,0.04086460431063039,0.9998541057240138,0.1901450501005025,0.8909735618259936,0.9058717147295428,0.8909735618259936,0.1901450501005025,0.9058717147295428,0.9913859984840471,0.3656465603640152,0.04086460431063039]
@test netresult1[2] ≈ 6.21966321647047 # z stat, uncorrected
@test netresult1[3] != 1.0 # sigma
hldir = filter!(f -> startswith(f, "jl_"), readdir())
@test length(hldir) == 1
@test length(readdir(hldir[1])) == 1
rm(hldir[1], recursive=true)
# similar test, but with multiple processors and other options, and interesting seed:
# for s in 1:200 netresult1 = quarnetGoFtest!(net3,d,false; seed=s, nsim=5);
#                if netresult1[3] > 3.0; @show (s, netresult1); break; end;end;
Distributed.addprocs(2)
# start with: julia -p 2 --project
# or: using Distributed; @everywhere begin; using Pkg; Pkg.activate("."); using PhyloNetworks; end
@everywhere using QuartetNetworkGoodnessFit
netresult1 = quarnetGoFtest!(net3,d,false; seed=182, nsim=5);
Distributed.rmprocs(workers())
@test netresult1[4] ≈ [0.0024449826689709165,0.01496306673600063,0.01496306673600063,0.0024449826689709165,0.04086460431063039,0.9998541057240138,0.1901450501005025,0.8909735618259936,0.9058717147295428,0.8909735618259936,0.1901450501005025,0.9058717147295428,0.9913859984840471,0.3656465603640152,0.04086460431063039]
@test netresult1[2] ≈ 6.21966321647047 # z stat, uncorrected
@test netresult1[3] ≈ 3.405362128771355 # sigma
@test netresult1[6] ≈ vcat(7.4043609719886545, repeat([-0.8885233166386386],4))
end

end


text:
# NOTE: theses tests focus on speciall characters, html entities
# escaping things etc.

@testset "Backslashes" begin # see issue #205
    st = raw"""
        Hello \ blah \ end
        and `B \ c` end and

        ```
        A \ b
        ```

        done
        """

    steps = explore_md_steps(st)
    tokens, = steps[:tokenization]

    # the first two backspaces are detected
    @test tokens[1].ss == "\\" && tokens[1].name == :CHAR_BACKSPACE
    @test tokens[2].ss == "\\" && tokens[2].name == :CHAR_BACKSPACE
    # the third one also
    @test tokens[5].ss == "\\" && tokens[5].name == :CHAR_BACKSPACE

    sp_chars, = steps[:spchars]

    # there's only two tokens left which are the backspaces NOT in the code env
    sp_chars = F.find_special_chars(tokens)
    for i in 1:2
        @test sp_chars[i] isa F.HTML_SPCH
        @test sp_chars[i].ss == "\\"
        @test sp_chars[i].r == "\"
    end

    inter_html, = steps[:inter_html]

    @test isapproxstr(inter_html, "
Hello  ##FDINSERT##  blah  ##FDINSERT##  end and  ##FDINSERT##  end and

 ##FDINSERT##

done
")

    @test isapproxstr(st |> seval, raw"""

Hello \ blah \ end and B \ c end and

A \ b

done

                """)
end

@testset "Backslashes2" begin # see issue #205
    st = raw"""
        Hello \ blah \ end
        and `B \ c` end \\ and

        ```
        A \ b
        ```

        done
        """
    steps = explore_md_steps(st)
    tokens, = steps[:tokenization]
    @test tokens[7].name == :CHAR_LINEBREAK
    h = st |> seval
    @test isapproxstr(st |> seval, raw"""

Hello \ blah \ end
                        and B \ c end
 and

A \ b

done

                        """)
end

@testset "Backtick" begin # see issue #205
    st = raw"""Blah \` etc"""
    @test isapproxstr(st |> seval, "
Blah ` etc
")
end

@testset "HTMLEnts" begin # see issue #206
    st = raw"""Blah π etc"""
    @test isapproxstr(st |> seval, "
Blah π etc
")
    # but ill-formed ones (either deliberately or not) will be parsed
    st = raw"""AT&T"""
    @test isapproxstr(st |> seval, "
AT&T
")
end

@testset "DoubleTicks" begin # see issue #204
    st = raw"""A `single` B"""
    steps = explore_md_steps(st)
    tokens = steps[:tokenization].tokens
    @test tokens[1].name == :CODE_SINGLE
    @test tokens[2].name == :CODE_SINGLE

    st = raw"""A ``double`` B"""
    steps = explore_md_steps(st)
    tokens = steps[:tokenization].tokens
    @test tokens[1].name == :CODE_DOUBLE
    @test tokens[2].name == :CODE_DOUBLE

    st = raw"""A `single` and ``double`` B"""
    steps = explore_md_steps(st)
    tokens = steps[:tokenization].tokens
    @test tokens[1].name == :CODE_SINGLE
    @test tokens[2].name == :CODE_SINGLE
    @test tokens[3].name == :CODE_DOUBLE
    @test tokens[4].name == :CODE_DOUBLE

    st = raw"""A `single` and ``double ` double`` B"""
    steps = explore_md_steps(st)
    blocks, tokens = steps[:ocblocks]
    @test blocks[2].name == :CODE_INLINE
    @test F.content(blocks[2]) == "double ` double"
    @test blocks[1].name == :CODE_INLINE
    @test F.content(blocks[1]) == "single"

    st = raw"""A `single` and ``double ` double`` and ``` triple ``` B"""
    steps = explore_md_steps(st)
    tokens = steps[:tokenization].tokens
    @test tokens[1].name == :CODE_SINGLE
    @test tokens[2].name == :CODE_SINGLE
    @test tokens[3].name == :CODE_DOUBLE
    @test tokens[4].name == :CODE_SINGLE
    @test tokens[5].name == :CODE_DOUBLE
    @test tokens[6].name == :CODE_TRIPLE
    @test tokens[7].name == :CODE_TRIPLE
    blocks, tokens = steps[:ocblocks]
    @test blocks[3].name == :CODE_BLOCK
    @test F.content(blocks[3]) == " triple "
    @test blocks[2].name == :CODE_INLINE
    @test blocks[1].name == :CODE_INLINE

    st = raw"""A `single` and ``double ` double`` and ``` triple `` triple```
               and ```julia 1+1``` and `single again` done"""
    steps = explore_md_steps(st)
    blocks, _ = steps[:ocblocks]
    @test blocks[4].name == :CODE_BLOCK_LANG
    @test F.content(blocks[4]) == " 1+1"
    @test blocks[3].name == :CODE_BLOCK
    @test F.content(blocks[3]) == " triple `` triple"
    @test blocks[2].name == :CODE_INLINE
    @test F.content(blocks[2]) == "double ` double"
    @test blocks[1].name == :CODE_INLINE
    @test F.content(blocks[1]) == "single"
end

@testset "\\ and \`" begin # see issue 203
    st = raw"""The `"Hello\n"` after the `readall` command is a returned value, whereas the `Hello` after the `run` command is printed output."""
    st |> seval
    @test isapproxstr(st |> seval, raw"""

The "Hello\n" after
                        the readall command is a returned value,
                        whereas the Hello after the run
                        command is printed output.
""")
end

@testset "i198" begin
    st = raw"""
            Essentially three things are imitated from LaTeX
            1. you can introduce definitions using `\newcommand`
            1. you can use hyper-references with `\eqref`, `\cite`, ...
            1. you can show nice maths (via KaTeX)
            """

    @test isapproxstr(st |> seval, raw"""

Essentially three things are imitated from LaTeX

you can introduce definitions using \newcommand

you can use hyper-references with \eqref, \cite, ...

you can show nice maths (via KaTeX)

                    """)
end

@testset "fixlinks" begin
   st = raw"""
        A [link] and
        B [link 2] and
        C [Python][] and
        D [a link][1] and
        blah
        [link]: https://julialang.org/
        [link 2]: https://www.mozilla.org/
        [Python]: https://www.python.org/
        [1]: http://slashdot.org/
        end
        """
    @test isapproxstr(st |> seval, """

                        A link and
                        B link 2 and
                        C Python and
                        D a link and blah end
""")
end

@testset "fixlinks2" begin
    st = raw"""
        A [link] and
        B ![link][id] and
        blah
        [link]: https://julialang.org/
        [id]: ./path/to/img.png
        """

    @test isapproxstr(st |> seval, """

                          A link and
                          B  and
                          blah

""")
end

@testset "fixlinks3" begin
    st = raw"""
        A [link] and
        B [unknown] and
        C ![link][id] and
        D
        [link]: https://julialang.org/
        [id]: ./path/to/img.png
        [not]: https://www.mozilla.org/
        """

    @test isapproxstr(st |> seval, """

                        A link and
                        B [unknown] and
                        C  and
                        D

""")
end

@testset "IndCode" begin # issue 207
    st = raw"""
        @def indented_code = true
        A

            a = 1+1
            if a > 1
                @show a
            end
            b = 2
            @show a+b
        end
        """
    @test isapproxstr(st |> seval, raw"""

A

a = 1+1
            if a > 1
                @show a
            end
            b = 2
            @show a+b

end
""")

    st = raw"""
        @def indented_code = true
        A `single` and ```python blah``` and

            a = 1+1
        then
        * blah
            + blih
            + bloh
        end
        """
    @test isapproxstr(st |> seval, raw"""

A single and

blah

and

a = 1+1

then

blah

blih

bloh

end

                        """)

    # st = raw"""
    #     @def indented_code = true
    #     A
    #
    #         function foo()
    #
    #             return 2
    #
    #         end
    #
    #         function bar()
    #             return 3
    #         end
    #
    #     B
    #
    #         function baz()
    #             return 5
    #
    #         end
    #
    #     C
    #     """
    # isapproxstr(st |> seval, raw"""
    #
A

function foo()
    #
    #                             return 2
    #
    #                         end
    #
    #                         function bar()
    #                             return 3
    #                         end

    #                         B
function baz()
    #                             return 5
    #
    #                         end

    #                         C

    #                         """)
end

@testset "More ``" begin
    st = raw"""
         A ``blah``.
         """
    @test isapproxstr(st |> seval, """
A blah.
""")
end

@testset "Issue 266" begin
    s = """Blah [`hello`] and later
       [`hello`]: https://github.com/cormullion/
       """ |> fd2html_td
    @test isapproxstr(s, """

Blah
        hello
        and later

        """)
end

@testset "TOML lang" begin
    s = raw"""
       blah
       ```TOML
       socrates
       ```
       end
       """
    @test isapproxstr(s |> fd2html_td, """

blah

socrates

end
""")
    s = raw"""
       blah
       ```julia-repl
       socrates
       ```
       end
       """
    @test isapproxstr(s |> fd2html_td, """

blah

socrates

end
""")
end


text:
#!/usr/bin/env julia

import Pkg
Pkg.update()
Pkg.add("Primes")
Pkg.add("Combinatorics")
# Pkg.add("LaTeXStrings")
Pkg.precompile()

max_stars_repo_name: robireton/math-capstone

max_issues_repo_name: robireton/math-capstone

max_forks_repo_name: robireton/math-capstone


text:
module NeuralNetworkAnalysis

include("init.jl")
include("problem.jl")
include("nnops.jl")
include("setops.jl")
include("split.jl")
include("utils.jl")
include("simulate.jl")
include("solve.jl")

# problem types
export ControlledPlant,
       BlackBoxController

# splitters
export BoxSplitter, ZonotopeSplitter,
       IndexedSplitter,
       SignSplitter

# solvers
export solve, forward, simulate,
       SampledApprox, VertexSolver, BoxSolver, SplitSolver, BlackBoxSolver

# utility functions
export @modelpath, read_nnet_mat, read_nnet_yaml, read_nnet_sherlock,
       write_nnet_sherlock,
       print_timed

end


text:
module Swizzles

export loop
export Nil, nil
export arrayify
export Guard
export Swizzle, Focus
export Pour, Beam
export Reduce, Sum
export Drop, SumOut
export @_, @$
export Delay, Intercept

include("base.jl")
include("util.jl")
include("properties.jl")

using Swizzles.Properties

include("ScalarArrays.jl")
include("NullArrays.jl")

include("WrapperArrays.jl")
include("ShallowArrays.jl")

include("StylishArrays.jl")
include("ArrayifiedArrays.jl")

include("Virtuals.jl")

include("ExtrudedArrays.jl")
include("NamedArrays.jl")
include("ValArrays.jl")

include("SwizzledArrays.jl")

include("EachindexArrays.jl")
include("PermutedArrays.jl")

struct Swizzle{Op, mask} <: Swizzles.Intercept
    op::Op
end

"""
    `Swizzle(op, mask...)`
    `Swizzle(op, mask)`

Create a function, `S`, which creates lazily reduced arrays using the operator
op. `S(A)` should produce an object which represents the reduction of `A` into a
result array `R`. `S(Z, A)` represents the reduction of `A` into a result array
`R` which has been initialized as `R .= Z`. Dimension `i` of `R` corresponds to
dimension `mask[i]` of `A` (if it exists). Dimensions of `A` which do not appear
in `mask` are reduced out. If `Z` is unspecified, the `initial` function is used
to create a suitable initial value. If no such initial value is found, the
initial value is `nothing` and `op` is wrapped in a `Guard`.

`Swizzles` can be materialized with `copy` or `copyto!`, and will eagerly fuse
themselves into special broadcast syntax. Broadcasting a `Swizzle` over a
broadcast expression `bc` will instead apply the `Swizzle` directy to `bc`
without materializing `bc` first. Thus, the code:
```
   y = 3.0 .+ Swizzle(+, 1).(x .* 2))
```
will result in one big happy fused operation:
```
   using Base.Broadcast: materialize, broadcasted
   y = materialize(broadcasted(+, 3.0, Swizzle(+, 1)(broadcasted(*, x, 2)))
```

See also: [`reduce`](@ref), [`Guard`](@ref)

# Examples
```jldoctest
julia> A = [1 2; 3 4; 5 6; 7 8; 9 10]
5×2 Array{Int64,2}:
 1   2
 3   4
 5   6
 7   8
 9  10
julia> Swizzle(+, 1).(A)
5-element Array{Int64,1}:
  3
  7
 11
 15
 19
julia> Swizzle(+).(A)
55
julia> Swizzle(+, nil, 2).(parse.(Int, ["1", "2"]))
1x2-element Array{Int64,1}:
 1 2
```
"""
@inline Swizzle(op::Op, _mask...) where {Op} = Swizzle{Op, _mask}(op)
@inline Swizzle(op::Op, _mask::Tuple) where {Op} = Swizzle{Op, _mask}(op)
@inline Swizzle(op::Op, ::Val{_mask}) where {Op, _mask} = Swizzle{Op, _mask}(op)

@inline function Properties.initial(ctr::Swizzle, arg)
    init = Properties.initial(ctr.op, eltype(arg))
    if init === nothing
        return nothing
    end
    return Ref(something(init))
end

@inline (ctr::Swizzle)(arg) = ctr(arrayify(arg))
@inline function (ctr::Swizzle{Op, _mask})(arg::AbstractArray) where {Op, _mask}
    init = initial(ctr, eltype_bound(arg))
    if init === nothing
        return Swizzle{Guard{Op}, _mask}(Guard(ctr.op))(Ref(nothing), arg)
    end
    return ctr(init, arg)
end
@inline function (ctr::Swizzle{Nothing, _mask})(arg::AbstractArray) where {_mask}
    return Swizzle{Nothing, _mask}(nothing)(Ref(nothing), arg)
end
@inline (ctr::Swizzle)(init, arg) = ctr(arrayify(init), arrayify(arg))
@inline function(ctr::Swizzle{Op, _mask})(init::Init, arg::Arg) where {Op, _mask, Arg <: AbstractArray, Init <: AbstractArray}
    if @generated
        mask = map(d -> d > ndims(arg) ? nil : d, _mask)
        return quote
            arr = SwizzledArray{Any, $(length(mask)), Op, $mask, Init, Arg}(ctr.op, init, arg)
            return convert(SwizzledArray{eltype_bound(arr)}, arr)
        end
    else
        mask = map(d -> d > ndims(arg) ? nil : d, _mask)
        arr = SwizzledArray{Any, length(mask), Op, mask, Init, Arg}(ctr.op, init, arg)
        return convert(SwizzledArray{eltype_bound(arr)}, arr)
    end
end

@inline mask(::Swizzle{<:Any, _mask}) where {_mask} = _mask
@inline mask(::Type{<:Swizzle{<:Any, _mask}}) where {_mask} = _mask

"""
    `Focus(mask...)`
    `Focus(mask)`

Similar to `Swizzle`, but the resulting operator `Y` creates lazy
transformations which do not reduce. Dimensions of the argument which do not
appear in `mask` are asserted to have length `1`.

See also: [`Swizzle`](@ref), [`permutedims`](@ref).

# Examples
```jldoctest
julia> A = [1 2 3 4 5]
1×5 Array{Int64,2}:
 1  2  3  4  5
julia> Focus(nil, nil, 2).(A)
1×1×5 Array{Int64,3}:
[:, :, 1] =
 1
[:, :, 2] =
 2
[:, :, 3] =
 3
[:, :, 4] =
 4
[:, :, 5] =
 5
```
"""
@inline Focus(_mask...) = Swizzle(nothing, _mask...)

struct Pour{Op, _imask} <: Swizzles.Intercept
    op::Op
end

"""
    `Pour(op, imask...)`
    `Pour(op, imask)`

Similar to [`Swizzle`](@ref), but the mask is "inverted". Creates an operator
which lazily reduces the input `A`. Dimension `i` of `A` is associated with
dimension `imask[i]` of the result. If `imask[i]` is `nil` (an instance of the
singleton type `Nil`) this dimension of `A` is reduced over. If
`i > length(imask)`, `imask[i]` is assumed to be `nil`.

See also: [`Swizzle`](@ref)

# Examples
```jldoctest
julia> A = [1 2; 3 4; 5 6; 7 8; 9 10]
5×2 Array{Int64,2}:
 1   2
 3   4
 5   6
 7   8
 9  10
julia> Pour(+, 2).(A)
1x5 Array{Int64,2}:
 3  7  11  15  19
julia> Pour(+, nil, 3).(A)
1×1×2 Array{Int64,3}:
[:, :, 1] =
 25
[:, :, 2] =
 30
```
"""
@inline Pour(op::Op, _imask...) where {Op} = Pour{Op, _imask}(op)
@inline Pour(op::Op, _imask::Tuple) where {Op} = Pour{Op, _imask}(op)
@inline Pour(op::Op, ::Val{_imask}) where {Op, _imask} = Pour{Op, _imask}(op)

@inline function parse_beam_mask(arr, _imask::Tuple{Vararg{Union{Int, Nil}}})
    N = max(0, _imask[1:min(ndims(arr), length(_imask))]...)
    return imasktuple(d->nil, identity, _imask, N)
end
@inline (ctr::Pour)(arg) = ctr(arrayify(arg))
@inline (ctr::Pour)(init, arg) = ctr(arrayify(init), arrayify(arg))
@inline function(ctr::Pour{Op, _imask})(arg::Arg) where {Op, _imask, Arg <: AbstractArray}
    if @generated
        mask = parse_beam_mask(Arg, _imask)
        return :(return Swizzle{Op, $mask}(ctr.op)(arg))
    else
        mask = parse_beam_mask(arg, _imask)
        return Swizzle(ctr.op, mask)(arg)
    end
end
@inline function(ctr::Pour{Op, _imask})(init::Init, arg::Arg) where {Op, _imask, Arg <: AbstractArray, Init <: AbstractArray}
    if @generated
        mask = parse_beam_mask(Arg, _imask)
        return :(return Swizzle{Op, $mask}(ctr.op)(init, arg))
    else
        mask = parse_beam_mask(arg, _imask)
        return Swizzle(ctr.op, mask)(init, arg)
    end
end

"""
    `Beam(imask...)`
    `Beam(imask)`

Similar to `Pour`, but the resulting operator `Y` creates lazy
transformations which do not reduce. Dimensions of the argument which do not
have integer destinations in `imask` are asserted to have length `1`.

See also: [`Pour`](@ref).

# Examples
```jldoctest
julia> A = [1 2 3 4 5]
1×5 Array{Int64,2}:
 1  2  3  4  5
julia> Beam(nil, 3).(A)
1×1×5 Array{Int64,3}:
[:, :, 1] =
 1
[:, :, 2] =
 2
[:, :, 3] =
 3
[:, :, 4] =
 4
[:, :, 5] =
 5
```
"""
@inline Beam(_imask...) = Pour(nothing, _imask...)

struct Drop{Op, dims} <: Swizzles.Intercept
    op::Op
end

"""
    `Drop(op, dims...)`
    `Drop(op, dims)`

Similar to `Swizzle`, but the resulting operator reduces only over the
dimensions listed in `dims`, dropping those dimensions. `Drop(op)` produces an
operator which reduces over all dimensions.

See also: [`Swizzle`](@ref), [`reduce`](@ref), [`dropdims`](@ref).

# Examples
```jldoctest
julia> A = [1 2; 3 4; 5 6; 7 8; 9 10]
5×2 Array{Int64,2}:
 1   2
 3   4
 5   6
 7   8
 9  10
julia> Drop(+, 2).(A)
3-element Array{Int64,1}:
 3
 7
 11
 15
 19
```
"""
@inline Drop(op::Op, dims...) where {Op} = Drop{Op, dims}(op)
@inline Drop(op) = Drop(op, Colon())
@inline Drop(op::Op, dims::Colon) where {Op} = Drop{Op, dims}(op)
@inline Drop(op::Op, dims::Tuple) where {Op} = Drop{Op, dims}(op)
@inline Drop(op::Op, ::Val{dims}) where {Op, dims} = Drop{Op, dims}(op)

@inline (ctr::Drop)(arg) = ctr(arrayify(arg))
@inline function parse_drop_mask(arr, dims::Tuple{Vararg{Int}})
    return (setdiff(1:ndims(arr), dims)...,)
end
@inline parse_drop_mask(arr, ::Colon) = ()
@inline function(ctr::Drop{Op, dims})(arg::Arg) where {Op, dims, Arg <: AbstractArray}
    if @generated
        mask = parse_drop_mask(Arg, dims)
        return :(return Swizzle{Op, $mask}(ctr.op)(arg))
    else
        mask = parse_drop_mask(arg, dims)
        return Swizzle(ctr.op, mask)(arg)
    end
end
@inline (ctr::Drop)(init, arg) = ctr(arrayify(init), arrayify(arg))
@inline function(ctr::Drop{Op, dims})(init::Init, arg::Arg) where {Op, dims, Arg <: AbstractArray, Init <: AbstractArray}
    if @generated
        mask = parse_drop_mask(Arg, dims)
        return :(return Swizzle{Op, $mask}(ctr.op)(init, arg))
    else
        mask = parse_drop_mask(arg, dims)
        return Swizzle(ctr.op, mask)(init, arg)
    end
end

"""
    `SumOut(dims...)`
    `SumOut(dims)`

Similar to [`Drop`](@ref), but `op` is set to `+`.

See also: [`Drop`](@ref).
"""
@inline SumOut(dims...) = Drop(Base.FastMath.add_fast, dims...)

struct Reduce{Op, dims} <: Swizzles.Intercept
    op::Op
end

"""
    `Reduce(op, dims...)`
    `Reduce(op, dims)`

Similar to `Swizzle`, but the resulting operator reduces only over the
dimensions listed in `dims`, reduceping those dimensions. `Reduce(op)` produces an
operator which reduces over all dimensions.

See also: [`Swizzle`](@ref), [`reduce`](@ref).

# Examples
```jldoctest
julia> A = [1 2; 3 4; 5 6; 7 8; 9 10]
5×2 Array{Int64,2}:
 1   2
 3   4
 5   6
 7   8
 9  10
julia> Reduce(+, 2).(A)
3-element Array{Int64,1}:
 3
 7
 11
 15
 19
```
"""
@inline Reduce(op::Op, dims...) where {Op} = Reduce{Op, dims}(op)
@inline Reduce(op) = Reduce(op, Colon())
@inline Reduce(op::Op, dims::Colon) where {Op} = Reduce{Op, dims}(op)
@inline Reduce(op::Op, dims::Tuple) where {Op} = Reduce{Op, dims}(op)
@inline Reduce(op::Op, ::Val{dims}) where {Op, dims} = Reduce{Op, dims}(op)

@inline (ctr::Reduce)(arg) = ctr(arrayify(arg))
@inline function parse_reduce_mask(arr, dims::Tuple{Vararg{Int}})
    return ntuple(d -> d in dims ? nil : d, Val(ndims(arr)))
end
@inline parse_reduce_mask(arr, ::Colon) = ()
@inline function(ctr::Reduce{Op, dims})(arg::Arg) where {Op, dims, Arg <: AbstractArray}
    if @generated
        mask = parse_reduce_mask(Arg, dims)
        return :(return Swizzle{Op, $mask}(ctr.op)(arg))
    else
        mask = parse_reduce_mask(arg, dims)
        return Swizzle(ctr.op, mask)(arg)
    end
end
@inline (ctr::Reduce)(init, arg) = ctr(arrayify(init), arrayify(arg))
@inline function(ctr::Reduce{Op, dims})(init::Init, arg::Arg) where {Op, dims, Arg <: AbstractArray, Init <: AbstractArray}
    if @generated
        mask = parse_reduce_mask(Arg, dims)
        return :(return Swizzle{Op, $mask}(ctr.op)(init, arg))
    else
        mask = parse_reduce_mask(arg, dims)
        return Swizzle(ctr.op, mask)(init, arg)
    end
end

"""
    `Sum(dims...)`
    `Sum(dims)`

Similar to [`Reduce`](@ref), but `op` is set to `+`.

See also: [`Reduce`](@ref).
"""
@inline Sum(dims...) = Reduce(Base.FastMath.add_fast, dims...)

include("sugar.jl")

end


text:
using MassCommits
using Test

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


text:
"""
    VectorFold{T,V <: AbstractVector{T}}
A mutable structure for folded vector that extends the methods of AbstractVector. Compared to `IVectorFold`, this tructure is about 20% faster using iterators.
Note that this structure keep an active pointer to the `current` *unfolded* pattern. However, its external behavior is similar to `IVectorFold`.
"""
mutable struct VectorFold{T,V<:AbstractVector{T}} <: AbstractVectorFold{T}
    pattern::V
    gap::T
    folds::Int
    current::Int
end

VectorFold(p, g, f; c=1) = VectorFold(p, g, f, c)

pattern_length(mvf::VectorFold) = length(mvf.pattern)

function pattern(mvf::VectorFold)
    distortion = gap(mvf) * (mvf.current - 1)
    return mvf.pattern .- distortion
end

pattern(mvf::VectorFold, index) = pattern(mvf)[index]

"""
    set_fold!(mvf::VectorFold, new_fold = mvf.current + 1)
Set the *unfolded* pattern to `new_fold`. By default move the next *fold* after `current`.
"""
function set_fold!(mvf, new_fold=mvf.current + 1)
    if new_fold != mvf.current && 0 < new_fold ≤ mvf.folds
        distortion = gap(mvf) * (new_fold - mvf.current)
        mvf.pattern .+= distortion
        mvf.current = new_fold
    end
end

# Base case iterate method
function Base.iterate(iter::VectorFold)
    reset_pattern!(iter)
    return pattern(iter, 1), 1
end

# "Induction" iterate method
function Base.iterate(iter::VectorFold, state::Int)
    state ≥ length(iter) && return nothing

    next_state = state + 1
    pl = pattern_length(iter)

    pattern_counter = mod1(next_state, pl)
    elem = iter.pattern[pattern_counter]
    pattern_counter == pl && set_fold!(iter)

    return elem, next_state
end

# Reverse iterate method
function Base.iterate(
    r_iter::Base.Iterators.Reverse{VectorFold{T,V}}, state::Int=length(r_iter.itr)
) where {T,V}
    state < 1 && return nothing

    iter = r_iter.itr
    state == length(iter) && set_fold!(iter, iter.folds)

    next_state = state - 1
    pl = pattern_length(iter)

    pattern_counter = mod(next_state, pl) + 1
    elem = iter.pattern[pattern_counter]
    pattern_counter == 1 && set_fold!(iter, iter.current - 1)

    return elem, next_state
end

# Specific dispatch for MVectorFold
function Base.rand(mvf::VectorFold)
    return Base.rand(mvf.pattern) +
           Base.rand((1 - mvf.current):(folds(mvf) - mvf.current)) * gap(mvf)
end

Base.rand(mvf::VectorFold, n::Int) = map(_ -> rand(mvf), 1:n)

# Specific dispatch for MVectorFold
function unfold(mvf::VectorFold; from=1, to=folds(mvf))
    pl = pattern_length(mvf)
    ul = (to - from + 1) * pl
    v = typeof(mvf.pattern)(undef, ul)

    count = 0
    for fold in from:to
        set_fold!(mvf, fold)
        for i in 1:pl
            v[count * pl + i] = mvf.pattern[i]
        end
        count += 1
    end
    return v
end

make_vector_fold(pattern, gap, fold, ::Val{:mutable}) = VectorFold(pattern, gap, fold)


text:
using JSON
using JuMP
using NLopt
#using Winston
using ThermalTraj
using PyCall
@pyimport helper

OBJSHIFT = [0.0,0.0]
OBJ = [0.0,0.0]

smclient = connect(ip"127.0.0.1",5512)
println("Connected to StateMachine server")

# server on 5512 is for state machine connection
@async begin
	while true
		println("waiting for new message")
		smmsg = readline(smclient)
		smdata = JSON.parse(smmsg)

		OBJSHIFT[1] = -smdata["SOBJX"]
		OBJSHIFT[2] = -smdata["SOBJY"]
		println("SET SHFTABS ", OBJSHIFT[1], ", ", OBJSHIFT[2])
	end
end

#create connection to Python Autopilot
client = connect(ip"127.0.0.1",5510)
println("Connected")

# and immediately request the current position from Python
# Commands to send to Python Autopilot
channel = "REQUEST_DATA"

my_dict = ["CHANNEL"=> channel]
reqPosMsg = json(my_dict)

# placeholder for trajectory tape
nodes = 35
x = zeros(5,nodes);
u = zeros(nodes);

Re = 6378137

dt = 0.5;
t = linspace(0,dt*nodes,nodes);

# run behaviors
tapeMsg = ""

#request state from python
println("Requesting data from Python controller")
print(client,reqPosMsg)
#block on reply from python
msgdata = readline(client)

### known position of thermal only using for initial testing, shiftabs used to move objective center
# Thermal at KSO tower
#<latitude>37.61633</latitude>
#<longitude>-122.38334</longitude>
#<strength-fps>13.49</strength-fps>
#<diameter-ft>4800.0</diameter-ft>
#Tlat = 37.61633
#Tlon = -122.38334

msg = JSON.parse(msgdata)
println("Now prepare for first trajectory optimization")
#dlat = (msg["LAT"]*180.0/pi - Tlat)*pi/180.0
#dlon = (msg["LON"]*180.0/pi - Tlon)*pi/180.0

x0 = zeros(5)
println("Heading received ", msg["HEADING"])
x0[3] = heading2XY(msg["HEADING"]*180.0/pi)
#x0[1] = Re*dlon
#x0[2] = Re*dlat
dx, dy = helper.lla2flatearth(msg["LAT"], msg["LON"])
x0[1] = dx
x0[2] = dy

if (x0[1]>100000.0)
	println("ERROR, x0 is way to large for trajectory optimization")
end

lastSendT = 0.0
nextSendDelay = 0.0

while true
	println("x0: ", x0)

	#solve new opt problem
	solveTraj!(x0,u,x)

	#pack new json tape
	#x[4,:] = (sign(u).*max(abs(x[4,:]'),abs(u)))'
	tapeX = x[1,:] - OBJ[1]
	tapeY = x[2,:] - OBJ[2]
	#OBJ[1] = 0.0
	#OBJ[2] = 0.0
	tape = ["CHANNEL"=> "NEW_TAPE", "t"=> t, "u"=>-x[4,:], "x"=>tapeX, "y"=>tapeY]
	tapeMsg = json(tape)
	println("Ready with new tape")

	if lastSendT != 0.0
		for j in 1:200
			#println("Requesting new data")
			print(client,reqPosMsg)
			msgdata = readline(client)
			chkNewTapeMsg = JSON.parse(msgdata)
			if chkNewTapeMsg["CUR_INDEX"] >= (pp-1)
				println("PP ", pp, " | ", nodes)
				break
			else
				sleep( dt/2.0 )
			end

			if j>= 199
				println("Failed to reach the PP new tape goal in a reasonable amount of time, sending new tape.")
			end
		end
	end

	#send tape
	println("Sending new tape")
	print(client,tapeMsg)
	lastSendT = time()

	pp = floor(83.0*nodes/100.0)
	#nextSendDelay = t[pp]
	rpx0 = x[:,pp]
	#println(rpx0)
	if rpx0[3]<-pi
		rpx0[3] += 2pi
	end
	if rpx0[3]>pi
		rpx0[3] -= 2pi
	end
	#println(rpx0)
	x0 = rpx0
	# incorporate shift from state machine logic
	OBJ[1] = OBJ[1] + OBJSHIFT[1]
	OBJ[2] = OBJ[2] + OBJSHIFT[2]
	println("SHFTABS ", -OBJ[1], ", ", -OBJ[2])
	x0[1] = x0[1] + OBJSHIFT[1]
	x0[2] = x0[2] + OBJSHIFT[2]
	OBJSHIFT[1] = 0.0
	OBJSHIFT[2] = 0.0
	# reset shift after it has been incorporated into the current plan
end

close(client)
close(smclient)

max_stars_repo_name: dehann/ThermalSoaring

max_issues_repo_name: dehann/ThermalSoaring

max_forks_repo_name: dehann/ThermalSoaring


text:
"""
    load_station_metadata() -> DataFrame

Loads the station metadata detailed in section 4 of [1], returning it in a `DataFrame` with
columns named according the the `Variable` column below:

```
------------------------------
Variable   Columns   Type
------------------------------
ID            1-11   Character
LATITUDE     13-20   Real
LONGITUDE    22-30   Real
ELEVATION    32-37   Real
STATE        39-40   Character
NAME         42-71   Character
GSN FLAG     73-75   Character
HCN/CRN FLAG 77-79   Character
WMO ID       81-85   Character
------------------------------
```

[1] - https://www1.ncdc.noaa.gov/pub/data/ghcn/daily/readme.txt
"""
function load_station_metadata()
    raw_data = CSV.read(
        datadep"ghcn-stations" * "/" * station_metadata_fname, DataFrame;
        header=0, delim='\t',
    ).Column1
    N_stations = length(raw_data)

    parsed_data = (
        ID = Vector{String}(undef, N_stations),
        LATITUDE = Vector{Float32}(undef, N_stations),
        LONGITUDE = Vector{Float32}(undef, N_stations),
        ELEVATION = Vector{Float32}(undef, N_stations),
        STATE = Vector{String}(undef, N_stations),
        NAME = Vector{String}(undef, N_stations),
        GSN_FLAG = Vector{String}(undef, N_stations),
        HCN_CRN_FLAG = Vector{String}(undef, N_stations),
        WMO_ID = Vector{String}(undef, N_stations),
    )

    for (n, row) in enumerate(raw_data)
        parsed_data.ID[n] = row[1:11]
        parsed_data.LATITUDE[n] = parse(Float32, row[13:20])
        parsed_data.LONGITUDE[n] = parse(Float32, row[22:30])
        parsed_data.ELEVATION[n] = parse(Float32, row[32:37])
        parsed_data.STATE[n] = strip(row[39:40])
        parsed_data.NAME[n] = strip(row[42:71])
        parsed_data.GSN_FLAG[n] = strip(row[73:75])
        parsed_data.HCN_CRN_FLAG[n] = strip(row[77:79])
        parsed_data.WMO_ID[n] = strip(row[81:85])
    end

    return DataFrame(parsed_data)
end

"""
    load_countries_metadata() -> DataFrame

Loads the countries metadata detailed in section 5 of [1], returning the result as a
`DataFrame` with columns `CODE` and `NAME`.

[1] - https://www1.ncdc.noaa.gov/pub/data/ghcn/daily/readme.txt
"""
function load_countries_metadata()
    raw_data = CSV.read(
        datadep"ghcn-countries" * "/" * countries_fname, DataFrame;
        header=0, delim='\t',
    ).Column1
    N_countries = length(raw_data)

    parsed_data = (
        CODE = Vector{String}(undef, N_countries),
        NAME = Vector{String}(undef, N_countries),
    )

    for (n, row) in enumerate(raw_data)
        parsed_data.CODE[n] = row[1:2]
        parsed_data.NAME[n] = strip(row[4:end])
    end
    return DataFrame(parsed_data)
end

"""
    load_inventories() -> DataFrame

Loads the inventory file detailed in section 7 of [1], and returns a `DataFrame` with
columns named according to the variables listed in the file-format information below.

```
------------------------------
Variable   Columns   Type
------------------------------
ID            1-11   Character
LATITUDE     13-20   Real
LONGITUDE    22-30   Real
ELEMENT      32-35   Character
FIRSTYEAR    37-40   Integer
LASTYEAR     42-45   Integer
------------------------------
```

[1] - https://www1.ncdc.noaa.gov/pub/data/ghcn/daily/readme.txt
"""
function load_inventories()
    raw_data = CSV.read(
        datadep"ghcn-inventories" * "/" * inventory_fname, DataFrame;
        header=0, delim='\t',
    ).Column1
    N_stations = length(raw_data)

    parsed_data = (
        ID = Vector{String}(undef, N_stations),
        LATITUDE = Vector{Float32}(undef, N_stations),
        LONGITUDE = Vector{Float32}(undef, N_stations),
        ELEMENT = Vector{String}(undef, N_stations),
        FIRSTYEAR = Vector{Int}(undef, N_stations),
        LASTYEAR = Vector{Int}(undef, N_stations),
    )

    for (n, row) in enumerate(raw_data)
        parsed_data.ID[n] = row[1:11]
        parsed_data.LATITUDE[n] = parse(Float32, row[13:20])
        parsed_data.LONGITUDE[n] = parse(Float32, row[22:30])
        parsed_data.ELEMENT[n] = row[32:35]
        parsed_data.FIRSTYEAR[n] = parse(Int, row[37:40])
        parsed_data.LASTYEAR[n] = parse(Int, row[42:45])
    end

    return DataFrame(parsed_data)
end

"""
    load_data_file(station_id::String) -> DataFrame

See section 3 of [1].

```
------------------------------
Variable   Columns   Type
------------------------------
ID            1-11   Character
YEAR         12-15   Integer
MONTH        16-17   Integer
ELEMENT      18-21   Character
VALUE1       22-26   Integer
MFLAG1       27-27   Character
QFLAG1       28-28   Character
SFLAG1       29-29   Character
VALUE2       30-34   Integer
MFLAG2       35-35   Character
QFLAG2       36-36   Character
SFLAG2       37-37   Character
  .           .          .
  .           .          .
  .           .          .
VALUE31    262-266   Integer
MFLAG31    267-267   Character
QFLAG31    268-268   Character
SFLAG31    269-269   Character
------------------------------
```

[1] - https://www1.ncdc.noaa.gov/pub/data/ghcn/daily/readme.txt
"""
function load_data_file(station_id::String)
    fdir = datadep"ghcn-data" * "/" * ghcnd_data_dir
    file_name = joinpath(fdir, station_id * ".dly")
    raw_data = CSV.read(
        file_name, DataFrame;
        header=0, delim='\t',
    ).Column1
    N_months = length(raw_data)

    # Construct a `NamedTuple` to store the parsed data.
    column_names = vcat(
        [:ID, :YEAR, :MONTH, :ELEMENT],
        map(1:31) do n
            [Symbol("VALUE$n"), Symbol("MFLAG$n"), Symbol("QFLAG$n"), Symbol("SFLAG$n")]
        end...,
    )
    columns = vcat(
        [
            Vector{String}(undef, N_months),
            Vector{Int}(undef, N_months),
            Vector{Int}(undef, N_months),
            Vector{String}(undef, N_months),
        ],
        map(1:31) do n
            [
                Vector{Union{Int, Missing}}(undef, N_months),
                Vector{String}(undef, N_months),
                Vector{String}(undef, N_months),
                Vector{String}(undef, N_months),
            ]
        end...,
    )

    # Specify the beginning and end location of each column of the data.
    column_locations = vcat(
        [(1, 11), (12, 15), (16, 17), (18, 21)],
        map(1:31) do n
            start_pos = 22 + 8 * (n - 1)
            [
                (start_pos, start_pos + 4),
                (start_pos + 5, start_pos + 5),
                (start_pos + 6, start_pos + 6),
                (start_pos + 7, start_pos + 7),
            ]
        end...,
    )

    # Iterate over the data.
    for (n, row) in enumerate(raw_data)
        for col in eachindex(columns)
            start_pos = first(column_locations[col])
            end_pos = last(column_locations[col])
            columns[col][n] = _to(eltype(columns[col]), row[start_pos:end_pos])
        end
    end

    return DataFrame(columns, column_names)
end

function _to(::Type{Union{Missing, T}}, s::String) where {T<:Number}
    if s == "-9999"
        return missing
    else
        return parse(T, s)
    end
end
_to(::Type{String}, s::String) = s
_to(::Type{Int}, s::String) = parse(Int, s)


text:
## header to provide surface and context
using Cairo
c = CairoRGBSurface(256,256);
cr = CairoContext(c);

save(cr);
set_source_rgb(cr,0.8,0.8,0.8);    # light gray
rectangle(cr,0.0,0.0,256.0,256.0); # background
fill(cr);
restore(cr);

save(cr);
## original example, taken from cairo/test/pattern-mesh.c
function draw(cr, width::Int, height::Int)

	PAT_WIDTH  = 170
	PAT_HEIGHT = 170
	SIZE = PAT_WIDTH
	PAD = 2
	WIDTH = (PAD + SIZE + PAD)
	HEIGHT = WIDTH

	pattern = CairoPatternMesh();

	#cairo_test_paint_checkered (cr); ?

	translate(cr, PAD, PAD);
	translate(cr, 10, 10);

	mesh_pattern_begin_patch(pattern);

	mesh_pattern_move_to(pattern, 0, 0);
	mesh_pattern_curve_to(pattern, 30, -30,  60,  30, 100, 0);
	mesh_pattern_curve_to(pattern, 60,  30, 130,  60, 100, 100);
	mesh_pattern_curve_to(pattern, 60,  70,  30, 130,   0, 100);
	mesh_pattern_curve_to(pattern, 30,  70, -30,  30,   0, 0);

	mesh_pattern_set_corner_color_rgb(pattern, 0, 1, 0, 0);
	mesh_pattern_set_corner_color_rgb(pattern, 1, 0, 1, 0);
	mesh_pattern_set_corner_color_rgb(pattern, 2, 0, 0, 1);
	mesh_pattern_set_corner_color_rgb(pattern, 3, 1, 1, 0);

	mesh_pattern_end_patch(pattern);

	mesh_pattern_begin_patch(pattern);

	mesh_pattern_move_to(pattern, 50, 50);
	mesh_pattern_curve_to(pattern,  80,  20, 110,  80, 150, 50);
	mesh_pattern_curve_to(pattern, 110,  80, 180, 110, 150, 150);
	mesh_pattern_curve_to(pattern, 110, 120,  80, 180,  50, 150);
	mesh_pattern_curve_to(pattern,  80, 120,  20,  80,  50, 50);

	mesh_pattern_set_corner_color_rgba(pattern, 0, 1, 0, 0, 0.3);
	mesh_pattern_set_corner_color_rgb(pattern, 1, 0, 1, 0);
	mesh_pattern_set_corner_color_rgba(pattern, 2, 0, 0, 1, 0.3);
	mesh_pattern_set_corner_color_rgb(pattern, 3, 1, 1, 0);

	mesh_pattern_end_patch(pattern);

	set_source(cr, pattern);
	paint(cr);
	#pattern_destroy (pattern);

end

draw(cr,170,170);

## mark picture with current date
restore(cr);
move_to(cr,0.0,12.0);
set_source_rgb(cr, 0,0,0);
show_text(cr,Libc.strftime(time()));
write_to_png(c,"sample_meshpattern.png");


text:
struct ProcessIterativeInitialization{T <: ProcessingMethod} <: ProcessingMethod
    stop_criteria::StopCriteria
    internal_process_method::T
    n_steps::Int
end
function ProcessIterativeInitialization(ϵ, problem, process_method)
    nx, ny = problem.NX, problem.NY

    ρ = zeros(nx, ny)
    ρ_old = zeros(nx, ny)

    ProcessIterativeInitialization(DensityConvergence(ϵ, ρ, ρ_old), process_method, 100)
end
function next!(process_method::ProcessIterativeInitialization{T}, q, f_in, t) where {T}
    if mod(t, max(10, round(Int, process_method.n_steps / 25))) == 0
        # next!(process_method.internal_process_method, q, f_in, 0)
    end

    if (should_stop!(process_method.stop_criteria, q, f_in))
        return true
    end

    return false
end


text:
function setup_rhs(tmp::EvolVars, bulkconstrains::BulkPartition{Nsys},
                   bulkderivs::BulkPartition{Nsys},
                   cache::HorizonCache, systems::SystemPartition,
                   integration::Integration) where {Nsys}
    # function to solve the nested system
    nested = Nested(systems, bulkconstrains, bulkderivs)

    diss_cache = similar(tmp)

    bulkevols_cache   = getbulkevolvedpartition(diss_cache)
    boundary_cache    = getboundary(diss_cache)
    gauge_cache       = getgauge(diss_cache)

    function (ff_t::EvolVars, ff::EvolVars, evoleq::EvolutionEquations, t)
        bulkevols_t = getbulkevolvedpartition(ff_t)
        boundary_t  = getboundary(ff_t)
        gauge_t     = getgauge(ff_t)

        bulkevols   = getbulkevolvedpartition(ff)
        boundary    = getboundary(ff)
        gauge       = getgauge(ff)

        # filter state vector after each integration (sub)step.
        vprint("INFO: filtering")

        # Note: one problem with this approach is that it always filters the
        # state vector upon calling, ie, if this function is called twice with
        # the *same* state vector (say, for error estimate purposes), the state
        # vector will be filtered *twice*.
        if t > 0 && integration.filter_poststep
            @inbounds for aa in 1:Nsys
                sys            = systems[aa]
                bulkevol       = bulkevols[aa]
                bulkevol_cache = bulkevols_cache[aa]

                apply_dissipation!(bulkevol, bulkevol_cache, sys)

                # exponential filter
                sys.filters(bulkevols[aa])
            end
            apply_dissipation!(boundary, boundary_cache, systems[1])
            apply_dissipation!(gauge, gauge_cache,  systems[Nsys])
        end

        vprint("INFO: compute_boundary_t")
        compute_boundary_t!(boundary_t, bulkevols[1], boundary, gauge, systems[1], evoleq)

        # solve nested system for the constrained variables
        vprint("INFO: nested system")
        nested(bulkevols, boundary, gauge, evoleq)

        vprint("INFO: compute_xi_t")
        compute_xi_t!(gauge_t, bulkconstrains[Nsys], bulkevols[Nsys], bulkderivs[Nsys],
                      gauge, cache, systems[Nsys], evoleq.gaugecondition)

        vprint("INFO: bulkevolved_t")
        # TODO: check if this loop is thread-safe
        # @inbounds @threads for aa in 1:Nsys
        @inbounds for aa in 1:Nsys
            sys           = systems[aa]
            bulkevol_t    = bulkevols_t[aa]
            bulkevol      = bulkevols[aa]
            bulkconstrain = bulkconstrains[aa]

            compute_bulkevolved_t!(bulkevol_t, bulkconstrain, gauge_t, bulkevol,
                                   boundary, gauge, sys, evoleq)
        end
        sync_bulkevolved!(bulkevols_t, bulkconstrains, gauge_t, systems, evoleq)

        nothing
    end
end


text:
include(joinpath(pwd(), "models/visualize.jl"))

function visualize!(vis, model,
    q_body, q_f1, q_f2, q_f3, q_f4;
    q_body_ref = q_body,
    ref_tl = 0.5,
    r = 0.025,
    l_torso = 0.367,
    w_torso = 0.267,
    Δt = 0.1)

	default_background!(vis)

	setobject!(vis["torso"],
    	Rect(Vec(-0.5 * l_torso, -0.5 * w_torso, -0.025),Vec(l_torso, w_torso, 0.05)),
    	MeshPhongMaterial(color = RGBA(0.0, 0.0, 0.0, 1.0)))

    setobject!(vis["torso_ref"],
    	Rect(Vec(-0.5 * l_torso, -0.5 * w_torso, -0.025),Vec(l_torso, w_torso, 0.05)),
    	MeshPhongMaterial(color = RGBA(0.0, 0.0, 0.0, ref_tl)))

	feet1 = setobject!(vis["feet1"], Sphere(Point3f0(0),
        convert(Float32, r)),
        MeshPhongMaterial(color = RGBA(1.0, 165.0 / 255.0, 0, 1.0)))

	feet2 = setobject!(vis["feet2"], Sphere(Point3f0(0),
		convert(Float32, r)),
		MeshPhongMaterial(color = RGBA(1.0, 165.0 / 255.0, 0, 1.0)))

	feet3 = setobject!(vis["feet3"], Sphere(Point3f0(0),
		convert(Float32, r)),
		MeshPhongMaterial(color = RGBA(1.0, 165.0 / 255.0, 0, 1.0)))

	feet4 = setobject!(vis["feet4"], Sphere(Point3f0(0),
        convert(Float32, r)),
        MeshPhongMaterial(color = RGBA(1.0, 165.0 / 255.0, 0, 1.0)))

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

	T = length(q_body)
	p_shift = [0.0, 0.0, r]

	for t = 1:T
		MeshCat.atframe(anim, t) do
			rot = MRP(q_body[t][4:6]...)
            rot_ref = MRP(q_body_ref[t][4:6]...)

			p_torso = q_body[t][1:3] + p_shift
            p_torso_ref = q_body_ref[t][1:3] + p_shift
			p_foot1 = q_f1[t] + p_shift
			p_foot2 = q_f2[t] + p_shift
			p_foot3 = q_f3[t] + p_shift
			p_foot4 = q_f4[t] + p_shift

			settransform!(vis["torso"], compose(Translation(p_torso), LinearMap(rot)))
            settransform!(vis["torso_ref"], compose(Translation(p_torso_ref), LinearMap(rot_ref)))
			settransform!(vis["feet1"], Translation(p_foot1))
			settransform!(vis["feet2"], Translation(p_foot2))
			settransform!(vis["feet3"], Translation(p_foot3))
			settransform!(vis["feet4"], Translation(p_foot4))
		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 foot_vis!(vis, i, t, contact_modes;
    h = 0.005,
    r = 0.035,
    tl = 0.5)

    if convert(Bool, contact_modes[t][i])
        setobject!(vis["f$(i)_$(t)"],
            Cylinder(Point(0.0, 0.0, 0.0), Point(0.0, 0.0, h), r),
            MeshPhongMaterial(color = RGBA(0.0, 1.0, 0.0, tl)))
        settransform!(vis["f$(i)_$(t)"], Translation(eval(Symbol("pf$(i)_ref"))[t]))
    end
    # setvisible!(vis["f$(i)_$(j)_$(t)"], convert(Bool, contact_modes[t][i]))
end

function cone_vis!(vis, i, t, contact_modes;
    h = 0.05,
    l = h * 1.0,
    w = h * 1.45,
    tl = 0.5,
    n = 10)

    pyramid = Pyramid(Point3(0.0, 0.0, 0.0), l, w)

    if convert(Bool, contact_modes[t][i])
        for p = 1:n
            setobject!(vis["pyramid$(i)_$(p)_$(t)"], pyramid,
                MeshPhongMaterial(
                # color = RGBA(1,153/255,51/255, 1.0))
                color = RGBA(51/255,1,1, tl))
                )
            settransform!(vis["pyramid$(i)_$(p)_$(t)"],
                compose(Translation(eval(Symbol("pf$(i)_ref"))[t][1], eval(Symbol("pf$(i)_ref"))[t][2], l),
                    LinearMap(RotX(π) * RotZ(π * p / n))),
                    )
            setvisible!(vis["pyramid$(i)_$(p)_$(t)"], convert(Bool, contact_modes[t][i]))
        end
    end
end

function force_vis!(vis, u, i, t;
    μ_friction = 1.0,
    f_max = 1.0,
    cone_vis_scale = 0.1)

    u_idx = (i - 1) * 3 .+ (1:3)
    ui = u[t][u_idx]
    us = cone_vis_scale * ui / f_max
    us[3] *= μ_friction
    p1 = eval(Symbol("pf$(i)_ref"))[t]

    force_vis = ArrowVisualizer(vis["force_$(i)_$(t)"])
    setobject!(force_vis, MeshPhongMaterial(color=RGBA(1.0, 0.0, 0.0, 1.0)))

    settransform!(force_vis,
    			Point(p1[1] - 0.0 * us[1], p1[2] - 0.0 * us[2], p1[3] - 0.0 * us[3]),
    			Vec(us[1], us[2], us[3]),
    			shaft_radius=0.01,
    			max_head_radius=0.025)
end

max_issues_repo_name: thowell/DirectMotionPlanning


text:
mutable struct Computer
    name::String
    board::DynamicalModel
    software::Vector{DynamicalModel}
    Computer(name = "computer",
             board = DynamicalModel("null_board"),
             software = Vector{DynamicalModel}()) = new(name, board, software)
end


text:
push!(LOAD_PATH, "./src")

# "Packages"
using Revise # reduce need for recompile
using Plots
using Documenter
using LinearAlgebra
using SparseArrays

# "User defined modules"
using CommonUtils
using Basis1D
using Basis2DQuad
using Basis3DHex
using UniformHexMesh

using SetupDG
using UnPack

N = 3
K1D = 8
CFL = .75;
T = .5 # endtime

VX,VY,VZ,EToV = uniform_hex_mesh(K1D,K1D,K1D)

rd = init_reference_hex(N,gauss_lobatto_quad(0,0,N))
md = init_mesh((VX,VY,VZ),EToV,rd)

# make domain periodic
@unpack Nfaces = rd
@unpack xf,yf,zf,K,mapM,mapP,mapB = md
LX = 2; LY = 2; LZ = 2
mapPB = build_periodic_boundary_maps(xf,yf,zf,LX,LY,LZ,Nfaces*K,mapM,mapP,mapB)
mapP[mapB] = mapPB
@pack! md = mapP

"Time integration"
rk4a,rk4b,rk4c = rk45_coeffs()
CN = (N+1)*(N+2)*3/2  # estimated trace constant
dt = CFL * 2 / (CN*K1D)
Nsteps = convert(Int,ceil(T/dt))
dt = T/Nsteps

function rhs(u,rd::RefElemData,md::MeshData,params...)

    @unpack Dr,Ds,Dt,LIFT,Vf = rd
    @unpack rxJ,sxJ,txJ,ryJ,syJ,tyJ,rzJ,szJ,tzJ,J,nxJ,nyJ,nzJ,sJ = md
    @unpack mapP,mapB = md

    uM = Vf*u
    du = uM[mapP]-uM
    uflux = @. .5*(du*nxJ - abs(nxJ)*du)

    ur,us,ut = (A->A*u).((Dr,Ds,Dt))
    ux = @. rxJ*ur + sxJ*us + txJ*ut
    rhsu = ux + LIFT*uflux

    return @. -rhsu/J
end

# set initial conditions
@unpack x,y,z = md
u = @. sin(pi*x)
u = @. exp(-25*(x^2+y^2))
resQ = zeros(size(x))
for i = 1:Nsteps
    for INTRK = 1:5
        rhsQ = rhs(u,rd,md)
        @. resQ = rk4a[INTRK]*resQ + dt*rhsQ
        @. u    = u + rk4b[INTRK]*resQ
    end

    if i%10==0 || i==Nsteps
        println("Number of time steps: ", i, " out of ", Nsteps)
    end
end

"plotting nodes"
gr(aspect_ratio=1,legend=false,
   markerstrokewidth=0,markersize=2)

@unpack Vp = rd
xp,yp,zp = (x->Vp*x).((x,y,z))
vv = Vp*u
ids = map(x->x[1],findall(@. abs(zp[:])<1e-10))
(xp,yp,zp,vv) = (x->x[ids]).((xp,yp,zp,vv))
scatter(xp,yp,vv,zcolor=vv,camera=(0,90))

max_stars_repo_name: yiminllin/ESDG-Fourier

max_issues_repo_name: yiminllin/ESDG-Fourier

max_forks_repo_name: yiminllin/ESDG-Fourier


text:
@testset "Tagging_Schemes" begin
    @testset "BIO1 and BIO2" begin
        tags_BIO1 = ["I-LOC", "O", "I-PER", "B-MISC", "I-MISC", "I-ORG"]
        tags_BIO2 = ["B-LOC", "O", "B-PER", "B-MISC", "I-MISC", "B-ORG"]

        output_tags = deepcopy(tags_BIO1)
        tag_scheme!(tags_BIO1, "BIO1", "BIO2")
        @test tags_BIO1 == tags_BIO2

        tag_scheme!(tags_BIO1, "BIO2", "BIO1")
        @test tags_BIO1 == output_tags
    end

    @testset "BIO1 and BIOES" begin
        tags_BIO1 = ["I-LOC", "O", "I-PER", "B-MISC", "I-MISC", "B-PER",
                                                        "I-PER", "I-PER"]
        tags_BIOES = ["S-LOC", "O", "S-PER", "B-MISC", "E-MISC", "B-PER",
                                                        "I-PER", "E-PER"]

        output_tags = deepcopy(tags_BIO1)
        tag_scheme!(tags_BIO1, "BIO1", "BIOES")
        @test tags_BIO1 == tags_BIOES

        tag_scheme!(tags_BIO1, "BIOES", "BIO1")
        @test tags_BIO1 == output_tags
    end

    @testset "BIO2 and BIOES" begin
        tags_BIO2 = ["B-LOC", "O", "B-PER", "B-MISC", "I-MISC", "B-PER",
                                                        "I-PER", "I-PER"]
        tags_BIOES = ["S-LOC", "O", "S-PER", "B-MISC", "E-MISC", "B-PER",
                                                        "I-PER", "E-PER"]

        output_tags = deepcopy(tags_BIO2)
        tag_scheme!(tags_BIO2, "BIO2", "BIOES")
        @test tags_BIO2 == tags_BIOES

        tag_scheme!(tags_BIO2, "BIOES", "BIO2")
        @test tags_BIO2 == output_tags
    end
end


text:
if VERSION < v"1.0.0"
    error("VS Code julia language server only works with julia 1.0.0+")
end

@info "Starting the Julia Language Server"

using InteractiveUtils, Sockets

include("../error_handler.jl")

struct LSPrecompileFailure <: Exception
    msg::AbstractString
end

function Base.showerror(io::IO, ex::LSPrecompileFailure)
    print(io, ex.msg)
end

try
    if length(Base.ARGS) != 5
        error("Invalid number of arguments passed to julia language server.")
    end

    conn = stdout
    (outRead, outWrite) = redirect_stdout()

    if Base.ARGS[2] == "--debug=yes"
        ENV["JULIA_DEBUG"] = "all"
    elseif Base.ARGS[2] != "--debug=no"
        error("Invalid argument passed.")
    end

    try
        using LanguageServer, SymbolServer
    catch err
        if err isa ErrorException && startswith(err.msg, "Failed to precompile")
            throw(LSPrecompileFailure(err.msg))
        else
            rethrow(err)
        end
    end

    symserver_store_path = joinpath(ARGS[5], "symbolstorev4")

    if !ispath(symserver_store_path)
        mkpath(symserver_store_path)
    end

    @info "Symbol server store is at '$symserver_store_path'."

    server = LanguageServerInstance(
        stdin,
        conn,
        Base.ARGS[1],
        Base.ARGS[4],
        (err, bt)->global_err_handler(err, bt, Base.ARGS[3], "Language Server"),
        symserver_store_path
    )
    run(server)
catch err
    global_err_handler(err, catch_backtrace(), Base.ARGS[3], "Language Server")
end

max_stars_repo_name: pgunnink/julia-vscode

max_issues_repo_name: pgunnink/julia-vscode

max_forks_repo_name: pgunnink/julia-vscode


text:
"""
Adaptively preconditioned optimisation algorithm following
J. Duchi, E. Hazan, Y. Singer; Adaptive Subgradient Methods for Online Learning and Stochastic Optimization, 12(Jul):2121−2159, 2011.

Fields:
    -x: position
    -stepcount: counts the steps
    -mass: runing componentwise average of grad.^2
    -stepsize: stepsize
    -minmass: prevent division by zero and limit preconditioner

Constructors:
    - construct using `AdagradState(x::Array{Float64,1};stepsize=.001,mass=ones(size(x)),stepcount=0,minmass=1e-10)`
"""

type AdagradState <: SamplerState
  x::Array{Float64,1}
  mass::Array{Float64,1}
  stepcount::Int32

  stepsize::Float64
  minmass::Float64

  function AdagradState(x::Array{Float64,1};
    stepsize=.001,
    mass=ones(size(x)),
    stepcount=0,
    minmass=1e-10)

    new(x,mass*stepcount,stepcount,stepsize,minmass)
  end
end

function sample!(s::AdagradState, grad::Function)
  g = grad(s.x)
  s.stepcount += 1
  if s.stepcount == 1
    s.mass = g.*g
  else
    s.mass += g.*g
  end
    mm = sqrt(max(s.minmass,s.mass))

  s.x[:] += s.stepsize./mm .* g
end


text:
using Plots, Test
pgfplotsx()

function create_plot( args...; kwargs... )
    pgfx_plot = plot(args...; kwargs...)
    return pgfx_plot, repr("application/x-tex", pgfx_plot)
end

function create_plot!( args...; kwargs... )
    pgfx_plot = plot!(args...; kwargs...)
    return pgfx_plot, repr("application/x-tex", pgfx_plot)
end

@testset "PGFPlotsX" begin
    pgfx_plot = plot(1:5)
    Plots._update_plot_object(pgfx_plot)
    @test pgfx_plot.o.the_plot isa PGFPlotsX.TikzDocument
    @test pgfx_plot.series_list[1].plotattributes[:quiver] === nothing

    @testset "3D docs example" begin
        n = 100
        ts = range(0, stop=8π, length=n)
        x = ts .* map(cos, ts)
        y = (0.1ts) .* map(sin, ts)
        z = 1:n
        pl = plot(x, y, z, zcolor=reverse(z), m=(10, 0.8, :blues, Plots.stroke(0)), leg=false, cbar=true, w=5)
        pgfx_plot = plot!(pl, zeros(n), zeros(n), 1:n, w=10)
        Plots._update_plot_object(pgfx_plot)
        if @test_nowarn(haskey(Plots.pgfx_axes(pgfx_plot.o)[1].options.dict, "colorbar") == true)
            @test Plots.pgfx_axes(pgfx_plot.o)[1]["colorbar"] === nothing
        end
     end # testset
     @testset "Color docs example" begin
        y = rand(100)
        plot(0:10:100, rand(11, 4), lab="lines", w=3, palette=:grays, fill=0, α=0.6)
        pl = scatter!(y, zcolor=abs.(y .- 0.5), m=(:heat, 0.8, Plots.stroke(1, :green)), ms=10 * abs.(y .- 0.5) .+ 4, lab="grad")
        Plots._update_plot_object(pl)
        axis = Plots.pgfx_axes(pl.o)[1]
        @test count( x->x isa PGFPlotsX.LegendEntry, axis.contents ) == 5
        @test count( x->x isa PGFPlotsX.Plot, axis.contents ) == 108 # each marker is its own plot, fillranges create 2 plot-objects
        marker = axis.contents[15]
        @test marker isa PGFPlotsX.Plot
        @test marker.options["mark"] == "*"
        @test marker.options["mark options"]["color"] == RGBA{Float64}( colorant"green", 0.8)
        @test marker.options["mark options"]["line width"] == 1
     end # testset
     @testset "Plot in pieces" begin
        plot(rand(100) / 3, reg=true, fill=(0, :green))
        scatter!(rand(100), markersize=6, c=:orange)
     end # testset
     @testset "Marker types" begin
        markers = filter((m->begin
           m in Plots.supported_markers()
           end), Plots._shape_keys)
        markers = reshape(markers, 1, length(markers))
        n = length(markers)
        x = (range(0, stop=10, length=n + 2))[2:end - 1]
        y = repeat(reshape(reverse(x), 1, :), n, 1)
        scatter(x, y, m=(8, :auto), lab=map(string, markers), bg=:linen, xlim=(0, 10), ylim=(0, 10))
     end # testset
     @testset "Layout" begin
        plot(Plots.fakedata(100, 10), layout=4, palette=[:grays :blues :heat :lightrainbow], bg_inside=[:orange :pink :darkblue :black])
     end # testset
     @testset "Polar plots" begin
        Θ = range(0, stop=1.5π, length=100)
        r = abs.(0.1 * randn(100) + sin.(3Θ))
        plot(Θ, r, proj=:polar, m=2)
     end # testset
     @testset "Drawing shapes" begin
        verts = [(-1.0, 1.0), (-1.28, 0.6), (-0.2, -1.4), (0.2, -1.4), (1.28, 0.6), (1.0, 1.0), (-1.0, 1.0), (-0.2, -0.6), (0.0, -0.2), (-0.4, 0.6), (1.28, 0.6), (0.2, -1.4), (-0.2, -1.4), (0.6, 0.2), (-0.2, 0.2), (0.0, -0.2), (0.2, 0.2), (-0.2, -0.6)]
         x = 0.1:0.2:0.9
         y = 0.7 * rand(5) .+ 0.15
         plot(x, y, line=(3, :dash, :lightblue), marker=(Shape(verts), 30, RGBA(0, 0, 0, 0.2)), bg=:pink, fg=:darkblue, xlim=(0, 1), ylim=(0, 1), leg=false)
     end # testset
     @testset "Histogram 2D" begin
        histogram2d(randn(10000), randn(10000), nbins=20)
     end # testset
     @testset "Heatmap-like" begin
        xs = [string("x", i) for i = 1:10]
        ys = [string("y", i) for i = 1:4]
        z = float((1:4) * reshape(1:10, 1, :))
        pgfx_plot = heatmap(xs, ys, z, aspect_ratio=1)
        Plots._update_plot_object(pgfx_plot)
        if @test_nowarn(haskey(Plots.pgfx_axes(pgfx_plot.o)[1].options.dict, "colorbar") == true)
           @test Plots.pgfx_axes(pgfx_plot.o)[1]["colorbar"] === nothing
           @test Plots.pgfx_axes(pgfx_plot.o)[1]["colormap name"] == "plots1"
        end

        pgfx_plot = wireframe(xs, ys, z, aspect_ratio=1)
        # TODO: clims are wrong
     end # testset
     @testset "Contours" begin
        x = 1:0.5:20
        y = 1:0.5:10
        f(x, y) = begin
           (3x + y ^ 2) * abs(sin(x) + cos(y))
        end
        X = repeat(reshape(x, 1, :), length(y), 1)
        Y = repeat(y, 1, length(x))
        Z = map(f, X, Y)
        p2 = contour(x, y, Z)
        p1 = contour(x, y, f, fill=true)
        plot(p1, p2)
        # TODO: colorbar for filled contours
     end # testset
     @testset "Varying colors" begin
        t = range(0, stop=1, length=100)
         θ = (6π) .* t
         x = t .* cos.(θ)
         y = t .* sin.(θ)
         p1 = plot(x, y, line_z=t, linewidth=3, legend=false)
         p2 = scatter(x, y, marker_z=((x, y)->begin
                             x + y
                         end), color=:bluesreds, legend=false)
         plot(p1, p2)
     end # testset
     @testset "Framestyles" begin
        scatter(fill(randn(10), 6), fill(randn(10), 6), framestyle=[:box :semi :origin :zerolines :grid :none], title=[":box" ":semi" ":origin" ":zerolines" ":grid" ":none"], color=permutedims(1:6), layout=6, label="", markerstrokewidth=0, ticks=-2:2)
        # TODO: support :semi
     end # testset
     @testset "Quiver" begin
        x = -2pi:0.2:2*pi
        y = sin.(x)

        u = ones(length(x))
        v = cos.(x)
        plot( x, y, quiver = (u, v), arrow = true )
        # TODO: could adjust limits to fit arrows if too long, but how?
     end # testset
     @testset "Annotations" begin
        y = rand(10)
        plot(y, annotations=(3, y[3], Plots.text("this is \\#3", :left)), leg=false)
        annotate!([(5, y[5], Plots.text("this is \\#5", 16, :red, :center)), (10, y[10], Plots.text("this is \\#10", :right, 20, "courier"))])
        scatter!(range(2, stop=8, length=6), rand(6), marker=(50, 0.2, :orange), series_annotations=["series", "annotations", "map", "to", "series", Plots.text("data", :green)])
     end # testset
end # testset


text:
using ArrayInterface

function ArrayInterface.restructure(x::AbstractArray{T1, 0}, y::AbstractArray{T2, 0}) where {T1, T2}
    out = similar(x, eltype(y))
    out .= y
    out
end


text:
const PATCHES = PatchStore(Dict(), Dict())

default_patch_store() = PATCHES

"""
    find(store::PatchStore, args...)

Find a patch from the `store` given the argument types.
"""
function find(store::PatchStore, args...)
    @debug "finding patch" dctkeys=collect(keys(store.dct))[1] args
    patch = get(store.dct, args, nothing)
    patch !== nothing && @debug "found patch" store.dct patch
    return patch
end

"""
    register(store::PatchStore, f::Callable, args)

Register a patch `f` in the `store` for the respective argument types.
"""
function register(store::PatchStore, f::Callable, args)
    # @info "Register" f args
    store.dct[args] = f
    return nothing
end

"""
    preserve(store::PatchStore)

Make a copy of the current state in the store.
"""
function preserve(store::PatchStore)
    empty!(store.prev)
    copy!(store.prev, store.dct)
end

"""
    restore(store::PatchStore)

Restore the previous state in the store.
"""
function restore(store::PatchStore)
    empty!(store.dct)
    copy!(store.dct, store.prev)
end