text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

using FMI
using Flux
using DifferentialEquations: Tsit5

import Random
Random.seed!(1234);

modelFMUPath = joinpath(dirname(@__FILE__), "..", "model", "SpringPendulum1D.fmu")
realFMUPath = joinpath(dirname(@__FILE__), "..", "model", "SpringFrictionPendulum1D.fmu")

t_start = 0.0
t_step = 0.01
t_stop = 5.0
tData = t_start:t_step:t_stop

# generate training data
realFMU = fmiLoad(realFMUPath)
fmiInstantiate!(realFMU; loggingOn=false)
fmiSetupExperiment(realFMU, t_start, t_stop)
fmiEnterInitializationMode(realFMU)
fmiExitInitializationMode(realFMU)
x0 = fmiGetContinuousStates(realFMU)
_, realSimData = fmiSimulateCS(realFMU, t_start, t_stop; recordValues=["mass.s", "mass.v"], setup=false, reset=false, saveat=tData)

# load FMU for NeuralFMU
myFMU = fmiLoad(modelFMUPath)
fmiInstantiate!(myFMU; loggingOn=false)
fmiSetupExperiment(myFMU, t_start, t_stop)
fmiEnterInitializationMode(myFMU)
fmiExitInitializationMode(myFMU)

# setup traing data
posData = collect(data[1] for data in realSimData.saveval)
#velData = collect(data[2] for data in realSimData.saveval)

# loss function for training
function losssum()
    global problem, x0, posData
    solution = problem(x0)

    posNet = collect(data[2] for data in solution.u)

    Flux.Losses.mse(posNet, posData)
end

# callback function for training
global iterCB = 0
global lastLoss = 0.0
function callb()
    global iterCB += 1
    global lastLoss

    if iterCB % 30 == 0
        loss = losssum()
        @info "Loss: $loss"
        @test loss < lastLoss
        lastLoss = loss
    end
end

vr = fmi2StringToValueReference(myFMU, "mass.m")

numStates = fmiGetNumberOfStates(myFMU)

# some NeuralFMU setups
nets = []

1 # default ME-NeuralFMU (learn dynamics and states, almost-neutral setup, parameter count << 100)
net = Chain(Dense( [1.0 0.0; 0.0 1.0] + rand(numStates,numStates)*0.01, zeros(numStates), tanh),
            states ->  fmiEvaluateME(myFMU, states),
            Dense( [1.0 0.0; 0.0 1.0] + rand(numStates,numStates)*0.01, zeros(numStates), identity))
push!(nets, net)

# 2 # default ME-NeuralFMU (learn dynamics)
net = Chain(states ->  fmiEvaluateME(myFMU, states),
            Dense(numStates, 16, tanh),
            Dense(16, 16, tanh),
            Dense(16, numStates))
push!(nets, net)

# 3 # default ME-NeuralFMU (learn states)
net = Chain(Dense(numStates, 16, leakyrelu),
            Dense(16, 16, leakyrelu),
            Dense(16, numStates),
            states -> fmiEvaluateME(myFMU, states))
push!(nets, net)

# 4 # default ME-NeuralFMU (learn dynamics and states)
net = Chain(Dense(numStates, 16, leakyrelu),
            Dense(16, 16, leakyrelu),
            Dense(16, numStates),
            states -> fmiEvaluateME(myFMU, states),
            Dense(numStates, 16, tanh),
            Dense(16, 16, tanh),
            Dense(16, numStates))
push!(nets, net)

# 5 # NeuralFMU with hard setting time to 0.0
net = Chain(states ->  fmiEvaluateME(myFMU, states), # not supported by this FMU:   states ->  fmiEvaluateME(myFMU, states, 0.0),
            Dense(numStates, 8, tanh),
            Dense(8, 16, tanh),
            Dense(16, numStates))
push!(nets, net)

# 6 # NeuralFMU with additional getter
getVRs = [fmi2StringToValueReference(myFMU, "mass.m")]
numGetVRs = length(getVRs)
net = Chain(states ->  fmiEvaluateME(myFMU, states, -1.0, fmi2ValueReference[], Real[], getVRs),
            Dense(numStates+numGetVRs, 8, tanh),
            Dense(8, 16, tanh),
            Dense(16, numStates))
push!(nets, net)

# 7 # NeuralFMU with additional setter
setVRs = [fmi2StringToValueReference(myFMU, "mass.m")]
numSetVRs = length(setVRs)
net = Chain(states ->  fmiEvaluateME(myFMU, states, -1.0, setVRs, [1.1]),
            Dense(numStates, 8, tanh),
            Dense(8, 16, tanh),
            Dense(16, numStates))
push!(nets, net)

# 8 # NeuralFMU with additional setter and getter
net = Chain(states ->  fmiEvaluateME(myFMU, states, -1.0, setVRs, [1.1], getVRs),
            Dense(numStates+numGetVRs, 8, tanh),
            Dense(8, 16, tanh),
            Dense(16, numStates))
push!(nets, net)

optim = ADAM(1e-4)
for i in 1:length(nets)
    @testset "Net setup #$i" begin
        global nets, problem, lastLoss, iterCB
        net = nets[i]
        problem = ME_NeuralFMU(myFMU, net, (t_start, t_stop), Tsit5(); saveat=tData)
        @test problem != nothing

        solutionBefore = problem(x0)
        @test length(solutionBefore.t) == length(tData)
        @test solutionBefore.t[1] == t_start
        @test solutionBefore.t[end] == t_stop

        # train it ...
        p_net = Flux.params(problem)

        iterCB = 0
        lastLoss = losssum()
        @info "Start-Loss for net #$i: $lastLoss"
        Flux.train!(losssum, p_net, Iterators.repeated((), 60), optim; cb=callb)

        # check results
        solutionAfter = problem(x0)
        @test length(solutionAfter.t) == length(tData)
        @test solutionAfter.t[1] == t_start
        @test solutionAfter.t[end] == t_stop
    end
end

fmiUnload(realFMU)
fmiUnload(myFMU)


text:
#  MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization
#  Copyright (C) 2020-2021, UChicago Argonne, LLC. All rights reserved.
#  Released under the modified BSD license. See COPYING.md for more details.

using Distributed
using JLD2

function LearningSolver(
    optimizer_factory;
    components = nothing,
    mode::AbstractString = "exact",
    simulate_perfect::Bool = false,
    solve_lp::Bool = true,
    extract_sa::Bool = true,
    extract_lhs::Bool = true,
)::LearningSolver
    return LearningSolver(
        miplearn.LearningSolver(
            solver = JuMPSolver(optimizer_factory),
            mode = mode,
            solve_lp = solve_lp,
            simulate_perfect = simulate_perfect,
            components = components,
            extract_lhs = extract_lhs,
            extract_sa = extract_sa,
        ),
        optimizer_factory,
    )
end

function solve!(
    solver::LearningSolver,
    instance::Instance;
    tee::Bool = false,
    discard_output::Bool = false,
)
    return @python_call solver.py.solve(
        instance.py,
        tee = tee,
        discard_output = discard_output,
    )
end

function fit!(solver::LearningSolver, instances::Vector{<:Instance})
    @python_call solver.py.fit([instance.py for instance in instances])
    return
end

function save(filename::AbstractString, solver::LearningSolver)
    internal_solver = solver.py.internal_solver
    internal_solver_prototype = solver.py.internal_solver_prototype
    solver.py.internal_solver = nothing
    solver.py.internal_solver_prototype = nothing
    solver_py_filename = tempname()
    miplearn.write_pickle_gz(solver.py, solver_py_filename)
    solver_py = read(solver_py_filename)
    solver.py.internal_solver = internal_solver
    solver.py.internal_solver_prototype = internal_solver_prototype
    jldsave(
        filename;
        miplearn_version = "0.2",
        solver_py = solver_py,
        optimizer_factory = solver.optimizer_factory,
    )
    return
end

function load_solver(filename::AbstractString)::LearningSolver
    jldopen(filename, "r") do file
        solve_py_filename = tempname()
        write(solve_py_filename, file["solver_py"])
        solver_py = miplearn.read_pickle_gz(solve_py_filename)
        internal_solver = JuMPSolver(file["optimizer_factory"])
        solver_py.internal_solver_prototype = internal_solver
        return LearningSolver(solver_py, file["optimizer_factory"])
    end
end

export Instance, LearningSolver, solve!, fit!, parallel_solve!, save, load_solver

ext: jl


text:
# GTAS but with negative spikes !

abstract type AntiKernel end

# negative correlations for pairs only!

struct pAGTAS{R,I,J<:Jittering,AK<:AntiKernel}
  gtas::GTAS{R,I,J}
  antimarkings::Vector{Tuple{I,I}}
  antiprobabilities::Vector{R}
  antikernels::Vector{AK}
end

function pAGTAS(rate_parent::R,markings::Vector{Vector{I}},
    marking_probs::Vector{R},jitters::Vector{J},
    antimarkings::Vector{Tuple{I,I}},
    antiprobabilities::Vector{R},
    antikernels::Vector{AK}) where {R<:Real,I<:Integer,J<:Jittering,AK<:AntiKernel}
  @assert length(antimarkings) == length(antiprobabilities)
  @assert length(antikernels) == length(antiprobabilities)
  @assert all(<=(1.0),antiprobabilities)
  gtas = GTAS(rate_parent,markings,marking_probs,jitters)
  return pAGTAS{R,I,J,AK}(gtas,antimarkings,antiprobabilities,antikernels)
end

n_units(agtas::pAGTAS) = agtas.gtas.n
n_units(gtas::GTAS) = gtas.n

struct AntiExponential <: AntiKernel
  τ::Float64
end

function compute_forward_cutprob(t_now::R,antiker::AntiExponential,
    train::Vector{R},Nmax::Integer,Cmax::R;dowarn::Bool=true) where R
  idxstart = searchsortedfirst(train,t_now+eps())
  idx = idxstart
  ret = Vector{Float64}(undef,Nmax)
  d = Exponential(antiker.τ)
  for k in 1:Nmax
    _cdf = cdf(d,train[idx]-t_now)
    ret[k] = _cdf
    idx+=1
    if _cdf >= Cmax
      break
    end
    if ! checkbounds(Bool,train,idx)
      @error "Not enough spikes left ! More spikes needed!"
      break
    end
  end
  nkeep = idx - idxstart
  keepat!(ret,1:nkeep)
  if (nkeep == Nmax) && dowarn
    @warn "the cdf is $(ret[end]) and not $Cmax !"
  end
  return (idxstart,ret)
end

function apply_antispike!(idx_start::Integer,cdfs::Vector{R},train::Vector{R}) where R
  umax = cdfs[end]
  u = rand()*umax
  idx_kill = idx_start
  for cdf in cdfs
    if u<=cdf
      break
    else
      idx_kill +=1
    end
  end
  deleteat!(train,idx_kill)
  return nothing
end

function apply_antispike!(t_now::R,antiker::AntiExponential,
      train::Vector{R};Nmax::Integer=200,Cmax::Float64=0.999) where R
  idx_start,cdfs = compute_forward_cutprob(t_now,antiker,train,Nmax,Cmax)
  apply_antispike!(idx_start,cdfs,train)
  return nothing
end

function make_samples(g::pAGTAS,t_tot::Real)
  trains = make_samples(g.gtas,1.5*t_tot)
  for ((pre,post),p,antiker) in zip(g.antimarkings,g.antiprobabilities,g.antikernels)
    trainpre = trains[pre]
    trainpost = trains[post]
    for t_now in trainpre
      if t_now > t_tot
        break
      elseif rand() < p
        apply_antispike2!(t_now,antiker,trainpost)
      end
    end
  end
  # remove excess time
  for train in trains
    idx = searchsortedfirst(train,t_tot)
    keepat!(train,1:idx-1)
  end
  return trains
end

function compute_forward_cutprobabilities(t_now::R,antiker::AntiKernel,
    train::Vector{R},eps_prob::R) where R
  idxstart = searchsortedfirst(train,t_now+eps())
  t_last = t_now + antikernel_horizon(antiker,eps_prob)
  idxend = searchsortedfirst(train,t_last)-1
  if idxend == length(train)-1
    @error "Not enough spikes left ! More spikes needed!"
  end
  ret = compute_forward_cutprobabilities(t_now,view(train,idxstart:idxend),antiker)
  return (idxstart,ret)
end

@inline function antikernel_horizon(antiker::AntiExponential,eps_prob::Float64)
  τ = antiker.τ
  return - τ * (log(τ)+log(eps_prob))
end

@inline function compute_forward_cutprobabilities(t_now::Real,train::SubArray{Float64},antiker::AntiExponential)
  d = Exponential(antiker.τ)
  ret = map(t -> pdf(d,t-t_now),train)
  ret ./= sum(ret)
  return ret
end

function apply_antispike2!(t_now::R,antiker::AntiExponential,
      train::Vector{R};eps_prob::Float64=1E-4) where R
  idx_start,cutprobs = compute_forward_cutprobabilities(t_now,antiker,train,eps_prob)
  idx_cut = rand(Categorical(cutprobs))
  deleteat!(train,idx_start+idx_cut-1)
  return nothing
end

##

# sequential antispikes

abstract type AbstractMarking end
abstract type AbstractAntiJittering <: Jittering end

struct Marking <: AbstractMarking
  vals::Vector{Int64}
end
struct AntiMarking <: AbstractMarking
  vals::Vector{Int64}
end

struct AntiJitterExp <: AbstractAntiJittering
  τ::Float64
end
struct AntiJitterExpSequential <: AbstractAntiJittering
  τ::Float64
end
struct AntiJitterStepSequential <: AbstractAntiJittering
  τ::Float64
end

struct sAGTAS{N,NTR<:NTuple{N,Float64},
        NTM<:NTuple{N,AbstractMarking},
        NTJ<:NTuple{N,Jittering}}
  n::Int64
  rate_markings::NTR
  markings::NTM
  jitters::NTJ
end
nmarkings(agta::sAGTAS) = length(agta.markings)

function sAGTAS(n::Int64,
    rate_markings::Vector{Float64},
    markings::Vector{M},
    jitters::Vector{J}) where {M<:AbstractMarking,J<:Jittering}
  @assert length(markings) == length(jitters)
  @assert length(rate_markings) == length(jitters)
  return sAGTAS(n,Tuple(rate_markings),Tuple(markings),Tuple(jitters))
end

function make_samples(g::sAGTAS,t_tot::R) where R<:Real
  # initialize trains
  trains = [Vector{Float64}(undef,0) for _ in 1:g.n]
  # one train for each marking
  # (saved to keep antimarkings for backwards move)
  t_forw = t_tot*1.1+10.0
  trains_markings = map(g.rate_markings) do rat
    make_poisson_samples(rat,t_forw)
  end
  # forward move
  for (train_marking,marking,jitt) in zip(trains_markings,g.markings,g.jitters)
    add_marking_to_trains!(trains,train_marking,marking,jitt)
  end
  # clean up trains: sorted, and > 0
  sort!.(trains)
  tmin = minimum(first.(trains))
  if tmin < 0
    for tr in trains
      tr .-= tmin
    end
  end
  # backwards now!
  for (train_marking,marking,jitt) in zip(trains_markings,g.markings,g.jitters)
    remove_antimarking_from_trains!(trains,train_marking,marking,jitt,t_tot)
  end
  # stop time right before t_tot
  for train in trains
    # filter!(<=(t_tot),train)
    k = searchsortedfirst(train,t_tot)
    keepat!(train,1:(k-1))
  end
  return trains
end

# do nothing for antimarking
function add_marking_to_trains!(::Vector,
    ::Vector,::AntiMarking,jitter::Jittering)
  return nothing
end

# do the usual for markings
function add_marking_to_trains!(trains::Vector{Vector{R}},
    trainmark::Vector{R},marking::Marking,jitter::Jittering) where R<:Real
  mark = marking.vals
  n = length(mark)
  for t_k in trainmark
    spikes_k = jitter!(fill(t_k,n),jitter)
    for (tkk,kk) in zip(spikes_k,mark)
     push!(trains[kk],tkk)
    end
  end
  return nothing
end

# do nothing for markings
function remove_antimarking_from_trains!(::Vector,
    ::Vector,::Marking,::Jittering,::Real)
  return nothing
end

@inline function antijitter_horizon(
    antijitt::Union{AntiJitterExp,AntiJitterExpSequential},
    eps_prob::Float64)
  τ = antijitt.τ
  return - τ * (log(τ)+log(eps_prob))
end

@inline function antijitter_horizon(antijitt::AntiJitterStepSequential,::Float64)
  return antijitter_horizon(antijitt)
end
@inline function antijitter_horizon(antijitt::AntiJitterStepSequential)
  return antijitt.τ
end

function remove_antimarking_from_trains!(trains::Vector{Vector{R}},
    trainmark::Vector{R},anti::AntiMarking,
    antijitter::Union{AntiJitterExpSequential,AntiJitterStepSequential},
    t_tot::R) where R<:Real
  mark = anti.vals
  # remove excess time from mark train
  kmax = searchsortedfirst(trainmark,t_tot)
  keepat!(trainmark,1:(kmax-1))
  for t_k in trainmark
    # find and remove next spike in first train
    train1 = trains[mark[1]]
    idx_cut = searchsortedfirst(train1,t_k)
    deleteat!(train1,idx_cut)
    if length(mark) > 1
      tnow = train1[idx_cut]
      t_incr =  antijitter_horizon(antijitter,1E-5)
      # in following trains (if any) remove with exp probability
      for ms in mark[2:end]
        trainh = trains[ms]
        tmax = tnow + t_incr
        idx_start,cutprobs =
          compute_forward_killprobabilities(tnow,tmax,trainh,antijitter)
        idx_cut = idx_start + rand(Categorical(cutprobs)) - 1
        tnow = trainh[idx_cut]
        deleteat!(trainh,idx_cut)
      end
    end
  end
  return nothing
end

function remove_antimarking_from_trains!(trains::Vector{Vector{R}},
    trainmark::Vector{R},anti::AntiMarking,
    antijitter::AntiJitterExp,t_tot::R) where R<:Real
  mark = anti.vals
  # remove excess time from mark train
  kmax = searchsortedfirst(trainmark,t_tot)
  keepat!(trainmark,1:(kmax-1))
  # time horizon
  t_incr =  antijitter_horizon(antijitter,1E-5)
  for t_k in trainmark
    tnow = t_k
    tmax = t_k + t_incr
    for ms in mark
      trainh = trains[ms]
      idx_start,cutprobs =
          compute_forward_killprobabilities(tnow,tmax,trainh,antijitter)
      idx_cut = idx_start + rand(Categorical(cutprobs)) - 1
      deleteat!(trainh,idx_cut)
    end
  end
  return nothing
end

# returns idx of closest time >= t_start and probabilities of killing spike
function compute_forward_killprobabilities(t_start::R,t_end::R,train::Vector{R},
     antijitter::Union{AntiJitterExpSequential,AntiJitterExp}) where R<:Real
  idx_start =  searchsortedfirst(train,t_start)
  idx_end = searchsortedfirst(train,t_end)
  ret = map(t -> exp(-(t-t_start)/antijitter.τ), view(train,idx_start:idx_end) )
  ret ./= sum(ret)
  return idx_start,ret
end

function compute_forward_killprobabilities(t_start::R,t_end::R,train::Vector{R},
     ::AntiJitterStepSequential) where R<:Real
  idx_start =  searchsortedfirst(train,t_start)
  idx_end = searchsortedfirst(train,t_end)-1 # smaller than t_end
  n_spikes = max(idx_end-idx_start+1,1) # if no spikes in interval, just erase next one :-(
  return idx_start,fill(inv(n_spikes),n_spikes)
end

function remove_next!(train::Vector{Float64},antitrain::Vector{Float64})
  for i in eachindex(antitrain)
    remove_next!(train,antitrain[i])
  end
  return nothing
end

function remove_next!(train::Vector{Float64},t::Float64)
  if  t < train[end]
    idx_cut = searchsortedfirst(train,t)
    deleteat!(train,idx_cut)
  end
  return nothing
end


text:
module PkgServer

using Pkg
using HTTP
using Base.Threads: Event, @spawn
using Random
using LibGit2

include("resource.jl")
include("meta.jl")

function __init__()
    # Set default HTTP useragent
    HTTP.setuseragent!("PkgServer (HTTP.jl)")
end

function start(;host="127.0.0.1", port=8000)
    mkpath("cache")
    update_registries()
    @sync begin
        @spawn while true
            sleep(1)
            forget_failures()
            update_registries()
        end
        @info "server listening on $(host):$(port)"
        HTTP.listen(host, port) do http
            resource = http.message.target
            # If the user is asking for `/meta`, generate the requisite JSON object and send it back
            if occursin(meta_re, resource)
                serve_meta(http)
                return
            end

            # If the user asked for soemthing that is an actual resource, send it directly
            if occursin(resource_re, resource)
                path = fetch(resource)
                if path !== nothing
                    serve_file(http, path)
                    return
                end
            end
            HTTP.setstatus(http, 404)
            startwrite(http)
        end
    end
end

end # module


text:
"""
    FEMMDeforLinearModule

Module for operations on interiors of domains to construct system matrices and
system vectors for linear deformation models.
"""
module FEMMDeforLinearModule

__precompile__(true)

using FinEtools.FTypesModule: FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict
import FinEtools.FENodeSetModule: FENodeSet
import FinEtools.FESetModule: AbstractFESet, manifdim
import FinEtools.IntegDomainModule: IntegDomain
import FinEtoolsDeforLinear.FEMMDeforLinearBaseModule: AbstractFEMMDeforLinear
import FinEtoolsDeforLinear.DeforModelRedModule: AbstractDeforModelRed, DeforModelRed2DAxisymm
import FinEtoolsDeforLinear.MatDeforLinearElasticModule: AbstractMatDeforLinearElastic
import FinEtools.CSysModule: CSys

"""
    FEMMDeforLinear{MR<:AbstractDeforModelRed,  S<:AbstractFESet, F<:Function, M<:AbstractMatDeforLinearElastic} <: AbstractFEMMDeforLinear

Class for linear deformation finite element modeling machine.
"""
mutable struct FEMMDeforLinear{MR<:AbstractDeforModelRed,  S<:AbstractFESet, F<:Function, M<:AbstractMatDeforLinearElastic} <: AbstractFEMMDeforLinear
    mr::Type{MR} # model reduction type
    integdomain::IntegDomain{S, F} # integration domain data
    mcsys::CSys # updater of the material orientation matrix
    material::M # material object
end
"""
    FEMMDeforLinear(mr::Type{MR}, integdomain::IntegDomain{S, F}, material::M) where {MR<:AbstractDeforModelRed, S<:AbstractFESet, F<:Function, M<:AbstractMatDeforLinearElastic}

Constructor of linear deformation finite element modeling machine.
"""
function FEMMDeforLinear(mr::Type{MR}, integdomain::IntegDomain{S, F}, material::M) where {MR<:AbstractDeforModelRed, S<:AbstractFESet, F<:Function, M<:AbstractMatDeforLinearElastic}
    @assert mr == material.mr "Model reduction is mismatched"
    @assert (integdomain.axisymmetric) || (mr != DeforModelRed2DAxisymm) "Axially symmetric requires axisymmetric to be true"
    return FEMMDeforLinear(mr, integdomain, CSys(manifdim(integdomain.fes)), material)
end

end


text:
using NewtonsMethod
using Test

#Test method 1
f(x)=x^2-4.0
f′(x)=2x
@test newtonroot(f,f′,x₀=0.2)[1]≈2.0

f(x)=log(x)-20
f′(x)=1/x
@test newtonroot(f,f′,x₀=0.2)[1]≈4.851651954097909e8

f(x)=3x^2-5x+1
f′(x)=6x-5
@test newtonroot(f,f′,x₀=0.2)[1]≈2.32408120756001784480129727405038536389084291259026547774006206530426035584922e-01

f(x)=-10x^3-5x^2+20
f′(x)=-20x^2-10x
@test newtonroot(f,f′,x₀=0.2)[1]≈1.113386205687796447877673190361708447749219094960359587626412934254462034979434

#Test method 2
f(x)=x^2-4.0
@test newtonroot(f,x₀=0.2)[1]≈2.0

f(x)=log(x)-20
@test newtonroot(f,x₀=0.2)[1]≈4.851651954097909e8

f(x)=3x^2-5x+1
@test newtonroot(f,x₀=0.2)[1]≈2.32408120756001784480129727405038536389084291259026547774006206530426035584922e-01

f(x)=-10x^3-5x^2+20
@test newtonroot(f,x₀=0.2)[1]≈1.113386205687796447877673190361708447749219094960359587626412934254462034979434

#Test Big Float (Not really sure what to do here)
#(I set all parameters to be BigFloat)
f(x)=3x^2-5x+1
@test newtonroot(f,x₀= BigFloat(0.2),tolerance = BigFloat(1E-100), maxiter = BigFloat(1000))[1]  ≈  2.32408120756001784480129727405038536389084291259026547774006206530426035584922e-01

#Test non-convergence (return nothing)
f(x)=2+x^2
@test newtonroot(f,x₀=0.2)==nothing

#Test maxiter (first should return the value, second should return nothing)
f(x)=log(x)-20
a=newtonroot(f,x₀=0.2)[1] #Algorithm needs 17 iterations in this case
b=newtonroot(f,x₀=0.2,maxiter=5)
@testset "maxiter" begin
  @test a≈4.851651954097909e8
  @test b==nothing
 end;

#Test tolerance (As tolerance parameter increases, the code is less accurate)
f(x)=3x^2-5x+1
a=newtonroot(f,x₀=1)[1]
b=newtonroot(f,x₀=1, tolerance=0.0005)[1]
c=newtonroot(f,x₀=1, tolerance=0.05)[1]
@test f(a)<f(b)<f(c)


text:
"""

A structure containing the
"""
mutable struct SubSolvers{Q<:MOI.AbstractOptimizer, S<:MOI.AbstractOptimizer, T<:ExtensionType}
    relaxed_optimizer::Q
    upper_optimizer::S
    ext::T
end
SubSolvers{Q,S,T}(; r::Q = Cbc.Optimizer(), u::S = Ipopt.Optimizer(), t::T = DefaultExt()) where {Q,S,T} = SubSolvers{Q,S,T}(r, u, t)
SubSolvers(; r::Q = Cbc.Optimizer(), u::S = Ipopt.Optimizer(), t::T = DefaultExt()) where {Q,S,T} = SubSolvers{Q,S,T}(r,u,t)

function _relaxed_optimizer(d::SubSolvers{Q,S,T}) where {Q <: MOI.AbstractOptimizer,
                                                         S <: MOI.AbstractOptimizer,
                                                         T <: ExtensionType}
    return d.relaxed_optimizer
end

function _upper_optimizer(d::SubSolvers{Q,S,T}) where {Q <: MOI.AbstractOptimizer,
                                                       S <: MOI.AbstractOptimizer,
                                                       T <: ExtensionType}
    return d.upper_optimizer
end

function _ext(d::SubSolvers{Q,S,T}) where {Q <: MOI.AbstractOptimizer,
                                                S <: MOI.AbstractOptimizer,
                                                T <: ExtensionType}
    return d.ext
end

function isempty(d::SubSolvers{Q,S,T}) where {Q,S,T}
    MOI.is_empty(d.relaxed_optimizer) && MOI.is_empty(d.upper_optimizer)
end
function MOI.empty!(d::SubSolvers{Q,S,T}) where {Q,S,T}
    MOI.empty!(d.relaxed_optimizer)
    MOI.empty!(d.upper_optimizer)
    return
end

max_stars_repo_name: MatthewStuber/EAGO

max_issues_repo_name: MatthewStuber/EAGO

max_forks_repo_name: MatthewStuber/EAGO


text:
minutes = 105

hours = minutes/60 # 1.75
println(hours)

# Floor Division
hours = minutes ÷ 60 # 1
println(hours)

# Modulus
remainder = minutes % 60 # 45
println(remainder)

typeof(true)
typeof(Bool)

x = 5
if 0 < x < 10
    println("This is cool")
end

# Keyboard Input
line = readline()
typeof(line)


text:
module MolProperties

include("penalized-logp.jl")

end # module

max_stars_repo_name: Soyukke/MolProperties

max_issues_repo_name: Soyukke/MolProperties

max_forks_repo_name: Soyukke/MolProperties


text:
"Implementa a análise estrutural."
module FEA

using LinearAlgebra,
      SparseArrays,
      NearestNeighbors

export No,
       Material,
       Elemento,
       Estrutura,
       comprimento,
       num_gls_estrutura,
       vetor_forcas,
       volume,
       deslocamentos,
       gls,
       k,
       tem_sobreposicao,
       carregar_apoiar_no!

"Número de graus de liberdade de um nó."
const num_gls_no = 2

# =======================================================================
"""
    No(id::Integer,
       x::Real,
       y::Real;
       fx::Real = 0,
       fy::Real = 0,
       rx::Bool = false,
       ry::Bool = false)

Implementa as propriedades de um nó.

# Argumentos
- `id::Integer`: Identificação do nó.
- `x::Real`: Coordenada x.
- `y::Real`: Coordenada y.
- `fx::Real`: Força na direção x.
- `fy::Real`: Força na direção y.
- `rx::Bool`: Apoio na direção x. true para impedido e 'false' para livre.
- `ry::Bool`: Apoio na direção y. true para impedido e 'false' para livre.
"""
mutable struct No
    id::Int64
    x::Float64
    y::Float64
    fx::Float64
    fy::Float64
    rx::Bool
    ry::Bool

    function No(id::Integer,
                x::Real,
                y::Real;
                fx::Real = 0,
                fy::Real = 0,
                rx::Bool = false,
                ry::Bool = false)

        if id < 1
            throw(ArgumentError("O número de identificação do nó deve ser maior que 1!"))
        end
        new(id, x, y, fx, fy, rx, ry)
    end
end

"""
    distancia(no₁::No, no₂::No)::Float64

Retorna a distância entre `no₁` e `no₂`.
"""
distancia(no₁::No, no₂::No)::Float64 = norm([(no₂.x - no₁.x), (no₂.y - no₁.y)])

"""
    gls(no::No)::Vector{Int64}

Retorna os graus de liberdade do nó.
"""
gls(no::No)::Vector{Int64} = [num_gls_no * no.id - 1, num_gls_no * no.id]

"""
    vetor_forcas(no::No)::Vector{Float64}

Retorna as forças nodais.
"""
vetor_forcas(no::No)::Vector{Float64} = [no.fx, no.fy]

"""
    vetor_apoios(no::No)::Vector{Bool}

    Retorna os apoios nodais.
"""
vetor_apoios(no::No)::Vector{Bool} = [no.rx, no.ry]

"""
    gls_livres(no::No)::Vector{Int64}

Retorna os graus de liberdade livres do nó.
"""
gls_livres(no::No)::Vector{Int64} = [i for (i, j) in zip(gls(no), vetor_apoios(no)) if j == false]

"""
    sao_colineares(no₁::No, no₂::No, no₃::No)::Bool

Verifica se três nós são colineares.

É feito o cálculo da área do triângulo formado pelos três nós. Caso essa área seja nula,
os nós são colineares.
"""
function sao_colineares(no₁::No, no₂::No, no₃::No)::Bool
    m = [no₁.x no₁.y 1
        no₂.x no₂.y 1
        no₃.x no₃.y 1]

    return isapprox(det(m), 0; atol=1e-9)
end

# =====================================================================================
"""
Define as propriedades do material que compõe os elementos.

# Argumentos
- `v::Integer`: Coeficiente de Poisson.
- `E::AbstractFloat`: Módulo de Elasticidade.
"""
struct Material
    id::Int64
    E::Float64

    function Material(id::Integer, E::AbstractFloat)
        if id < 1
            throw(ArgumentError("O número de identificação do material não pode ser menor que 1!"))
        end
        if !(E > 0)
            throw(ArgumentError("O Módulo de Elasticidade deve ser maior que 0!"))
        end
        new(id, E)
    end
end

# =====================================================================================
"""
    Elemento(id::Integer, no₁::No, no₂::No, area::Real, material::Material)

Implementa as propriedades de uma elemento de treliça.

# Argumentos
- `id::Integer`: Número de identificação do elemento.
- `no₁::No`: Nó inicial.
- `no₂::No`: Nó final.
- `area::Real`: Área da seção transversal.
- `material::Material`: Material utilizado no elemento.
"""
mutable struct Elemento
    id::Int64
    no₁::No
    no₂::No
    area::Float64
    material::Material

    function Elemento(id::Integer, no₁::No, no₂::No, area::Real, material::Material)
        if id < 1
            throw(ArgumentError("O número de identificação do elemento não pode ser menor que 1!"))
        end
        if area < 0
            throw(ArgumentError("A área da seção não pode ser menor que 0!"))
        end
        new(id, no₁, no₂, area, material)
    end
end

"""
    comprimento(elemento::Elemento)::Float64

Retorna o comprimento do elemento.
"""
comprimento(elemento::Elemento)::Float64 = distancia(elemento.no₁, elemento.no₂)

"""
    gls(elemento::Elemento)::Vector{Int64}

Retorna os graus de liberdade do elemento.
"""
gls(elemento::Elemento)::Vector{Int64} = vcat(gls(elemento.no₁), gls(elemento.no₂))

"""
    vetor_forcas(elemento::Elemento)::Vector{Float64}

Retorna o vetor de forças atuantes nos graus de liberdade do elemento.
"""
vetor_forcas(elemento::Elemento)::Vector{Float64} = vcat(vetor_forcas(elemento.no₁), vetor_forcas(elemento.no₂))

"""
    gls_livres(elemento::Elemento)::Vector{Int64}

Retorna os graus de liberdade livres do elemento.
"""
gls_livres(elemento::Elemento)::Vector{Int64} = vcat(gls_livres(elemento.no₁), gls_livres(elemento.no₂))

"""
    angulo(elemento::Elemento)::Float64

Retorna o ângulo de inclinação do elemento em relação a horizontal.

Segue as convenções do círculo trigonométrico.
"""
function angulo(elemento::Elemento)::Float64
    dx = elemento.no₂.x - elemento.no₁.x
    dy = elemento.no₂.y - elemento.no₁.y

    return atan(dy, dx)
end

"""
    matriz_rotacao(elemento::Elemento)::Matrix{Float64}

Retorna a matriz de rotação do elemento.
"""
function matriz_rotacao(elemento::Elemento)::Matrix{Float64}
    θ = angulo(elemento)
    s = sin(θ)
    c = cos(θ)

    return [c s 0 0
           -s c 0 0
            0 0 c s
            0 0 -s c]
end

"""
    k_local(elemento::Elemento)::Matrix{Float64}

Retorna a matriz de rigidez local do elemento.
"""
function k_local(elemento::Elemento)::Matrix{Float64}
    L = comprimento(elemento)
    E = elemento.material.E
    A = elemento.area

    return (E * A / L) * [1 0 -1 0
                          0 0 0 0
                         -1 0 1 0
                          0 0 0 0]
end

"""
    k(elemento::Elemento)::Matrix{Float64}

Retorna a matriz de rigidez do elemento rotacionada no sistema global.
"""
function k(elemento::Elemento)::Matrix{Float64}
    k = k_local(elemento)
    r = matriz_rotacao(elemento)

    return r' * k * r
end

"""
    sao_colineares(elemento₁::Elemento, elemento₂::Elemento)::Bool

Verifica se dois elementos são colineares.
"""
function sao_colineares(elemento₁::Elemento, elemento₂::Elemento)::Bool
    el1 = elemento₁
    el2 = elemento₂
    if sao_colineares(el1.no₁, el1.no₂, el2.no₁) && sao_colineares(el1.no₁, el1.no₂, el2.no₂)
        return true
    else
        return false
    end
end

"""
    tem_sobreposicao(elemento₁::Elemento, elemento₂::Elemento)::Bool

Verifica se há sobreposição entre dois elementos.
"""
function tem_sobreposicao(elemento₁::Elemento, elemento₂::Elemento)::Bool
    # Verificar se os nós dos elementos são colineares.
    el1 = elemento₁
    el2 = elemento₂
    L1 = comprimento(el1)
    L2 = comprimento(el2)
    # Se os dois elementos possuírem os mesmos nós, há sobreposição.
    if el1.no₁ in [el2.no₁, el2.no₂] && el1.no₂ in [el2.no₁, el2.no₂]
        return true
    else
        if sao_colineares(el1, el2)
            ls = [distancia(el1.no₁, el2.no₁),
                  distancia(el1.no₁, el2.no₂),
                  distancia(el1.no₂, el2.no₁),
                  distancia(el1.no₂, el2.no₂)]

            # Distância entre os nós extremos dos dois elementos de referência.
            dmax = maximum(ls)

            if dmax ≤ (L1 + L2) && !isapprox(dmax, L1 + L2; atol=1e-9)
                return true
            else
                return false
            end
        else
            return false
        end
    end
end

# =====================================================================================
"""
    Estrutura(nos::Vector{No}, elementos::Vector{Elemento}, λ::AbstractFloat=1e-9)

Implementa as propriedades de uma Estrutura.

# Argumentos
- `nos::Array{No}`: Nós que compõem a estrutura.
- `elementos::Array{Elemento}`: Elementos que compõem a estrutura.
- `λ::Float64`: Parâmetro de Tikhonov.
"""
struct Estrutura
    nos::Vector{No}
    elementos::Vector{Elemento}
    λ::Float64

    function Estrutura(nos::Vector{No}, elementos::Vector{Elemento}, λ::AbstractFloat=1e-9)
        new(nos, elementos, λ)
    end
end

"""
    carregar_apoiar_no!(nos::Vector{No}, x, y; fx=0, fy=0, rx=false, ry=false)

Adiciona as cargas `fx` e `fy` ao nó mais próximo do ponto de coordenadas `x` e `y`.
"""
function carregar_apoiar_no!(nos::Vector{No}, x, y; fx=0, fy=0, rx=false, ry=false)
    # Matriz contendo os dados de todos os nós.
    m = Matrix{Float64}(undef, (2, length(nos)))
    for i=eachindex(nos)
        m[1, i] = nos[i].x
        m[2, i] = nos[i].y
    end
    kd = KDTree(m)
    no_p = knn(kd, [x, y], 1)[1][1]

    nos[no_p].fx = fx
    nos[no_p].fy = fy
    nos[no_p].rx = rx
    nos[no_p].ry = ry
end

"""
    num_gls_estrutura(estrutura::Estrutura)::Int64 = num_gls_no * length(estrutura.nos)

Retorna o número de graus de liberdade da estrutura.
"""
num_gls_estrutura(estrutura::Estrutura)::Int64 = num_gls_no * length(estrutura.nos)

"""
    gls_livres(estrutura::Estrutura)::Vector{Int64}

Retorna os graus de liberdade livres da estrutura.
"""
gls_livres(estrutura::Estrutura)::Vector{Int64} = vcat([gls_livres(el) for el in estrutura.elementos]...) |> unique |> sort

"""
    volume(estrutura::Estrutura)::Float64

Retorna o volume da estrutura.
"""
volume(estrutura::Estrutura)::Float64 = sum([el.area * comprimento(el) for el in estrutura.elementos])

"""
    vetor_forcas(estrutura::Estrutura, incluir_gls_impedidos::Bool = true)::Vector{Float64}

Retorna o vetor de forças da estrutura.
"""
function vetor_forcas(estrutura::Estrutura, incluir_gls_impedidos::Bool = true)::Vector{Float64}
    forcas = Array{Float64}(undef, num_gls_estrutura(estrutura))

    for no in estrutura.nos
        forcas[gls(no)] = [no.fx, no.fy]
    end

    if incluir_gls_impedidos
        return forcas
    else
        return forcas[gls_livres(estrutura)]
    end
end

"""
    vetor_apoios(estrutura::Estrutura)::Vector{Int64}

Retorna os graus de liberdade impedidos.
"""
function vetor_apoios(estrutura::Estrutura)::Vector{Int64}
    apoios::Array{Int64} = []

    for no in estrutura.nos
        vaps = vetor_apoios(no)

        if any(vaps)
            gl = gls(no)
            append!(apoios, [gl[i] for (i, j) in enumerate(vaps) if j])
        end
    end
    return apoios
end

"""
    k(estrutura::Estrutura)::SparseMatrixCSC

Retorna a matriz de rigidez esparsa da estrutura.

Procedimento de montagem com 3 vetores separados.
"""
function k(estrutura::Estrutura)::SparseMatrixCSC
    # Número de termos em cada matriz de rigidez de elementos.
    num_gle = (num_gls_no * 2)
    apoios = vetor_apoios(estrutura)

    num_it = length(estrutura.elementos) * num_gle^2

    linhas = ones(Int64, num_it)
    colunas = ones(Int64, num_it)
    termos = zeros(num_it)

    c = 1

    for elem in estrutura.elementos
        gle = gls(elem)
        k_glob = k(elem)

        # Posições livres do elemento.
        pos_livres = [i for (i, gl) in enumerate(gle) if gl ∉ apoios]

        if length(pos_livres) > 0
            for i in pos_livres
                lin = gle[i] - count(x -> x < gle[i], apoios)
                for j in pos_livres
                    linhas[c] = lin
                    colunas[c] = gle[j] - count(x -> x < gle[j], apoios)
                    termos[c] = i != j ? k_glob[i, j] : k_glob[i, j] + estrutura.λ
                    c += 1
                end
            end
        end
    end
    ks = sparse(linhas, colunas, termos)
    dropzeros!(ks)

    return ks
end

"""
    deslocamentos(estrutura::Estrutura, incluir_gls_impedidos::Bool = true)::Array{Float64}

Retorna os deslocamentos nodais.
"""
function deslocamentos(estrutura::Estrutura, incluir_gls_impedidos::Bool = true)::Array{Float64}
    u = k(estrutura) \ vetor_forcas(estrutura, false)
    if !incluir_gls_impedidos
        return u
    else
        uf = zeros(num_gls_estrutura(estrutura))
        uf[gls_livres(estrutura)] = u
        return uf
    end
end

"""
    nos_deslocados(estrutura::Estrutura, escala::Int64 = 1)::Vector{No}

Retorna um vetor contendo os nós em suas novas posições com a estrutura deformada.
"""
function nos_deslocados(estrutura::Estrutura, escala::Int64 = 1)::Vector{No}
    num_nos = length(estrutura.nos)
    u_pontos = reshape(deslocamentos(estrutura), (2, num_nos))'
    nos_desloc = Array{No}(undef, length(estrutura.nos))

    for i in 1:num_nos
        nos_desloc[i] = No(i, estrutura.nos[i].x + escala * u_pontos[i, 1],
            estrutura.nos[i].y + escala * u_pontos[i, 2])
    end
    return nos_desloc
end

"""
    forcas_internas(estrutura::Estrutura)::Vector{Float64}

Retorna os valores dos esforços normais atuantes em cada elemento.
"""
function forcas_internas(estrutura::Estrutura)::Vector{Float64}
    forcas = Array{Float64}(undef, length(estrutura.elementos))
    nos_def = nos_deslocados(estrutura)

    for (i, elem) in enumerate(estrutura.elementos)
        l₀ = comprimento(elem)
        l₁ = distancia(nos_def[elem.no₁.id], nos_def[elem.no₂.id])
        deform = (l₁ - l₀) / l₀
        forcas[i] = elem.area * elem.material.E * deform
    end

    return forcas
end

end

max_stars_repo_name: lorranfoliveira/ge

max_issues_repo_name: lorranfoliveira/ge

max_forks_repo_name: lorranfoliveira/ge


text:
const PointTypes = Union{VoronoiDelaunay.AbstractPoint2D, GeometryBasics.Point2}

function mean(pts::Vector{T}) where T <: PointTypes
	ax = 0.0
	ay = 0.0

	@simd for p in pts
		ax += getx(p)
		ay += gety(p)
	end

	Np = length(pts)
	T(ax/Np, ay/Np)
end

for name in [:sort, :sort!, :issorted, :sortperm]
	@eval begin
		function Base.$name(pts::Vector{T}) where T <: PointTypes
			center = mean(pts)
			$name(pts, by = p -> p - center, rev = true)
		end
	end
end

# http://stackoverflow.com/questions/6989100/sort-points-in-clockwise-order
function Base.isless(p::T, q::T) where T <: PointTypes
	if getx(p) >= 0.0 && getx(q) < 0.0
		return true
	elseif getx(p) < 0.0 && getx(q) >= 0.0
		return false
	elseif getx(p) == getx(q) == 0.0
		if gety(p) >= 0.0 || gety(q) >= 0.0
			return gety(p) > gety(q)
		else
			return gety(p) < gety(q)
		end
	end

	det = getx(p)*gety(q) - getx(q)*gety(p)
	if det < 0.0
		return true
	elseif det > 0.0
		return false
	end

	# p and q are on the same line from the center; check which one is
	# closer to the origin
	origin = T(0.0, 0.0)
	abs2(p, origin) > abs2(q, origin)
end


text:
import Base: extrema
using Optim

have_calculus = true
try
	Pkg.installed("Calculus")
catch
	have_calculus = false
end

if have_calculus
	import Calculus.derivative
end

create_tol = eps(Float64) # The max error tolerated when using a spline to approximate something else

function set_create_tol(newtol::Real)
	global create_tol = newtol
end

min_knot_dist_factor = 10000 #knots must be separated by at leas 10000*eps(T) to be considered unique.

function forwardDifference{T<:Number}(x::AbstractArray{T, 1})
    length = size(x, 1)
    return [x[i]-x[i-1] for i in 2:length]
end

type Spline{T<:Number}
    knots::Array{T, 1}
    values::Array{T, 1}
    a::Array{T, 1}
    b::Array{T, 1}
    c::Array{T, 1}
    d::Array{T, 1}

    function Spline(knots::AbstractArray{T, 1}, values::AbstractArray{T, 1})
	nKnots = length(knots)
    	if nKnots <= 1
		      throw(ArgumentError("There must be at least two knots"))
	end

	if length(values) != nKnots
		    throw(ArgumentError("Values is of wrong size"))
	end

	if any(isnan(knots))
		throw(ArgumentError("Knots cannot be NaN"))
	end

	if any(isnan(values))
		throw(ArgumentError("Values cannot be NaN"))
	end

        if !issorted(knots)
		throw(ArgumentError("Knots must be sorted"))
	end

	todelete = Array{Int, 1}(0)
	for knot in 2:nKnots
		if knots[knot] == knots[knot-1]
			throw(ArgumentError("Knots must all be unique"))
		end
		if knots[knot]-knots[knot-1] < min_knot_dist_factor*eps(T) #then replace the two knots & values by their averages
			knots[knot] = (knots[knot]+knots[knot-1])/2
			values[knot] = (values[knot]+values[knot-1])/2
			push!(todelete, knot-1)
		end
	end

	reverse!(todelete) #reverse the list of knots to delete so that removing one does not affect subsequent removals
	for i in todelete
	    deleteat!(knots, i)
	    deleteat!(values, i)
	end
	nKnots = length(knots)

	a = values
	h = forwardDifference(knots)
	df = forwardDifference(values) ./ h
	# Note: due the way that the gtsv! solver works, when solving A*X = B, this matrix can be X and then be overwritten By B. Therefore, the name of c is logical
	c = Array{Float64, 1}([0])
	append!(c, 3*forwardDifference(df))
	append!(c, [0])

	lowerDiag::Array{Float64, 1} = h[1:nKnots-2]
	append!(lowerDiag, [0.0]);

	diag = Array{Float64}([1])
	append!(diag, [2*(h[i]+h[i-1]) for i in 2:nKnots-1])
	append!(diag, [1.0])

	upperDiag = Array{Float64, 1}([0.0])
	append!(upperDiag, h[2:nKnots-1])

	#LAPACK.gtsv!([h[1:nKnots-2]; 0], [1; [2*(h[i]+h[i-1]) for i in 2:nKnots-1]; 1], [0; h[2:nKnots-1]], c;)
	LAPACK.gtsv!(lowerDiag, diag, upperDiag, c)
	d = [(c[i+1]-c[i])/(3*h[i]) for i in 1:nKnots-1]
	b = [df[i]-((h[i]/3)*((2*c[i])+c[i+1])) for i in 1:nKnots-1]
	new(knots, values, a[1:nKnots-1], b, c[1:nKnots-1], d)
    end
end

"""Recalculate a spline with new values, assuming that the corresponding knots are already in s.knots"""
function recalculate!{T}(s::Spline{T}, values::Array{T, 1})
  @assert issorted(s.knots) "Knots are not sorted. This is probably a bug in the Splines library"
	nKnots = length(s.knots)
	nvalues = length(values)

	if nvalues != nKnots
		throw(ArgumentError("Values is of wrong size"))
	end

	if any(isnan(values))
		throw(ArgumentError("values cannot contain NaN"))
	end

	h = forwardDifference(s.knots)
	todelete = Array{Int, 1}(0)
	if minimum(h) < min_knot_dist_factor*eps(T)
		for (knot, knotdiff) in enumerate(h)
			if knotdiff == zero(T)
				throw(ArgumentError("s.knots must all be unique"))
			end
			if knotdiff < min_knot_dist_factor*eps(T) #then replace the two s.knots & values by their averages
				s.knots[knot+1] = (s.knots[knot+1]+s.knots[knot])/2
				values[knot+1] = (values[knot+1]+values[knot])/2
				push!(todelete, knot)
			end
		end
	end
	reverse!(todelete) #reverse the list of knots to delete so that removing one does not affect subsequent removals

	for i in todelete
		deleteat!(s.knots, i)
		deleteat!(values, i)
	end
	nKnots = length(s.knots)

	a = values
	if length(todelete) > 0
		h = forwardDifference(s.knots)
	end
	df = forwardDifference(values) ./ h
	# Note: due the way that the gtsv! solver works, when solving A*X = B, this matrix can be X and then be overwritten By B. Therefore, the name of c is logical
	c = Array{Float64, 1}([0])
	append!(c, 3*forwardDifference(df))
	append!(c, [0])

	lowerDiag::Array{Float64, 1} = h[1:nKnots-2]
	append!(lowerDiag, [0.0]);

	diag = Array{Float64}([1])
	append!(diag, [2*(h[i]+h[i-1]) for i in 2:nKnots-1])
	append!(diag, [1.0])

	upperDiag = Array{Float64, 1}([0.0])
	append!(upperDiag, h[2:nKnots-1])

	#LAPACK.gtsv!([h[1:nKnots-2]; 0], [1; [2*(h[i]+h[i-1]) for i in 2:nKnots-1]; 1], [0; h[2:nKnots-1]], c;)
	LAPACK.gtsv!(lowerDiag, diag, upperDiag, c)
	d = [(c[i+1]-c[i])/(3*h[i]) for i in 1:nKnots-1]
	b = [df[i]-((h[i]/3)*((2*c[i])+c[i+1])) for i in 1:nKnots-1]
	s.a = a[1:nKnots-1]
	s.b = b
	s.c = c[1:nKnots-1]
	s.d = d
	s.values = values
#	new(knots, a[1:nKnots-1], b, c[1:nKnots-1], d)
	return nothing
end

function indomain{T}(s::Spline{T}, x::T)
	return (x>=s.knots[1] && x<=s.knots[end])
end

function domain{T}(s::Spline{T})
	return (xmin(s), xmax(s))
end

Spline{T<:Number}(knots::AbstractArray{T, 1}, values::AbstractArray{T, 1}) = Spline{T}(knots, values)

#A helper function to find the max error of s(t)-f(t) for a spline s and a callable object f. If f is not continuous and differentiable, this
#function WILL NOT WORK.
function find_maxerr{T}(s::Spline{T}, f, x_lower::Real=domain(s)[1], x_upper::Real=domain(s)[2])
	min = optimize(t->1/((f(t)-s(t))^2+1), x_lower, x_upper)

	return min.minimum
end

# Tries to construct a spline for anything that has a method of call() defined for it. If f is not smooth and continuous, this will not yield good results
function Spline(f, xlower::Real, xupper::Real, tol=create_tol)
	#Start with a really simple set of knots
	knots = [xlower, xlower+xupper/2, xupper]
	values = Float64[f(i) for i in knots]
	s = Spline(knots, values)

	maxerr = find_maxerr(s, f, xlower, xupper)
	err_at_maxerr = abs(f(maxerr)-s(maxerr))

	loopcount = 0
	while err_at_maxerr > tol
		newknot = maxerr
		newvalue = f(newknot)
		insert!(s, newknot, newvalue)

		maxerr = find_maxerr(s, f, xlower, xupper)
		err_at_maxerr = abs(f(maxerr)-s(maxerr))

		loopcount += 1

		if loopcount > 10000
			warn("Creating a spline from a callable object failed after 10000 iterations")
			break
		end
	end

return s
end

function (s::Spline{T}){T}(x::T)
	if x == domain(s)[2]
		return s.values[end]
	end

	if !indomain(s, x)
		throw(DomainError())
		#error("x out of range")
	end

	bracket = searchsortedlast(s.knots, x)
	if bracket == length(s.knots)
		bracket -= 1
	end

	h = x-s.knots[bracket]

	return @evalpoly(h, s.a[bracket], s.b[bracket], s.c[bracket], s.d[bracket])::T
end

"""calling a spline _f_ with another spline _g_ as argument produces a spline approximation to f(g)"""
function (f::Spline{T}){T}(g::Spline{T})
  gmin, gmax = extrema(g)
  if !indomain(f, gmax) || !indomain(f, gmin)
    warn("f(g) is not defined; Range of g is [$gmin, $gmax], but domain of f is only [$(domain(f)[1]), $(domain(f)[2])]. Throwing DomainError")
    throw(DomainError())
  end

  # approixmating f(g) is best done by first approximating g with first-order polynomials, and then strechting or shrinking intervals of f appropriately
  g_discrete = adaptive_discretize(g)

  #However, it is very necessary that the knots/values in f be preserved in the transformation;
  #otherwise, the transformed spline might bear little resemblance to what it should.
  #The simplest way to do this is, for each knot $f_k$ in f to put a knot $g_x$ into g's list of knots such that g($g_k$) = $f_k$

  newknots = g_discrete[:x]
  newvalues = [f(x) for x in g_discrete[:y]]
  f_warped = Spline(newknots, newvalues)

  #We now have some approximation to f(g(t)), but it probably isn't very good.
  # So we begin going through repeatedly finding the location of the max error
  # and adding a knot there
  tol = .001 #FIXME: This is purely arbitrary and probably a bad value
  maxϵ_pos = find_maxerr(f_warped, t->f(g(t)))
  maxϵ = abs(f_warped(maxϵ_pos)-f(g(maxϵ_pos)))
  while maxϵ > tol
    insert!(f_warped, maxϵ_pos, f(g(maxϵ_pos)))
    maxϵ_pos = find_maxerr(f_warped, t->f(g(t)))
    maxϵ = abs(f_warped(maxϵ_pos)-f(g(maxϵ_pos)))
  end

  return f_warped
end

"""Add a knot to a spline"""
function insert!{T}(s::Spline{T}, knot::T, value::T)
	knots = s.knots
	values = [s(i)::T for i in knots] #FIXME: why allocate a whole new array? all but one of the values are already stored in s.a
	if knot in knots
		values[findin(knots, knot)[1]] == value #FIXME: this should be assignment
		recalculate!(s, values)
		return nothing
	end

	push!(knots, knot)
	push!(values, value)

	perm = sortperm(knots)
	s.knots = knots[perm]
	values = values[perm]

	recalculate!(s, values)
	return nothing
end

#"""Add a bunch of knots to a spline all at once"""
function insert!{T}(s::Spline{T}, knots::AbstractArray{T, 1}, values::AbstractArray{T, 1})
  if length(knots) != length(values)
    throw(ArgumentError("knots and values must be the same size"))
  end
  if unique(knots) != knots
    throw(ArgumentError("knots must all be unique"))
  end

  newvalues = copy(s.values)

  #If there are knots in _knots_ that are already in the spline ,then just change the value and add
  #that index to a list of knots and values to be deleted
  todelete = Array{Int, 1}(0)
  for (i, knot) in enumerate(knots)
    if knot in s.knots
      duplicate_pos = findin(s.knots, knot)[1]
      newvalues[duplicate_pos] = values[i] #change the value
      push!(todelete, i) #and add the index to a list of indecies to delete
    end
  end
  #delete whatever was in the list
  reverse!(todelete) #we have to delete starting at the end so that a deletion at one index does not throw off future deletions
  for i in todelete
    deleteat!(knots, i)
    deleteat!(values, i)
  end

  append!(s.knots, knots)
  append!(newvalues, values)

  #Now the knots must be sorted to be in increasing order and the values must re-arranged in the same manner
  perm = sortperm(s.knots)
  s.knots = s.knots[perm]
  newvalues = newvalues[perm]

  recalculate!(s, newvalues)
  return nothing
end

"""Move the knot at knotIndex to the location newlocation"""
function moveknot!{T}(s::Spline{T}, knotIndex::Int, newlocation::T)
  if knotIndex > length(s.knots) || knotIndex <1
    throw(BoundsError(s, knotIndex))
  end

  if !indomain(s, newlocation)
    throw(DomainError())
  end

  if newlocation in s.knots
    throw(ArgumentError("newlocation cannot be a pre-existing knot"))
  end

  s.knots[knotIndex] = newlocation

  recalculate!(s, s.values)
end

"""Move the knots at _knotIndecies_ to the new locations in _newlocations_"""
function moveknot!{T}(s::Spline{T}, knotIndecies::AbstractArray{Int, 1}, newlocations::AbstractArray{T, 1})
  if length(knotIndecies) != length(newlocations)
    throw(ArgumentError("knotIndecies and newlocations must have the same length, not $(length(knotIndecies)) and $(length(newlocations))"))
  end

  for index in knotIndecies
    if index > length(s.knots) || index < 1
      throw(BoundsError(s, index))
    end
  end

  for newloc in newlocations
    if !indomain(s, newloc)
      throw(DomainError)
    end
  end

  for (i, index) in enumerate(knotIndecies)
    s.knots[index] = newlocations[i]
  end

  recalculate!(s, s.values)
end

"""Delete a knot or knots from a spline"""
function deleteknotat!{T}(s::Spline{T}, knotIndex::Integer...)
	knots = s.knots
	values = Array{T, 1}(length(knots))
	for (i, x) in enumerate(knots)
		values[i] = s(x)
	end

	deleteat!(knots, knotIndex)
	deleteat!(values, knotIndex)

	recalculate!(s, values)
end

"""Trim the beginning (everything less than cutoff) off of a spline"""
function trimstart!{T}(s::Spline{T}, cutoff::T)
	if !indomain(s, cutoff)
		throw(DomainError())
	end

	newknot = cutoff
	newvalue = s(newknot)
	insert!(s, newknot, newvalue)

	indeciesToTrim = find(x->x < newknot, s.knots)
	deleteknotat!(s, indeciesToTrim...)
end

"""Trim the end (everything greater than cutoff) off of a spline"""
function trimend!{T}(s::Spline{T}, cutoff::T)
	if !indomain(s, cutoff)
		throw(DomainError())
	end

	newknot = cutoff
	newvalue = s(newknot)
	insert!(s, newknot, newvalue)

	indeciesToTrim = find(x->x > newknot, s.knots)
	deleteknotat!(s, indeciesToTrim...)
end

"""Trim the start and the end (everything less than startx or greater than endx) off of a spline"""
function trim!{T}(s::Spline{T}, startx::T, endx::T)
	insert!(s, startx, s(startx))
	insert!(s, endx, s(endx))
	indeciesToTrim = find(x->x<startx || x>endx, s.knots)
	deleteknotat!(s, indeciesToTrim...)
end

include("cubic_helpers.jl")

function maximum{T}(s::Spline{T})
	maxS = -Inf

	for i in 1:length(s.knots)-1
		maxS = max(maxS, maxcubic(s.knots[i+1]-s.knots[i], s.a[i], s.b[i], s.c[i], s.d[i])[1])
	end

	return maxS
end

function minimum{T}(s::Spline{T})
	minS = typemax(T)

	for i in 1:length(s.knots)-1
		minS = min(minS, mincubic(s.knots[i+1]-s.knots[i], s.a[i], s.b[i], s.c[i], s.d[i])[1])
	end

	return minS
end

function extrema{T}(s::Spline{T})
	minS = typemax(T)
	maxS = typemin(T)

	for i in 1:length(s.knots)-1
		minS = min(minS, mincubic(s.knots[i+1]-s.knots[i], s.a[i], s.b[i], s.c[i], s.d[i])[1])
		maxS = max(maxS, maxcubic(s.knots[i+1]-s.knots[i], s.a[i], s.b[i], s.c[i], s.d[i])[1])
	end

	return (minS, maxS)
end

"""Similar to indmax for arrays, return the x-position of the maximum of the given spline"""
function xmax{T}(s::Spline{T})
	theMax = maximum(s)
	maxPos = 0
	curMax = (0, 0)
	for i in 1:length(s.knots)-1
		curMax = maxcubic(s.knots[i+1]-s.knots[i], s.a[i], s.b[i], s.c[i], s.d[i])
		if curMax[1] == theMax
			maxPos = s.knots[i] + curMax[2]
			break
		end
	end

	return maxPos
end

function maxabs(s::Spline)
	return max(abs(maximum(s)), abs(minimum(s)))
end

"""Find the appropriate step size for discretizing a spline with a uniform mesh maximum absolute error _tol_ due to linear interpolation. If negative, this will be interpreted as a relative error from the maximum value"""
function find_stepsize{T}(s::Spline{T}, tol::Number)
	if tol < 0
		tol = maximum(s)*abs(tol)
	end
	minValues = 100
	maxa2d = maxabs2deriv(s)
	if maxa2d > 0
		stepSize = .9*min(sqrt(8*tol/maxa2d), (s.knots[end]-s.knots[1])/minValues)
	else
		stepSize = (s.knots[end]-s.knots[1])/minValues
	end

	return stepSize
end

"""
Create an evenly spaced mesh of points to discretize the spline so that the maximum absolute error due to linear interpolation is <= tol.
If this is too small, return a mesh with 100 points. Note, this may use horrendous amounts of memory. Consider using adaptive_discretize_mesh instead
"""
function uniform_discretize_mesh{T}(s::Spline{T}, tol::Number=-1)
	stepSize = find_stepsize(s, tol)
	mesh = s.knots[1]:stepSize:s.knots[end]
	return mesh
end

"""Turn a spline object into a discrete list of values on an automatically created uniform mesh"""
function uniform_discretize{T}(s::Spline{T}; tol::Number=-1) #default .001 absolute error tolerated
	mesh = uniform_discretize_mesh(s, tol)

	return discretize(s, mesh)
end

"""Turn a spline object _s_ into a discrete list of values, seperated on the domain of s by _step_"""
function uniform_discretize{T}(s::Spline{T}, step::T)
	mesh = xmin(s):step:xmax(s)
	return discretize(s, mesh)
end

"""Turn a spline object into a discrete list of values at the mesh points given by _mesh_"""
function discretize{T}(s::Spline{T}, mesh::AbstractArray{T})
	discretized = Array{T, 1}(length(mesh))
	for (i, x) in enumerate(mesh)
		discretized[i] = s(x)
	end

	return discretized
end

"""A helper function to find the maximum absolute error due to linear interpolation of
a cubic spline over a mesh _x_

Parameters:
 - *s* is the spline that is being discretized
 - *x* holds x-values that need to be used between the current knot (which is not easily visible to this function) and the one after it to get a reasonable approximation.
 - *a*, *b*, *c*, and *d* are the polynomial coefficients of the spline at the knot where this is being done.
"""
function find_maxϵ{T}(s::Spline{T}, x::Array{T, 1}, a::T, b::T, c::T, d::T)
  if any(isnan(x))
    throw(ArgumentError("NAN encountered in find_maxϵ"))
  end

  maxϵ = maxϵ_pos = maxϵ_index = typemin(T) #maxcubic(x[1], x[2], 0, s.b[i]-((s(s.knots[i]+x[2])-s(s.knots[i]+x[1]))/(x[2]-x[1])), s.c[i], s.d[i])
	for (j, xj) in take(enumerate(x), length(x)-1)
    lin_approx_slope = (s(x[j+1])-s(xj))/(x[j+1]-xj)
	  maxϵ_new, maxϵ_pos_new = maxabscubic(xj-x[1], x[j+1]-x[1], a-(s(xj)-lin_approx_slope*(xj-x[1])), b-lin_approx_slope, c, d)
    @assert maxϵ_pos != xj-x[1] "$xj $(a-(s(xj)-lin_approx_slope*(xj-x[1])))"

	  if maxϵ_new > maxϵ
	    maxϵ = maxϵ_new
		  maxϵ_pos = maxϵ_pos_new
		  maxϵ_index = j
    end
  end

  return maxϵ, maxϵ_pos, maxϵ_index
end

"""
A helper function for adaptive_discretize_mesh

Parameters:
 - *s* is the spline that is being discretized.
 - *x* holds x-values that need to be used between the current knot (which is not easily visible to this function) and the one after it to get a reasonable approximation.
 - *absTol* is the maximum absolute error that is tolerated
 - *a*, *b*, *c*, and *d* are the polynomial coefficients of the spline at the knot where this is being done.

Return:
 - The new error after adding zero or more values to the mesh *x*

Given _s_ and _x_, this function finds the maximum error associated with using linear
interpolation over the mesh _x_. If this error is too large, then it adds a new value
to _x_ where the max error was and returns the new error.
"""
function densify_mesh!{T}(s::Spline{T}, x::Array{T, 1}, absTol::T, a::T, b::T, c::T, d::T)
  maxϵ, maxϵ_pos, maxϵ_index = find_maxϵ(s, x, a, b, c, d)
  @assert !isnan(maxϵ_pos)
  if maxϵ > absTol
    if maxϵ_pos == x[maxϵ_index+1]
      insert!(x, maxϵ_index+1, (x[maxϵ_index+1]+x[maxϵ_index])/2)
    elseif maxϵ_pos in x
      error("duplicate value $maxϵ_pos x=$x")
    else
      insert!(x, maxϵ_index+1, x[1]+maxϵ_pos)
    end
    #but maxϵ has now changed, so we re-compute it
    maxϵ = find_maxϵ(s, x, a, b, c, d)[1]
  end

  return maxϵ
end

"""
Create an adaptive mesh of points for discretizing a spline _s_.

Parameters:
 - *s* A spline to discretize
 - *tol* The tolerance. If this is negative, then it is assumed to be a relative tolerance (with a few modifications to get around the problem of the spline being 0 somewhere)
 otherwise, it is assumed to be an absolute tolerance.
"""
function adaptive_discretize_mesh{T}(s::Spline{T}, tol::Number=-.01)

	if tol > 0
		absTol = tol
	end
	nKnots = length(s.knots)

	x = Array{Array{T, 1}, 1}(nKnots)

	for i in 1:nKnots-1
		x[i] = [s.knots[i], s.knots[i+1]]
		maxϵ = maxabscubic(0.0, x[i][2]-x[i][1], zero(T), s.b[i]-((s.values[i+1]-s.values[i])/(s.knots[i+1]-s.knots[i])), s.c[i], s.d[i])[1]

		if tol < 0
			absTol = abs(tol)*maxabscubic(zero(T), s.knots[i+1]-s.knots[i], s.a[i], s.b[i], s.c[i], s.d[i])[1]
		end

		while maxϵ > absTol
      maxϵ = densify_mesh!(s, x[i], absTol, s.a[i], s.b[i], s.c[i], s.d[i])
		end
		deleteat!(x[i], length(x[i]))
	end
	x[nKnots] = [s.knots[nKnots]]

	return vcat(x...)::Array{T, 1}
end

function adaptive_discretize{T}(s::Spline{T}, tol::Number=-.01)
	mesh = adaptive_discretize_mesh(s, tol)
	y = [s(x) for x in mesh]

	return Dict(:x => mesh, :y => y)
end

"""Scale a spline so that its maximum absolute value is __normalize_to__"""
function normalize_max!{T}(s::Spline{T}, normalize_to::T=one(T))
	s2 = s/maxabs(s)
	s.knots = s2.knots
	s.values = s2.values
	s.a = s2.a
	s.b = s2.b
	s.c = s2.c
	s.d = s2.d
	return nothing
end

"""Return the smallest value of x for which s(x) is defined"""
function xmin{T}(s::Spline{T})
	return s.knots[1]
end

"""Return the largest value of x for which s(x) is defined"""
function xmax{T}(s::Spline{T})
	return s.knots[end]
end

"""Shift the domain of a spline by amount _Δx_"""
function domain_shift!{T}(s::Spline{T}, Δx::T)
	for i in 1:length(s.knots)
		s.knots[i] += Δx
	end
end

export call
export insert!
export moveknot!
export deleteknotat!
export trimstart!
export trimend!
export trim!
export maximum
export minimum
export extrema
export xmax
export maxabs
export Spline
export discretize
export uniform_discretize
export adaptive_discretize
export indomain
export domain
export normalize_max!
export xmin
export xmax
export domain_shift!


text:
module EMLtypes
	using DanaTypes
	using DotPlusInheritance
	using Reexport
	include("types/coefficient.jl")
	include("types/constant.jl")
	include("types/positive.jl")
	include("types/negative.jl")
	include("types/fraction.jl")
	include("types/percent.jl")
	include("types/control_signal.jl")
	include("types/efficiency.jl")
	include("types/pressure.jl")
	include("types/press_delta.jl")
	include("types/head_mass.jl")
	include("types/head.jl")
	include("types/temperature.jl")
	include("types/temp_delta.jl")
	include("types/time_h.jl")
	include("types/time_min.jl")
	include("types/time_sec.jl")
	include("types/frequency.jl")
	include("types/angle.jl")
	include("types/area.jl")
	include("types/length.jl")
	include("types/length_delta.jl")
	include("types/volume.jl")
	include("types/volume_mol.jl")
	include("types/volume_mass.jl")
	include("types/current.jl")
	include("types/charge.jl")
	include("types/capacitance.jl")
	include("types/indutance.jl")
	include("types/voltage.jl")
	include("types/resistance.jl")
	include("types/potency.jl")
	include("types/currency.jl")
	include("types/mass.jl")
	include("types/mol.jl")
	include("types/kmol.jl")
	include("types/molweight.jl")
	include("types/molweight_inv.jl")
	include("types/dens_mol.jl")
	include("types/dens_mass.jl")
	include("types/conc_mol.jl")
	include("types/inv_conc_mol.jl")
	include("types/conc_mass.jl")
	include("types/inv_conc_mass.jl")
	include("types/reaction_mol.jl")
	include("types/reaction_mass.jl")
	include("types/cp_mass.jl")
	include("types/cp_mol.jl")
	include("types/cv_mol.jl")
	include("types/enth_mass.jl")
	include("types/enth_mol.jl")
	include("types/entr_mol.jl")
	include("types/entr_mass.jl")
	include("types/heat_reaction.jl")
	include("types/heat_rate.jl")
	include("types/heat_flux.jl")
	include("types/heat_trans_coeff.jl")
	include("types/energy.jl")
	include("types/energy_mass.jl")
	include("types/energy_mol.jl")
	include("types/power.jl")
	include("types/flow_mass.jl")
	include("types/flow_mass_delta.jl")
	include("types/flow_vol.jl")
	include("types/flow_vol_delta.jl")
	include("types/flow_mol.jl")
	include("types/flow_mol_delta.jl")
	include("types/flux_mol.jl")
	include("types/flux_mol_delta.jl")
	include("types/flux_mass.jl")
	include("types/flux_mass_delta.jl")
	include("types/flux_vol.jl")
	include("types/flux_vol_delta.jl")
	include("types/vel_angular.jl")
	include("types/rotation.jl")
	include("types/velocity.jl")
	include("types/velocity_delta.jl")
	include("types/acceleration.jl")
	include("types/fricfactor.jl")
	include("types/moment_inertia.jl")
	include("types/hookes_const.jl")
	include("types/conductivity.jl")
	include("types/diffusivity.jl")
	include("types/fugacity.jl")
	include("types/viscosity.jl")
	include("types/vol_mol.jl")
	include("types/vol_mass.jl")
	include("types/spec_surface_vol.jl")
	include("types/spec_surface_mass.jl")
	include("types/surf_tens.jl")
	include("types/act_coeff.jl")
	include("types/ph.jl")
end


text:
"""
This is the IRNN model, a recurrent net with relu activations whose
recurrent weights are initialized with the identity matrix.  From: Le,
Q. V., Jaitly, N., & Hinton, G. E. (2015). A Simple Way to Initialize
Recurrent Networks of Rectified Linear Units. arXiv preprint
arXiv:1504.00941.
"""

irnn(;n=1,std=0.01) = quote
    x1 = input()
    w1 = par($n,0; init=Gaussian(0,$std))
    x2 = dot(w1,x1)
    w6 = par($n,0; init=Identity())
    x3 = dot(w6,x6)
    x4 = add(x2,x3)
    w4 = par(0; init=Constant(0))
    x5 = add(w4,x4)
    x6 = relu(x5)
end

# IRNN(n; std=0.01)=Net(Mmul(n; init=randn!, initp=(0,std)),
#                       (Mmul(n; init=eye!), 5),
#                       Add2(), Bias(), Relu())

# eye!(a)=copy!(a, eye(eltype(a), size(a)...)) #


text:
function P(eos::AbstractEquationOfState, V, T) end
function A(eos::AbstractEquationOfState, ρ, T) end


text:
using DataFrames, CSV
#gr()
#pgfplots()
a = 2
#pyplot(guidefont=font(a*17),tickfont=font(a*15),legendfont=font(a*12),titlefont=font(a*20))
#a=0.25;
#default(guidefont=font(a*17),tickfont=font(a*15),legendfont=font(a*12),titlefont=font(a*20))
plotSettings(;(:polyPts=>false),(:simulate=>true),(:mpc_markers =>(:circle,:blueviolet,0.0,0.0)),(:plant=>true),(:plantOnly=>false),(:size=>(a*900,a*600)),(:format=>"png"));

resultsDir!(n)
dfs = case2dfs(c)
filename = "case"
CSV.write(string("results/",filename,".csv"), dfs; quotechar = ' ')

################################################
# save obstacle and vehicle data for ploting
pts = 4 # sample at several points in time
t = Vector(linspace(0,n.r.ip.plant[:t][end],pts)) # assuming it starts at 0!
#t[3] = 19.5 # for s9_D2 to line up crash
indmin(abs.(t[2]-n.r.ip.plant[:t]))

xv = zeros(pts)
yv = zeros(pts)
psi = zeros(pts)
dfs = DataFrame()
for j in 1:pts
  idx = indmin(abs.(t[j]-n.r.ip.plant[:t]))
  xv[j] = n.r.ip.plant[:x][idx]
  yv[j] = n.r.ip.plant[:y][idx]
  psi[j] = n.r.ip.plant[:psi][idx]
end
dfs[:idx] = 1:pts
dfs[:t] = round.(t,1,10) # round for graph
dfs[:xv] = xv
dfs[:yv] = yv
dfs[:psi] = psi
cd(n.r.resultsDir)
  CSV.write("veh.csv", dfs; quotechar = ' ')
cd(n.r.mainDir)

for i in 1:length(c["obstacle"]["x0"])
  xo = zeros(pts)
  yo = zeros(pts)
  ro = zeros(pts)
  clr = zeros(pts)
  dfs = DataFrame()
  for j in 1:pts
    xo[j] = c["obstacle"]["x0"][i] + c["obstacle"]["vx"][i]*t[j]
    yo[j] = c["obstacle"]["y0"][i] + c["obstacle"]["vy"][i]*t[j]
    ro[j] = c["obstacle"]["radius"][i]
    clr[j] = (j)/pts*100
  end
  dfs[:idx] = 1:pts
  dfs[:t] = round.(t,1,10) # round for graph
  dfs[:xo] = xo
  dfs[:yo] = yo
  dfs[:ro] = ro
  dfs[:xv] = xv
  dfs[:yv] = yv
  dfs[:clr] = clr
  cd(n.r.resultsDir)
    CSV.write(string("obs_",i,".csv"), dfs; quotechar = ' ')
  cd(n.r.mainDir)
end
# save obstacle and vehicle data for ploting
################################################

savePlantData!(n)
saveOptData(n)

if _pretty_defaults[:simulate]
  println("Plotting the Final Results!")
  mainSim(n;(:mode=>:intersection))
end

#optPlot(n)
#posterP(n)

# calculate the control effort
dfs = DataFrame()
n.r.ip.plant[:sa][isnan.(n.r.ip.plant[:sa])]=0 # remove any NaNs
n.r.ip.plant[:sr][isnan.(n.r.ip.plant[:sr])]=0
n.r.ip.plant[:jx][isnan.(n.r.ip.plant[:jx])]=0
dfs[:sa] = mean(n.r.ip.plant[:sa].^2)
dfs[:sr] = mean(n.r.ip.plant[:sr].^2)
dfs[:jx] = mean(n.r.ip.plant[:jx].^2)
dfs[:total] = dfs[:sa] + dfs[:sr] + dfs[:jx]
cd(n.r.resultsDir)
  CSV.write("ce.csv", dfs; quotechar = ' ')
cd(n.r.mainDir)


text:
module Groebner

debug() = false

PRIMES = 0
F4TIME = 0.0
RECTIME = 0.0
CORRTIME = 0.0

ADD_ROW_HIT  = 0
ADD_ROW_MISS = 0
ADD_ROW_COLLISION = 0

function printhash()
    global ADD_ROW_HIT
    global ADD_ROW_MISS
    global ADD_ROW_COLLISION

    println("ADD_ROW_HIT = $(ADD_ROW_HIT)")
    println("ADD_ROW_MISS = $(ADD_ROW_MISS)")
    println("ADD_ROW_COLLISION = $(ADD_ROW_COLLISION)")

    ADD_ROW_HIT  = 0
    ADD_ROW_MISS = 0
    ADD_ROW_COLLISION = 0
end
ADD_ROW_HIT  = 0
ADD_ROW_MISS = 0

function printall()
    global PRIMES
    global F4TIME
    global RECTIME
    global CORRTIME
    s = sum((F4TIME, RECTIME, CORRTIME))
    println("run success,")
    println("PRIMES = $PRIMES")
    println("F4TIME = $(F4TIME/s)")
    println("RECTIME = $(RECTIME/s)")
    println("CORRTIME = $(CORRTIME/s)")
    PRIMES = 0
    F4TIME = 0.0
    RECTIME = 0.0
    CORRTIME = 0.0
end

# most of these are not needed
# TODO
import AbstractAlgebra
import AbstractAlgebra.Generic: MPoly, GFElem, MPolyRing, Poly
import AbstractAlgebra: leading_term, QQ, PolynomialRing, terms,
                        coeff, divides, base_ring, elem_type,
                        rref, isconstant, leading_coefficient,
                        map_coefficients, monomials, degree,
                        degrees, isconstant, leading_monomial,
                        GF, gens, MatrixSpace, coefficients,
                        crt, ordering, exponent_vectors, lift,
                        MPolyBuildCtx, finish, push_term!, ZZ,
                        content, change_base_ring, exponent_vector,
                        lcm, monomial, RingElem, set_exponent_vector!,
                        nvars, data, characteristic, isdivisible_by,
                        divexact, symbols

# for example systems
import Combinatorics

import Primes
import Primes: nextprime

import Random

import Logging
import Logging: ConsoleLogger, LogLevel

import MultivariatePolynomials
import MultivariatePolynomials: AbstractPolynomial, AbstractPolynomialLike

# type aliases for internal objects
include("internaltypes.jl")

# so, what is that exactly?
include("utils.jl")

# some simple reference implementations
include("common.jl")

# computation parameters control
include("metainfo.jl")

# input-output conversions for polynomials
include("io.jl")

#= operations with numbers =#
# unsigned arithmetic for finite field computations
include("arithmetic/unsigned.jl")
# modular arithmetic for CRT and rational reconstruction
include("arithmetic/modular.jl")

#= generic f4 implementation =#
# the heart of this library
include("f4/structs.jl")
include("f4/symbolic.jl")
include("f4/hash.jl")
include("f4/linear.jl")
include("f4/sorting.jl")
include("f4/lucky.jl")
include("f4/coeffs.jl")
include("f4/f4.jl")

# f4 with tracing
include("f4/f4trace.jl")

include("f4/groebner.jl")
include("f4/isgroebner.jl")
include("f4/normalform.jl")
include("f4/correctness.jl")
# include("f4/statistics.jl")

#= generic fglm implementation =#
include("fglm/linear.jl")
include("fglm/fglm.jl")

# api
include("interface.jl")

# example systems definitions
include("testgens.jl")

export groebner
export isgroebner
export normalform

export fglm
export kbase

end


text: println("Hello Wolrd");


text:
"""
    mprism(x, y, z, xextent, yextent, zextent, M)

Compute the magnetic anomaly of a prism bounded by the axis-parallel planes defined
by `xextent`, `yextent`, and `zextent`.
Observation point is given by cartesian coordinates `x, y, z`.
The magnetization of the body is defined by the three components of `M` measured in units of A/m.
"""
function mprism(x::T, y::T, z::T, x1::Array{T}, y1::Array{T}, z1::Array{T}, M::Array{T}) where T<:Real
    dZ = zero(eltype(T))
    dX = zero(eltype(T))
    dY = zero(eltype(T))
    f = zero(eltype(T))
    g12 = 0.0
    g13 = 0.0
    g23 = 0.0
    for i in 1:2
        for j in 1:2
            for k in 1:2
                g12 = g1(x - x1[i], z - z1[k], y - y1[j])
                g13 = g1(x - x1[i], y - y1[j], z - z1[k])
                g23 = g1(y - y1[j], x - x1[i], z - z1[k])

                f = (-1.0)^(i+j+k)
                dX += f * (
                M[1] * g2(y - y1[j], z - z1[k], x - x1[i]) +
                M[2] * g12 + #g1(x - x1[i], z - z1[k], y - y1[j])
                M[3] * g13 #g1(x - x1[i], y - y1[j], z - z1[k])
                )

                dY += f * (
                M[1] * g12 + #g1(x - x1[i], z - z1[k], y - y1[j])
                M[2] * g2(x - x1[i], z - z1[k], y - y1[j]) +
                M[3] * g23 #g1(y - y1[j], x - x1[i], z - z1[k])
                )

                dZ += f * (
                M[1] * g13 + #g1(x - x1[i], y - y1[j], z - z1[k]) +
                M[2] * g23 + #g1(y - y1[j], x - x1[i], z - z1[k]) +
                M[3] * g2(x - x1[i], y - y1[j], z - z1[k]))
            end
        end
    end
    dX *= 1e-7
    dY *= 1e-7
    dZ *= 1e-7
    return [dX, dY, dZ]
end

function g1(u::T, v::T, w::T) where T<:Real
    return sign(v) * log(
        sqrt(u^2 + w^2) /
        (abs(v) + sqrt(u^2 + v^2 + w^2))
        )
end

function g2(u::T, v::T, w::T) where T<:Real
    return atan(u * v / (w * sqrt(u^2 + v^2 + w^2)))
end

"""
    mprism_(x, y, z, xextent, yextent, zextent, M)

Compute the magnetic anomaly B of a prism bounded by the axis-parallel planes defined
by `xextent`, `yextent`, and `zextent`.
Observation point is given by cartesian coordinates `x, y, z`.
The magnetization of the body is defined by the three components of `M` measured in units of A/m.
"""
function mprism_(x::T, y::T, z::T, x1::Array{T}, y1::Array{T}, z1::Array{T}, M::Array{T}) where T<:Real
	G11 = G2(y, z, x, y1, z1, x1)
	G12 = G1(x, z, y, x1, z1, y1)
	G13 = G1(x, y, z, x1, y1, z1)
	#G21 = G12
	G22 = G2(x, z, y, x1, z1, y1)
	G23 = G1(y, x, z, y1, x1, z1)
	#G31 = G13
	#G32 = G23
	G33 = G2(x, y, z, x1, y1, z1)
	G = Array{T}(undef, 3, 3)
	G .= [G11 G12 G13
		G12 G22 G23
		G13 G23 G33]

	B = G * M
end

"""
    mprism!(B, x, y, z, xextent, yextent, zextent, M)

Compute the magnetic anomaly B of a prism bounded by the axis-parallel planes defined
by `xextent`, `yextent`, and `zextent`.
Observation point is given by cartesian coordinates `x, y, z`.
The magnetization of the body is defined by the three components of `M` measured in units of A/m.

Note that this function overrides the value of the argument `B`.
"""
function mprism!(B::Array{T}, x::T, y::T, z::T, x1::Array{T}, y1::Array{T}, z1::Array{T}, M::Array{T}) where T<:Real
	G11 = G2(y, z, x, y1, z1, x1)
	G12 = G1(x, z, y, x1, z1, y1)
	G13 = G1(x, y, z, x1, y1, z1)
	#G21 = G12
	G22 = G2(x, z, y, x1, z1, y1)
	G23 = G1(y, x, z, y1, x1, z1)
	#G31 = G13
	#G32 = G23
	G33 = G2(x, y, z, x1, y1, z1)

	B .= [G11 G12 G13
		G12 G22 G23
		G13 G23 G33] * M
	return nothing
	# B .= G * M
end

function G2(x::T, y::T, z::T, x1::Array{T, 1}, y1::Array{T, 1}, z1::Array{T, 1}) where T<:Real

        # arg1 = (x - x1[1], y - y1[1], z - z1[1])
        # arg2 = (x - x1[2], y - y1[2], z - z1[2])

        return -g2(x - x1[1], y - y1[1], z - z1[1]) +
         g2(x - x1[1], y - y1[1], z - z1[2]) +
         g2(x - x1[1], y - y1[2], z - z1[1]) -
         g2(x - x1[1], y - y1[2], z - z1[2]) +
         g2(x - x1[2], y - y1[1], z - z1[1]) -
         g2(x - x1[2], y - y1[1], z - z1[2]) -
         g2(x - x1[2], y - y1[2], z - z1[1]) +
         g2(x - x1[2], y - y1[2], z - z1[2])
end

function G1(x::T, y::T, z::T, x1::Array{T, 1}, y1::Array{T, 1}, z1::Array{T, 1}) where T<:Real
        return -g1(x - x1[1], y - y1[1], z - z1[1]) +
         g1(x - x1[1], y - y1[1], z - z1[2]) +
         g1(x - x1[1], y - y1[2], z - z1[1]) -
         g1(x - x1[1], y - y1[2], z - z1[2]) +
         g1(x - x1[2], y - y1[1], z - z1[1]) -
         g1(x - x1[2], y - y1[1], z - z1[2]) -
         g1(x - x1[2], y - y1[2], z - z1[1]) +
         g1(x - x1[2], y - y1[2], z - z1[2])
end

max_stars_repo_name: ruboerner/GravMag

max_issues_repo_name: ruboerner/GravMag

max_forks_repo_name: ruboerner/GravMag


text:
N = 100_000

fit(Bernoulli, rand(Bernoulli(0.7), N))

fit(Beta, rand(Beta(1.3, 3.7), N))

fit(Binomial, rand(Binomial(N, 0.3)), N)

fit(DiscreteUniform, rand(DiscreteUniform(300_000, 700_000), N))

fit(Exponential, rand(Exponential(0.1), N))

fit(Gamma, rand(Gamma(7.9, 3.1), N))

fit(Geometric, rand(Geometric(0.1), N))

fit(Laplace, rand(Laplace(10.0, 3.0), N))

fit(Normal, rand(Normal(11.3, 5.7), N))

fit(Poisson, rand(Poisson(19.0), N))

fit(Uniform, rand(Uniform(1.1, 98.3), N))


text:
export RFFTSpectrumArray, resamplefreq, ntimeframes

# ----- RFFTSpectrumArray ---------------------------
struct RFFTSpectrumArray{T} <: AbstractSpectrumArray{T}
    data::Matrix{T} # non-interleaving frequencies x channels
    rate::Float64 # sampling frequency in Hz
    ntimeframes::Int # number of original signal frames or nothing if not known
    names::Vector{Symbol}

    function RFFTSpectrumArray{T}(X::Matrix{T}, rate::Float64, ntimeframes::Int, names::Vector{Symbol}) where T
        ntimeframes >> 1 + 1 == size(X, 1) || throw(DimensionMismatch("ntimeframes >> 1 + 1 = $(ntimeframes >> 1 + 1) != size(X, 1) = $(size(X, 1))"))
        _check_channel_names(names)
        length(names) == size(X, 2) || throw(DimensionMismatch("the number of names ($(length(names))) does not match the number of channels ($(size(X, 2)))!"))
        new{T}(X, rate, ntimeframes, names)
    end
end

RFFTSpectrumArray(X::AbstractMatrix{T}, rate::Frequency, ntimeframes::Integer, names::Union{Nothing,Vector{Symbol}}=nothing) where T =
    RFFTSpectrumArray{T}(X, toHz(rate), Int(ntimeframes), isnothing(names) ? _default_channel_names(size(X, 2)) : names)

function RFFTSpectrumArray(X::AbstractSpectrumArray{T}, nframes_::Integer, odd::Bool=true, extrapolation_bc=Line()) where T
    domX = domain(X)
    domY = range(0, rate(X) / 2; length=nframes_)
    Xu = unwrap(X)
    data_old = data(Xu)
    data_new = Matrix{T}(undef, nframes_, nchannels(X))
    for i in 1:nchannels(X)
        mags = LinearInterpolation(domX, abs.(data_old[:, i]); extrapolation_bc=extrapolation_bc)(domY)
        mags[mags .< 0] .= 0
        phis = LinearInterpolation(domX, angle.(data_old[:, i]); extrapolation_bc=extrapolation_bc)(domY)
        data_new[:, i] .= (MagPhase(m, p) for (m, p) in zip(mags, phis))
    end
    RFFTSpectrumArray{T}(data_new, rate(X), 2nframes_ - (odd ? 1 : 2), names(X))
end

ntimeframes(X::RFFTSpectrumArray) = X.ntimeframes
domain(X::RFFTSpectrumArray) = range(0, rate(X) / 2; length=nframes(X))
toindex(X::RFFTSpectrumArray{T}, t::Frequency) where T = round(Int, (nframes(X) - 1) * 2toHz(t) / rate(X)) + 1

data_no0(X::RFFTSpectrumArray) = @view data(X)[2:end, :]
domain_no0(X::RFFTSpectrumArray) = @view domain(X)[2:end]

function Base.similar(X::RFFTSpectrumArray, t::Type{T}, dims::Dims) where T
    # tries to copy channel names
    # if there are fever names in the source array use default ones
    # changing number of frames (frequencies) looses information on number of original timeframes.
    ns = dims[2] ≤ nchannels(X) ? X.names[1:dims[2]] : _default_channel_names(dims[2])
    if dims[1] != nframes(X)
        throw(ArgumentError("can not change the number of frequencies!"))
    else
        return RFFTSpectrumArray(similar(data(X), t, dims), rate(X), ntimeframes(X), ns)
    end
end

# redefined so frequencies & channel names are treated
Base.getindex(X::RFFTSpectrumArray{T}, I::R, J::S) where {T, R <: FrameIndex, S <: ChannelIndex} = throw(ArgumentError("frame index ($R) not allowed, use Colon!"))

function Base.getindex(X::RFFTSpectrumArray{T}, ::Colon, J::S) where {T, S <: ChannelIndex}
    J2 = tochannelidx(X, J)
    names_ = names(X)[J2]
    data_ = data(X)[:, J2]
    return RFFTSpectrumArray(data_, rate(X), ntimeframes(X), names_)
end

function Base.hcat(X::RFFTSpectrumArray...)
    length(unique(rate.(X))) == 1 || throw(ArgumentError("hcat: non-unique rates!"))
    length(unique(nframes.(X))) == 1 || throw(DimensionMismatch("hcat: non-unique number of frames/frequencies!"))
    length(unique(ntimeframes.(X))) == 1 || throw(DimensionMismatch("hcat: non-unique number of original frames!"))
    newnames = _unique_channel_names(X...)
    data_ = hcat(map(data, X)...)
    return eltype(X)(data_, rate(X[1]), ntimeframes(X[1]), newnames) # eltype gives common supertype
end

function Base.vcat(X::RFFTSpectrumArray...)
    throw(ErrorException("vcat: no support for $(eltype(X))"))
end

function Base.show(io::IO, ::MIME"text/plain", X::RFFTSpectrumArray{T}) where T
    d = domain(X)
    step = d[2] - d[1]
    nframesstr = " from $(ntimeframes(X)) frames"
    clist = join([":$(n)" for n in names(X)], ", ")
    print(io, "RFFTSpectrumArray{$T}: $(nchannels(X)) channels: $(clist)\n $(nframes(X)) freqs $(first(d)) - $(last(d)) Hz, step ≈ $(step) Hz sampled at $(rate(X)) Hz$(nframesstr):\n ")
    print(io, data(X))
end

function resamplefreq(X::RFFTSpectrumArray, n::Int, extrapolation_bc=Line())
    Y = similar(X, n, nchannels(X))
    Xu = unwrap(X)
    domX = domain(Xu)
    domY = domain(Y)

    for i in 1:nchannels(Xu)
        # mags = CubicSplineInterpolation(domX, abs.(data(Xu)[:, i]); extrapolation_bc=extrapolation_bc)(domY) # can cause negative amplitudes
        mags = LinearInterpolation(domX, abs.(data(Xu)[:, i]); extrapolation_bc=extrapolation_bc)(domY)
        phis = CubicSplineInterpolation(domX, angle.(data(Xu)[:, i]); extrapolation_bc=extrapolation_bc)(domY)
        data(Y)[:, i] .= (MagPhase(m, p) for (m, p) in zip(mags, phis))
    end
    Y
end


text:
mutable struct RandomOCClassifier <: OCClassifier
    data
    outlier_bias
    RandomOCClassifier(data) = new(data, 0.1)
    RandomOCClassifier(data, pools::Vector) = RandomOCClassifier(data)
    RandomOCClassifier(data, outlier_bias::Float64) = new(data, outlier_bias)
    RandomOCClassifier(data, outlier_bias::Float64, pools::Vector) = RandomOCClassifier(data, outlier_bias)
end

set_adjust_K!(model::RandomOCClassifier, adjust_K::Bool) = nothing

get_model_params(model::RandomOCClassifier) = (:Ϟ => 0.42, :param2 => 5)

function fit!(model::RandomOCClassifier, solver)
    debug(LOGGER, "[FIT] $(typeof(model)) always returns :Optimal.")
    return :Optimal
end

function predict(model::RandomOCClassifier, target::Array{T,2}) where T <: Real
    rand(size(target, 2)) * 2 .- model.outlier_bias * 2
end

initialize!(model::RandomOCClassifier, strategy::InitializationStrategy) = nothing

get_kernel(model::RandomOCClassifier) = nothing

invalidate_solution!(model::RandomOCClassifier) = nothing

set_pools!(model::RandomOCClassifier, pools::Vector{Symbol}) = nothing
set_pools!(model::RandomOCClassifier, pools::Dict{Symbol, Vector{Int}}) = nothing

function set_data!(model::RandomOCClassifier, data::Array{T, 2}) where T <: Real
    model.data = data
    return nothing
end


text:
##### Beginning of file

PredictMDExtra.import_all()
PredictMDExtra.import_all(Main)

PredictMDExtra._import_all_on_init()
PredictMDExtra._import_all_on_init(Main)

##### End of file


text:
# ---
# title: 762. Prime Number of Set Bits in Binary Representation
# id: problem762
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Bit Manipulation
# link: <https://leetcode.com/problems/prime-number-of-set-bits-in-binary-representation/description/>
# hidden: true
# ---
#
# Given two integers `L` and `R`, find the count of numbers in the range `[L,
# R]` (inclusive) having a prime number of set bits in their binary
# representation.
#
# (Recall that the number of set bits an integer has is the number of `1`s
# present when written in binary. For example, `21` written in binary is `10101`
# which has 3 set bits. Also, 1 is not a prime.)
#
# **Example 1:**
#
#
#
#     Input: L = 6, R = 10
#     Output: 4
#     Explanation:
#     6 -> 110 (2 set bits, 2 is prime)
#     7 -> 111 (3 set bits, 3 is prime)
#     9 -> 1001 (2 set bits , 2 is prime)
#     10->1010 (2 set bits , 2 is prime)
#
#
# **Example 2:**
#
#
#
#     Input: L = 10, R = 15
#     Output: 5
#     Explanation:
#     10 -> 1010 (2 set bits, 2 is prime)
#     11 -> 1011 (3 set bits, 3 is prime)
#     12 -> 1100 (2 set bits, 2 is prime)
#     13 -> 1101 (3 set bits, 3 is prime)
#     14 -> 1110 (3 set bits, 3 is prime)
#     15 -> 1111 (4 set bits, 4 is not prime)
#
#
# **Note:**
#
#   1. `L, R` will be integers `L <= R` in the range `[1, 10^6]`.
#   2. `R - L` will be at most 10000.
#
#
## @lc code=start
using LeetCode

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


text:
using Distributed
using Dates
import REPL
using Printf: @sprintf

# parse some command-line arguments
function extract_flag!(args, flag, default=nothing)
    for f in args
        if startswith(f, flag)
            # Check if it's just `--flag` or if it's `--flag=foo`
            if f != flag
                val = split(f, '=')[2]
                if default !== nothing && !(typeof(default) <: AbstractString)
                  val = parse(typeof(default), val)
                end
            else
                val = default
            end

            # Drop this value from our args
            filter!(x -> x != f, args)
            return (true, val)
        end
    end
    return (false, default)
end
do_help, _ = extract_flag!(ARGS, "--help")
if do_help
    println("""
        Usage: runtests.jl [--help] [--list] [--jobs=N] [TESTS...]

               --help             Show this text.
               --list             List all available tests.
               --quickfail        Fail the entire run as soon as a single test errored.
               --jobs=N           Launch `N` processes to perform tests (default: Threads.nthreads()).

               Remaining arguments filter the tests that will be executed.""")
    exit(0)
end
_, jobs = extract_flag!(ARGS, "--jobs", Threads.nthreads())
do_quickfail, _ = extract_flag!(ARGS, "--quickfail")

include("setup.jl")     # make sure everything is precompiled
@test oneAPI.functional()

# choose tests
const tests = []
const test_runners = Dict()
## files in the test folder
for (rootpath, dirs, files) in walkdir(@__DIR__)
  # find Julia files
  filter!(files) do file
    endswith(file, ".jl") && file !== "setup.jl" && file !== "runtests.jl"
  end
  isempty(files) && continue

  # strip extension
  files = map(files) do file
    file[1:end-3]
  end

  # prepend subdir
  subdir = relpath(rootpath, @__DIR__)
  if subdir != "."
    files = map(files) do file
      joinpath(subdir, file)
    end
  end

  append!(tests, files)
  for file in files
    test_runners[file] = ()->include("$(@__DIR__)/$file.jl")
  end
end
## GPUArrays testsuite
for name in keys(TestSuite.tests)
    push!(tests, "gpuarrays$(Base.Filesystem.path_separator)$name")
    test_runners["gpuarrays$(Base.Filesystem.path_separator)$name"] = ()->TestSuite.tests[name](oneArray)
end
unique!(tests)

# parse some more command-line arguments
## --list to list all available tests
do_list, _ = extract_flag!(ARGS, "--list")
if do_list
    println("Available tests:")
    for test in sort(tests)
        println(" - $test")
    end
    exit(0)
end
## no options should remain
optlike_args = filter(startswith("-"), ARGS)
if !isempty(optlike_args)
    error("Unknown test options `$(join(optlike_args, " "))` (try `--help` for usage instructions)")
end
## the remaining args filter tests
if !isempty(ARGS)
  filter!(tests) do test
    any(arg->startswith(test, arg), ARGS)
  end
end

# add workers
const test_exeflags = Base.julia_cmd()
filter!(test_exeflags.exec) do c
    return !(startswith(c, "--depwarn") || startswith(c, "--check-bounds"))
end
push!(test_exeflags.exec, "--check-bounds=yes")
push!(test_exeflags.exec, "--startup-file=no")
push!(test_exeflags.exec, "--depwarn=yes")
if Base.JLOptions().project != C_NULL
    push!(test_exeflags.exec, "--project=$(unsafe_string(Base.JLOptions().project))")
end
const test_exename = popfirst!(test_exeflags.exec)
function addworker(X; kwargs...)
    withenv("JULIA_NUM_THREADS" => 1, "OPENBLAS_NUM_THREADS" => 1) do
        procs = addprocs(X; exename=test_exename, exeflags=test_exeflags, kwargs...)
        @everywhere procs include($(joinpath(@__DIR__, "setup.jl")))
        procs
    end
end
addworker(min(jobs, length(tests)))

# pretty print information about gc and mem usage
testgroupheader = "Test"
workerheader = "(Worker)"
name_align        = maximum([textwidth(testgroupheader) + textwidth(" ") +
                             textwidth(workerheader); map(x -> textwidth(x) +
                             3 + ndigits(nworkers()), tests)])
elapsed_align     = textwidth("Time (s)")
gc_align      = textwidth("GC (s)")
percent_align = textwidth("GC %")
alloc_align   = textwidth("Alloc (MB)")
rss_align     = textwidth("RSS (MB)")
printstyled(" "^(name_align + textwidth(testgroupheader) - 3), " | ")
printstyled("         | ---------------- CPU ---------------- |\n", color=:white)
printstyled(testgroupheader, color=:white)
printstyled(lpad(workerheader, name_align - textwidth(testgroupheader) + 1), " | ", color=:white)
printstyled("Time (s) | GC (s) | GC % | Alloc (MB) | RSS (MB) |\n", color=:white)
print_lock = stdout isa Base.LibuvStream ? stdout.lock : ReentrantLock()
if stderr isa Base.LibuvStream
    stderr.lock = print_lock
end
function print_testworker_stats(test, wrkr, resp)
    @nospecialize resp
    lock(print_lock)
    try
        printstyled(test, color=:white)
        printstyled(lpad("($wrkr)", name_align - textwidth(test) + 1, " "), " | ", color=:white)
        time_str = @sprintf("%7.2f",resp[2])
        printstyled(lpad(time_str, elapsed_align, " "), " | ", color=:white)

        cpu_gc_str = @sprintf("%5.2f", resp[4])
        printstyled(lpad(cpu_gc_str, gc_align, " "), " | ", color=:white)
        # since there may be quite a few digits in the percentage,
        # the left-padding here is less to make sure everything fits
        cpu_percent_str = @sprintf("%4.1f", 100 * resp[4] / resp[2])
        printstyled(lpad(cpu_percent_str, percent_align, " "), " | ", color=:white)
        cpu_alloc_str = @sprintf("%5.2f", resp[3] / 2^20)
        printstyled(lpad(cpu_alloc_str, alloc_align, " "), " | ", color=:white)

        cpu_rss_str = @sprintf("%5.2f", resp[6] / 2^20)
        printstyled(lpad(cpu_rss_str, rss_align, " "), " |\n", color=:white)
    finally
        unlock(print_lock)
    end
end
global print_testworker_started = (name, wrkr)->begin
end
function print_testworker_errored(name, wrkr)
    lock(print_lock)
    try
        printstyled(name, color=:red)
        printstyled(lpad("($wrkr)", name_align - textwidth(name) + 1, " "), " |",
            " "^elapsed_align, " failed at $(now())\n", color=:red)
    finally
        unlock(print_lock)
    end
end

# run tasks
results = []
all_tasks = Task[]
all_tests = copy(tests)
try
    # Monitor stdin and kill this task on ^C
    # but don't do this on Windows, because it may deadlock in the kernel
    t = current_task()
    running_tests = Dict{String, DateTime}()
    if !Sys.iswindows() && isa(stdin, Base.TTY)
        stdin_monitor = @async begin
            term = REPL.Terminals.TTYTerminal("xterm", stdin, stdout, stderr)
            try
                REPL.Terminals.raw!(term, true)
                while true
                    c = read(term, Char)
                    if c == '\x3'
                        Base.throwto(t, InterruptException())
                        break
                    elseif c == '?'
                        println("Currently running: ")
                        tests = sort(collect(running_tests), by=x->x[2])
                        foreach(tests) do (test, date)
                            println(test, " (running for ", round(now()-date, Minute), ")")
                        end
                    end
                end
            catch e
                isa(e, InterruptException) || rethrow()
            finally
                REPL.Terminals.raw!(term, false)
            end
        end
    end
    @sync begin
        for p in workers()
            @async begin
                push!(all_tasks, current_task())
                while length(tests) > 0
                    test = popfirst!(tests)
                    local resp
                    wrkr = p

                    # run the test
                    running_tests[test] = now()
                    try
                        resp = remotecall_fetch(runtests, wrkr, test_runners[test], test)
                    catch e
                        isa(e, InterruptException) && return
                        resp = Any[e]
                    end
                    delete!(running_tests, test)
                    push!(results, (test, resp))

                    # act on the results
                    if resp[1] isa Exception
                        print_testworker_errored(test, wrkr)
                        do_quickfail && Base.throwto(t, InterruptException())

                        # the worker encountered some failure, recycle it
                        # so future tests get a fresh environment
                        rmprocs(wrkr, waitfor=30)
                        p = addworker(1)[1]
                    else
                        print_testworker_stats(test, wrkr, resp)
                    end
                end
                if p != 1
                    # Free up memory =)
                    rmprocs(p, waitfor=30)
                end
            end
        end
    end
catch e
    isa(e, InterruptException) || rethrow()
    # If the test suite was merely interrupted, still print the
    # summary, which can be useful to diagnose what's going on
    foreach(task -> begin
            istaskstarted(task) || return
            istaskdone(task) && return
            try
                schedule(task, InterruptException(); error=true)
            catch ex
                @error "InterruptException" exception=ex,catch_backtrace()
            end
        end, all_tasks)
    for t in all_tasks
        # NOTE: we can't just wait, but need to discard the exception,
        #       because the throwto for --quickfail also kills the worker.
        try
            wait(t)
        catch e
            showerror(stderr, e)
        end
    end
finally
    if @isdefined stdin_monitor
        schedule(stdin_monitor, InterruptException(); error=true)
    end
end

# construct a testset to render the test results
o_ts = Test.DefaultTestSet("Overall")
Test.push_testset(o_ts)
completed_tests = Set{String}()
for (testname, (resp,)) in results
    push!(completed_tests, testname)
    if isa(resp, Test.DefaultTestSet)
        Test.push_testset(resp)
        Test.record(o_ts, resp)
        Test.pop_testset()
    elseif isa(resp, Tuple{Int,Int})
        fake = Test.DefaultTestSet(testname)
        for i in 1:resp[1]
            Test.record(fake, Test.Pass(:test, nothing, nothing, nothing))
        end
        for i in 1:resp[2]
            Test.record(fake, Test.Broken(:test, nothing))
        end
        Test.push_testset(fake)
        Test.record(o_ts, fake)
        Test.pop_testset()
    elseif isa(resp, RemoteException) && isa(resp.captured.ex, Test.TestSetException)
        println("Worker $(resp.pid) failed running test $(testname):")
        Base.showerror(stdout, resp.captured)
        println()
        fake = Test.DefaultTestSet(testname)
        for i in 1:resp.captured.ex.pass
            Test.record(fake, Test.Pass(:test, nothing, nothing, nothing))
        end
        for i in 1:resp.captured.ex.broken
            Test.record(fake, Test.Broken(:test, nothing))
        end
        for t in resp.captured.ex.errors_and_fails
            Test.record(fake, t)
        end
        Test.push_testset(fake)
        Test.record(o_ts, fake)
        Test.pop_testset()
    else
        if !isa(resp, Exception)
            resp = ErrorException(string("Unknown result type : ", typeof(resp)))
        end
        # If this test raised an exception that is not a remote testset exception,
        # i.e. not a RemoteException capturing a TestSetException that means
        # the test runner itself had some problem, so we may have hit a segfault,
        # deserialization errors or something similar.  Record this testset as Errored.
        fake = Test.DefaultTestSet(testname)
        Test.record(fake, Test.Error(:test_error, testname, nothing, Any[(resp, [])], LineNumberNode(1)))
        Test.push_testset(fake)
        Test.record(o_ts, fake)
        Test.pop_testset()
    end
end
for test in all_tests
    (test in completed_tests) && continue
    fake = Test.DefaultTestSet(test)
    Test.record(fake, Test.Error(:test_interrupted, test, nothing,
                                    [("skipped", [])], LineNumberNode(1)))
    Test.push_testset(fake)
    Test.record(o_ts, fake)
    Test.pop_testset()
end
println()
Test.print_test_results(o_ts, 1)
if !o_ts.anynonpass
    println("    \033[32;1mSUCCESS\033[0m")
else
    println("    \033[31;1mFAILURE\033[0m\n")
    Test.print_test_errors(o_ts)
    throw(Test.FallbackTestSetException("Test run finished with errors"))
end