text:
#--------------------------------------------------
# Imports
#--------------------------------------------------

using JuliaPkgsList

using Test
using JSON

#--------------------------------------------------
# Aux values and functions
#--------------------------------------------------

# String → String
testFilePath(path :: String) = joinpath(@__DIR__, TEST_FILES_DIR, path)

const TEST_FILES_DIR = "test-files"
const PKGS_INFO_FILE = testFilePath("pkgs-info.json")
const SAMPLE_FILE    = testFilePath("sample-pkgs-info.json")

const JULIA_REPO = "https://github.com/JuliaLang/julia.git"
const JULIA_NAME = "julia"

tryrm(path :: AbstractString) =
    try rm(path); catch err @warn(err) end

#--------------------------------------------------
# Tests
#--------------------------------------------------

include("utils.jl")

@testset "JuliaPkgsList.jl :: main processing function" begin
    pkgs = loadPkgsData(SAMPLE_FILE)
    # main function
    namesSorted = getSortedPkgsInfo(pkgs, getName)
    @test namesSorted[1] == "julia"
    @test namesSorted[2] == "JuMP"
    reposSortedNoJulia = getSortedPkgsInfo(pkgs, getRepo; excludedRepos=[JULIA_REPO])
    @test reposSortedNoJulia[1] == "https://github.com/jump-dev/JuMP.jl.git"
    # aux repo/name function
    @test getSortedPkgs(pkgs)[3] == "https://github.com/FluxML/Zygote.jl.git"
    @test getSortedPkgs(pkgs, true; excludedRepos=[JULIA_REPO])[2] == "Zygote"
end

@testset "JuliaPkgsList.jl :: main entry function" begin
    # names with version, excluded packages
    fname1 = testFilePath("top-5-pkgs-list.txt")
    genTopPkgsList(5;
        pkgsListOutputFile=fname1, addPkgNum=false,
        showName=true, includeVersion=true,
        pkgsInfoFile=SAMPLE_FILE
    )
    @test isfile(fname1)
    @test read(fname1, String) ==
        "JuMP, 0.21.5\nZygote, 0.6.0\nJuliaDB, 0.13.1\nJuliaZH, 1.5.4\nJavis, 0.4.0"
    tryrm(fname1)
    # repos, empty list of excluded packages
    fname2 = testFilePath("pkgs-list.txt")
    genTopPkgsList(4;
        pkgsListOutputFile=fname2, pkgsInfoFile=SAMPLE_FILE,
        excludedReposFile=testFilePath("noexcluded.txt")
    )
    fname2 = testFilePath("4-pkgs-list.txt")
    @test isfile(fname2)
    @test read(fname2, String) == join([
        "https://github.com/JuliaLang/julia.git", "https://github.com/jump-dev/JuMP.jl.git",
        "https://github.com/FluxML/Zygote.jl.git", "https://github.com/JuliaData/JuliaDB.jl.git"
    ], "\n")
    tryrm(fname2)
end

include("pkgs-API.jl")


text:
using MakieThemes
using Test
using Makie
using CSV
using DataFrames
using StatsMakie

@testset "simple" begin
  for i in ggthemr_colorthemes()
        scene = show_ggthemr(i)
        save("img/$i.png", scene)
  end
end

@testset "original" begin
  for dataset ∈ (:www, :drivers, :mtcars, :diamonds)
    @eval const $(dataset) = CSV.read("data/"*$(string(dataset))*".tsv", delim = '\t', allowmissing = :none)
  end

  AbstractPlotting.set_theme!(ggthemr(:fresh))

  p1 = scatterlines(Data(www), :Minute, :Users,
    Group(color = :Measure, marker = :Measure),
    markersize = 6, marker = [:rect, :circle]);

  p2 = plot(density, Data(mtcars),
    :mpg, Group(color = :cyl));

  p3 = plot(Position.stack, histogram, Data(diamonds),
    :price, Group(color = :cut));

  p4 = boxplot(Data(drivers), :Year, :Deaths);

  vbox(hbox(p1, p2), hbox(p3, p4))
end

    driver_x <- scale_x_discrete(breaks = seq(1969, 1984, 3), label = function(x) paste0("'", substr(x, 3, 4)))
    display_theme <- list(theme(legend.position = 'none'), scale_y_continuous(label = function(x) format(x, big.mark = ",", scientific = FALSE)))

      ggplot(transform(mtcars, cyl = factor(cyl)), aes_string(x = 'mpg', fill = 'cyl', colour = 'cyl')) +
        geom_density(alpha = 0.75) +
        labs(fill = 'Cylinders', colour = 'Cylinders', x = 'MPG', y = 'Density') +
        legend_top(),

      ggplot(transform(diamonds, cut = factor(cut)), aes_string(x = 'price', fill = 'cut')) +
        geom_histogram(binwidth = 850) +
        display_theme +
        xlab('Price (USD)') +
        ylab('Count') +
        scale_x_continuous(label = function(x) paste0(x / 1000, 'k')),

      ggplot(drivers, aes_string(x = 'Year', y = 'Deaths')) +
        geom_boxplot(size = 0.25) +
        ylab('Monthly Deaths') +
        display_theme +
        driver_x

    )

    assign(".Random.seed", old_seed, envir=globalenv())
    do.call(grid.arrange,  display)

  }

ext: jl


text:
using DataDeps

register(DataDep(
  "lenses",
  "http://archive.ics.uci.edu/ml/datasets/Lenses",
  "http://archive.ics.uci.edu/ml/machine-learning-databases/lenses/lenses.data",
  "711f22d5a4bed062bb715d39d23d54d2667206fbd97c94c7fd013ede47fb96c2",
  post_fetch_method=(path -> begin
    UCIData.process_dataset(path,
        target_index=6,
        feature_indices=2:5,
        id_indices=1,
        categoric_indices=2:5,
        delim=' ',
        ignorerepeated=true,
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
@symbol_func function alpha_R_star(cur_reactor::AbstractReactor)
  cur_scaling = cur_reactor.mode_scaling

  cur_alpha_star = cur_scaling[:R]

  cur_alpha_star += cur_scaling[:a]

  cur_alpha_star += cur_scaling[:P]

  cur_alpha_star -= 2 * ( 1 + cur_scaling[:n] )

  cur_alpha_star
end


text:
export quasi_anechoic, quasi_anechoic_window, boundary_delay

function quasi_anechoic_window(a::SampleArray, fstart=irfirstbig; post, pre=125ms)
    @assert nchannels(a) == 1
    imax = fstart(a)
    ipre = toindexdelta(a, pre) # index deltas
    ipost = toindexdelta(a, post)
    n = ipost + ipre + 1
    # @show imax, n, ipre, ipost
    # TODO make the two-sided Tukey window properly
    win = zeros(Float64, n)
    win[1:ipre+1] .= tukey(ipre+1, 0.25)
    win[end-ipost:end] .= tukey(ipost+1, 0.25)
    win[ipre÷2:end-ipost÷2] .= 1.0
    fullwin = similar(a)
    fullwin[imax-ipre:imax+ipost, 1:1] .= win
    fullwin[1:imax-ipre-1] .= 0
    fullwin[imax+ipost+1:end] .= 0

    (win=fullwin, imax=imax, ipre=ipre, ipost=ipost)
end

function quasi_anechoic(a::SampleArray, fstart=irfirstbig; post, pre=125ms, cut=true)
    win, imax, ipre, ipost = quasi_anechoic_window(a, fstart; post=post, pre=pre)
    q = a .* win
    if cut
        return q[imax-ipre:imax+ipost, :]
    else
        return q
    end
end

# TODO make c a general constant
function boundary_delay(micheight, micdist; c=342.0)
    refdist = 2sqrt(micheight^2 + (micdist/2)^2) # distance of the reflection
    deltadist = refdist - micdist
    directtime = micdist/c
    reftime = refdist/c
    delay = deltadist/c
    (minfreq=1/delay, delay=delay, directtime=directtime, reftime=reftime, deltadist=deltadist, refdist=refdist)
end


text:
using Printf
function vtkWriteHex(filename,X,Connectivity,c,cNames)

precision = "2";
c .= c .* (abs.(c) .> 1.e-20);
open(filename, "w") do fid
  print(fid, "# vtk DataFile Version 2.0 \n");
  print(fid, "Unstructured Grid Example \n");
  print(fid, "ASCII \n");
  print(fid, "DATASET UNSTRUCTURED_GRID \n");
  @printf(fid, "POINTS %d float \n",size(X,2));

  for i in 1:size(X,2)
    @printf(fid, " %6e %6e %6e \n", X[:,i]...);
  end
  @printf(fid, "\n");
  @printf(fid, "CELLS %d %d \n",size(Connectivity,2),9*size(Connectivity,2));
  for i in 1:size(Connectivity, 2)
    @printf(fid, "8 %d %d %d %d %d %d %d %d\n",(Connectivity[:,i] .- 1)...);
  end
  @printf(fid, "CELL_TYPES %d \n",size(Connectivity,2));
  for i in 1:size(Connectivity,2)
    @printf(fid, "%d  \n",12);
  end
  @printf(fid, "Cell_DATA %d \n",size(Connectivity,2));
  for i=1:size(c,2)
    @printf(fid, "SCALARS %s double 1 \n",cNames[i]);
    print(fid, "LOOKUP_TABLE default \n");
    for j=1:4:size(c,1)-3
      @printf(fid, "%6e %6e %6e %6e \n",c[j:j+3,i]...);
    end
  end
end

end


text:
using FMMLIB2D
using Test

@testset "Laplace FMM 2D" begin
    include("laplace.jl")
end

@testset "Helmholtz FMM 2D" begin
    include("helmholtz.jl")
end

@testset "Complex FMM 2D" begin
    include("complex.jl")
end


text:
function request_handler(req::HTTP.Request)
    target = String(HTTP.Messages.getfield(req, :target))
    uri = HTTP.URI(target)
    query = uri.query
    urisstring = unescapeuri(query)
    params = HTTP.queryparams(urisstring) # create URI out of target String
    wewant = ["access_token", "token_type", "state", "expires_in"]
    ipa = filter(pa-> in(pa[1], wewant), params)
    if length(ipa) > 3
        if get(ipa, "state", "") == "987"
            if get(ipa, "token_type", "") == "Bearer"
                c = spotcred()::SpotifyCredentials
                c.ig_access_token = get(ipa, "access_token", "")
                c.token_type = get(ipa, "token_type", "")
                expi_in = parse(Int, get(ipa, "expires_in", 0))
                c.expires_at = string(Dates.now() + Dates.Second(expi_in))
                SPOTCRED[] = c
                # Feedback in 'wait_for_ig_access'
            else
                @error "Unexpected token_type"
            end
        else
            @error "Unexpected state"
        end
        HTTP.Response(202, repr(ipa))
    elseif length(ipa) > 0
        @error "Unexpected query string: $query"
        HTTP.Response(202, repr(ipa))
    else
        if target == "/favicon.ico"
            HTTP.Response(404)
        else
            fpath = joinpath(@__DIR__, "standardize_uri.html")
            HTTP.Response( 200, body = read(fpath ) )
        end
    end
end

function has_ig_access_token()
    c = spotcred()::SpotifyCredentials
    t = c.ig_access_token
    t != ""
end
function wait_for_ig_access(; stopwait::DateTime = now() + Dates.Second(20))
    while !has_ig_access_token() && now() < stopwait
        sleep(0.1)
    end
    now() >= stopwait  ? println("\n\tTimeout server") : println("\n\tImplicit grant token expires at $(SPOTCRED[].expires_at) - closing server")
end

function start_a_listening_server()
    stored = spotcred().redirect
    stripped = strip(stored, '/')
    u = parse(URI, stripped)
    host() = parse(IPAddr, u.host)
    port() = parse(Int, u.port)
    path() = u.path
    server = Sockets.listen(host(), port())
    close_server_when_ready_task = @async begin
        wait_for_ig_access()
        close(server)
    end
    println("\tListening for authorization on $(host()):$(port()) and path $(path())")
    listen_task = @async HTTP.serve(request_handler, host(), port(); server, stream =false)
    listen_task, close_server_when_ready_task
end

function launch_a_browser_that_asks_for_implicit_grant()
    vs =["user-read-private",
        "user-read-email",
        "user-follow-read"]
    scopes = join(vs, "%20")
    uri = OAUTH_AUTHORIZE_URL * "?" *
            "client_id=" * spotcred().client_id *
            "&redirect_uri=" * replace(spotcred().redirect, "/" => "%2F") *
            "&scope=" * scopes *
            "&show_dialog=true" *
            "&response_type=token" *
            "&state=987"
    println("\tLaunching a browser at: $uri")
    success = false
    #=
    browsercmd = get_working_browser_cmd()
    if browsercmd != ``
        @info "Browser command from .ini-file: $browsercmd"
    end
    =#
    browser = ""
    while !success && COUNTBROWSER.value < length(BROWSERS)
        #if browsercmd != ``
        printstyled("\tTrying to launch browser candidate: $(BROWSERS[COUNTBROWSER.value + 1])\n"; color = :green)
        success, browser = open_a_browser(url= uri)
        #else
        #    printstyled("\tBrowser command from .ini-file: $browsercmd\n"; color = :green)
        #    @async run(browsercmd)
        #    break
        #end
    end
    # Reset browser picker
    COUNTBROWSER.value = 0
    #=
    if browsercmd == "" && COUNTBROWSER.value <= length(BROWSERS)
        browsercmd = launch_command(browser; url = uri)
    end
    if browsercmd != ``
        printstyled("\tStoring browser command:$(browsercmd)\n", color = :green)
        set_working_browser_cmd(browsercmd)
    end
    =#
end

"""
If the user gives us access, that will be in the future.
Is it worth waiting for? If so, call `wait_for_ig_access()`
"""
function get_implicit_grant()
    listen_task, close_server_when_ready_task = start_a_listening_server()
    launch_a_browser_that_asks_for_implicit_grant()
    listen_task, close_server_when_ready_task
end


text:
addprocs(8)

include("daggen.jl")
using DagScheduler
using Base.Test

isdir(".mempool") && rm(".mempool"; recursive=true)
@everywhere begin
    DagScheduler.META_IMPL[:node] = ENV["NODE_META_IMPL"]
    DagScheduler.META_IMPL[:cluster] = ENV["CLUSTER_META_IMPL"]
    DagScheduler.META_IMPL[:map_num_entries] = 1024*100
    DagScheduler.META_IMPL[:map_entry_sz] = 1512
end

node1 = NodeEnv(2, getipaddr(), [3,4,5])
node2 = NodeEnv(6, getipaddr(), [7,8,9])
runenv = DagScheduler.Plugin.setrunenv(RunEnv(; nodes=[node1,node2]))
RC = (DagScheduler.META_IMPL[:cluster] == "DagScheduler.FdbMeta.FdbExecutorMeta") ? DagScheduler.FdbRC : DagScheduler.SimpleRC

@testset "deep dag" begin
    @everywhere begin
        RC = (DagScheduler.META_IMPL[:cluster] == "DagScheduler.FdbMeta.FdbExecutorMeta") ? DagScheduler.FdbRC : DagScheduler.SimpleRC
        DagScheduler.RefCounter[] = RC
    end

    info("Testing deep dag...")
    dag1 = gen_straight_dag(ones(Int, 6^4))
    result = collect(rundag(runenv, dag1))
    info("result = ", result)
    @test result == 1
    @everywhere MemPool.cleanup()
    DagScheduler.print_stats(runenv)
    @test isempty(RC.nz_refcounts())

    info("Testing cross connected dag...")
    dag3 = gen_cross_dag()
    result = collect(rundag(runenv, dag3))
    info("result = ", result)
    @test result == 84
    @everywhere MemPool.cleanup()
    DagScheduler.print_stats(runenv)
    @test isempty(RC.nz_refcounts())
end

@testset "sorting" begin
    @everywhere begin
        DagScheduler.RefCounter[] = DagScheduler.NoRC
    end

    info("Testing sorting...")

    for L in (10^6, 10^7)
        dag2 = gen_sort_dag(L, 40, 4, 1)
        DagScheduler.dref_to_fref!(dag2)
        result = collect(rundag(runenv, dag2))
        info("result = ", typeof(result), ", length: ", length(result), ", sorted: ", issorted(result))
        @test isa(result, Array{Float64,1})
        @test issorted(result)
        @test length(result) == L
        @everywhere MemPool.cleanup()
        DagScheduler.print_stats(runenv)

        # for cross dag
        dag4 = gen_sort_dag(L, 40, 4, 40)
        DagScheduler.dref_to_fref!(dag4)
        result = collect(rundag(runenv, dag4))
        info("result = ", typeof(result), ", length: ", length(result))
        fullresult = collect(Dagger.treereduce(delayed(vcat), result))
        @test isa(fullresult, Array{Float64,1})
        @test issorted(fullresult)
        @test length(fullresult) == L
        @everywhere MemPool.cleanup()
        DagScheduler.print_stats(runenv)
    end
end

@testset "meta" begin
    @everywhere begin
        RC = (DagScheduler.META_IMPL[:cluster] == "DagScheduler.FdbMeta.FdbExecutorMeta") ? DagScheduler.FdbRC : DagScheduler.SimpleRC
        DagScheduler.RefCounter[] = RC
    end

    info("Testing meta annotation...")
    x = [delayed(rand)(10) for i=1:10]
    y = delayed((c...) -> [c...]; meta=true)(x...)
    result = collect(rundag(runenv, y))
    @test isa(result, Vector{<:Dagger.Chunk})
    @test length(result) == 10
    @everywhere MemPool.cleanup()
    DagScheduler.print_stats(runenv)
    @test isempty(RC.nz_refcounts())
end

DagScheduler.cleanup(runenv)
isdir(".mempool") && rm(".mempool"; recursive=true)


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct Sku <: SwaggerModel
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    tier::Any # spec type: Union{ Nothing, String } # spec name: tier
    resourceType::Any # spec type: Union{ Nothing, String } # spec name: resourceType
    kind::Any # spec type: Union{ Nothing, String } # spec name: kind
    locations::Any # spec type: Union{ Nothing, Vector{String} } # spec name: locations
    capabilities::Any # spec type: Union{ Nothing, Vector{SKUCapability} } # spec name: capabilities
    restrictions::Any # spec type: Union{ Nothing, Vector{Restriction} } # spec name: restrictions

    function Sku(;name=nothing, tier=nothing, resourceType=nothing, kind=nothing, locations=nothing, capabilities=nothing, restrictions=nothing)
        o = new()
        validate_property(Sku, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(Sku, Symbol("tier"), tier)
        setfield!(o, Symbol("tier"), tier)
        validate_property(Sku, Symbol("resourceType"), resourceType)
        setfield!(o, Symbol("resourceType"), resourceType)
        validate_property(Sku, Symbol("kind"), kind)
        setfield!(o, Symbol("kind"), kind)
        validate_property(Sku, Symbol("locations"), locations)
        setfield!(o, Symbol("locations"), locations)
        validate_property(Sku, Symbol("capabilities"), capabilities)
        setfield!(o, Symbol("capabilities"), capabilities)
        validate_property(Sku, Symbol("restrictions"), restrictions)
        setfield!(o, Symbol("restrictions"), restrictions)
        o
    end
end # type Sku

const _property_map_Sku = Dict{Symbol,Symbol}(Symbol("name")=>Symbol("name"), Symbol("tier")=>Symbol("tier"), Symbol("resourceType")=>Symbol("resourceType"), Symbol("kind")=>Symbol("kind"), Symbol("locations")=>Symbol("locations"), Symbol("capabilities")=>Symbol("capabilities"), Symbol("restrictions")=>Symbol("restrictions"))
const _property_types_Sku = Dict{Symbol,String}(Symbol("name")=>"String", Symbol("tier")=>"String", Symbol("resourceType")=>"String", Symbol("kind")=>"String", Symbol("locations")=>"Vector{String}", Symbol("capabilities")=>"Vector{SKUCapability}", Symbol("restrictions")=>"Vector{Restriction}")
Base.propertynames(::Type{ Sku }) = collect(keys(_property_map_Sku))
Swagger.property_type(::Type{ Sku }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_Sku[name]))}
Swagger.field_name(::Type{ Sku }, property_name::Symbol) =  _property_map_Sku[property_name]

const _allowed_Sku_name = ["Standard_LRS", "Standard_GRS", "Standard_RAGRS", "Standard_ZRS", "Premium_LRS"]

const _allowed_Sku_tier = ["Standard", "Premium"]

const _allowed_Sku_kind = ["Storage", "BlobStorage"]

function check_required(o::Sku)
    true
end

function validate_property(::Type{ Sku }, name::Symbol, val)
    if name === Symbol("name")
        Swagger.validate_param(name, "Sku", :enum, val, _allowed_Sku_name)
    end
    if name === Symbol("tier")
        Swagger.validate_param(name, "Sku", :enum, val, _allowed_Sku_tier)
    end
    if name === Symbol("kind")
        Swagger.validate_param(name, "Sku", :enum, val, _allowed_Sku_kind)
    end
end


text:
function readmodel2d(ParaPath::String; mor = false)
    fid = open(ParaPath,"r")
    iflag, ext = read(fid,Int64,2)
    modelsize = read(fid,Int64,2) # Model size Nz and Nx with PML boundary

    if iflag == 2 # unlimited medium
        nDZ = modelsize[1] - 2*ext
        nHX = modelsize[2] - 2*ext
    elseif iflag == 1 # free surface
        nDZ = modelsize[1] - ext
        nHX = modelsize[2] - 2*ext
    end
    BDnDZ = modelsize[1]
    BDnHX = modelsize[2]

    pvel = read(fid,Float64,modelsize[1]*modelsize[2])
    pvel = reshape(pvel,modelsize[1],modelsize[2])
    svel = read(fid,Float64,modelsize[1]*modelsize[2])
    svel = reshape(svel,modelsize[1],modelsize[2])
    rho = read(fid,Float64,modelsize[1]*modelsize[2])
    rho = reshape(rho,modelsize[1],modelsize[2])
    lambda = pvel.^2.*rho - 2*svel.^2.*rho
    mu = svel.^2.*rho
    dx,dz,dt = read(fid,Float64,3)
    DZ = nDZ*dz
    HX = nHX*dx
    nT,sn = read(fid,Int64,2)
    T = nT*dt
    pkf = read(fid,Float64,1)
    medium = medium2d(pvel, svel, rho, lambda, mu, dx, dz, dt, DZ, HX,
        nDZ, nHX, BDnDZ, BDnHX,T,nT,pkf[1],ext,iflag)

    BDntpp = modelsize
    BDntxz = [BDntpp[1]-1, BDntpp[2]-1]
    BDnvx = [BDntpp[1], BDntpp[2]-1]
    BDnvz = [BDntpp[1]-1, BDntpp[2]]
    if iflag == 1
        ntpp = [BDntpp[1]-ext BDntpp[2]-2*ext]
        ntxz = [BDntxz[1]-ext, BDntxz[2]-2*ext]
        nvx = [BDnvx[1]-ext, BDnvx[2]-2*ext]
        nvz = [BDnvz[1]-ext, BDnvz[2]-2*ext]
    elseif iflag == 2
        ntpp = [BDntpp[1]-2*ext, BDntpp[2]-2*ext]
        ntxz = [BDntxz[1]-2*ext, BDntxz[2]-2*ext]
        nvx = [BDnvx[1]-2*ext, BDnvx[2]-2*ext]
        nvz = [BDnvz[1]-2*ext, BDnvz[2]-2*ext]
    end
    nwf = nwf2d(BDntpp,BDntxz,BDnvx,BDnvz,ntpp,ntxz,nvx,nvz)
    wf = initwf(nDZ,nHX,ext,iflag,mor)

    if mor == true
        FDC = FDCoeff(4)
        fd = fdmtx(nwf,medium,FDC,ext,iflag)
        return spmod2d(medium, wf, nwf, fd)
    elseif mor == false
        FDC = FDCoeff(4)
        fd = fdmtx(nwf,medium,FDC,ext)
        pml = BD(nwf, medium, ext, iflag)
        return nspmod2d(medium, wf, nwf, fd, pml)
    end
    close(fid)
end


text:
abstract type SteadyStateDiffEqAlgorithm <: DiffEqBase.AbstractSteadyStateAlgorithm end

struct SSRootfind{F} <: SteadyStateDiffEqAlgorithm
  nlsolve::F
end
SSRootfind(;nlsolve=(f,u0,abstol) -> (res=NLsolve.nlsolve(f,u0,ftol = abstol);res.zero)) = SSRootfind(nlsolve)

struct DynamicSS{A,AT,RT,TS} <: SteadyStateDiffEqAlgorithm
  alg::A
  abstol::AT
  reltol::RT
  tspan::TS
end
DynamicSS(alg; abstol = 1e-8, reltol = 1e-6, tspan = Inf) =
  DynamicSS(alg, abstol, reltol, tspan)

# Backward compatibility:
DynamicSS(alg, abstol, reltol) = DynamicSS(alg; abstol=abstol, reltol=reltol)

## SciMLBase Trait Definitions

SciMLBase.isadaptive(alg::SteadyStateDiffEqAlgorithm) = true

SciMLBase.isautodifferentiable(alg::SSRootfind) = true
SciMLBase.allows_arbitrary_number_types(alg::SSRootfind) = true
SciMLBase.allowscomplex(alg::SSRootfind) = true

SciMLBase.isautodifferentiable(alg::DynamicSS) = SciMLBase.isautodifferentiable(alg.alg)
SciMLBase.allows_arbitrary_number_types(alg::DynamicSS) = SciMLBase.allows_arbitrary_number_types(alg.alg)
SciMLBase.allowscomplex(alg::DynamicSS) = SciMLBase.allowscomplex(alg.alg)


text:
"""
    SoftSphere(nl_only)

The soft-sphere potential.
"""
struct SoftSphere{S, C} <: GeneralInteraction
    cutoff::C
    nl_only::Bool
end

SoftSphere{S}(cutoff, nl_only) where S =
    SoftSphere{S, typeof(cutoff)}(cutoff, nl_only)

SoftSphere(cutoff, nl_only) =
    SoftSphere{false, typeof(cutoff)}(cutoff, nl_only)

SoftSphere(nl_only=false) =
    SoftSphere(ShiftedPotentialCutoff(3.0), nl_only)

@inline @inbounds function force(inter::SoftSphere{S, C},
                                    coord_i,
                                    coord_j,
                                    atom_i,
                                    atom_j,
                                    box_size) where {S, C}
    dr = vector(coord_i, coord_j, box_size)
    r2 = sum(abs2, dr)

    if !S && iszero(atom_i.σ) || iszero(atom_j.σ)
        return zero(coord_i)
    end

    σ = sqrt(atom_i.σ * atom_j.σ)
    ϵ = sqrt(atom_i.ϵ * atom_j.ϵ)

    cutoff = inter.cutoff
    σ2 = σ^2
    params = (σ2, ϵ)

    if cutoff_points(C) == 0
        f = force_nocutoff(inter, r2, inv(r2), params)
    elseif cutoff_points(C) == 1
        sqdist_cutoff = cutoff.sqdist_cutoff * σ2
        r2 > sqdist_cutoff && return zero(coord_i)

        f = force_cutoff(cutoff, r2, inter, params)
    elseif cutoff_points(C) == 2
        sqdist_cutoff = cutoff.sqdist_cutoff * σ2
        activation_dist = cutoff.activation_dist * σ2

        r2 > sqdist_cutoff && return zero(coord_i)

        if r2 < activation_dist
            f = force_nocutoff(inter, r2, inv(r2), params)
        else
            f = force_cutoff(cutoff, r2, inter, params)
        end
    end

    return f * dr
end

@fastmath function force_nocutoff(::SoftSphere, r2, invr2, (σ2, ϵ))
    six_term = (σ2 * invr2) ^ 3

    return (24ϵ * invr2) * 2 * six_term ^ 2
end

@inbounds function potential_energy(inter::SoftSphere{S, C},
                                    s::Simulation,
                                    i::Integer,
                                    j::Integer) where {S, C}
    U = eltype(s.coords[i]) # this is not Unitful compatible
    i == j && return zero(T)

    dr = vector(s.coords[i], s.coords[j], s.box_size)
    r2 = sum(abs2, dr)

    if !S && iszero(s.atoms[i].σ) || iszero(s.atoms[j].σ)
        return zero(U)
    end

    σ = sqrt(s.atoms[i].σ * s.atoms[j].σ)
    ϵ = sqrt(s.atoms[i].ϵ * s.atoms[j].ϵ)

    cutoff = inter.cutoff
    σ2 = σ^2
    params = (σ2, ϵ)

    if cutoff_points(C) == 0
        potential(inter, r2, inv(r2), params)
    elseif cutoff_points(C) == 1
        sqdist_cutoff = cutoff.sqdist_cutoff * σ2
        r2 > sqdist_cutoff && return zero(U)

        potential_cutoff(cutoff, r2, inter, params)
    elseif cutoff_points(C) == 2
        r2 > sqdist_cutoff && return zero(U)

        if r2 < activation_dist
            potential(inter, r2, inv(r2), params)
        else
            potential_cutoff(cutoff, r2, inter, params)
        end
    end
end

@fastmath function potential(::SoftSphere, r2, invr2, (σ2, ϵ))
    six_term = (σ2 * invr2) ^ 3

    return 4ϵ * (six_term ^ 2)
end


text:
abstract type AbstractSampler{XT,YT} end

process_loaded_samples!(::AbstractSampler) = nothing

function load_samples!(s::AbstractSampler{<:Any,YT}) where YT
    if !isnothing(s.filename) && isfile(s.filename)
        @info "Loading previous calculation from $(s.filename)"
        data = readdlm(s.filename)
        n = size(data, 1)
        done = convert.(Int, data[:,1])
        data = data[:,2:end]
        if n > 0
            nc = size(data, 2)
            expected = YT <: Complex ? 3 : 2
            nc == expected ||
                throw(ArgumentError("Unexpected number of columns in $(s.filename), expected $(expected)"))
            xn = min(length(s.x), n)
            xn < n &&
                throw(ArgumentError("More data stored ($(n)) than sampler supports ($(xn))"))

            s.x[done] = data[:,1]
            s.y[done] = (YT <: Complex ?
                (data[:,2] + im*data[:,3]) :
                data[:,2])
            s.done[done] .= true
        end
    end
    process_loaded_samples!(s)
end

function save_samples!(s::AbstractSampler{<:Any,<:Real})
    isnothing(s.filename) && return
    sel = done(s)
    writedlm(s.filename, hcat(done(s), s.x[sel], s.y[sel]))
end

function save_samples!(s::AbstractSampler{<:Any,<:Complex})
    isnothing(s.filename) && return
    sel = done(s)
    y = s.y[sel]
    writedlm(s.filename, hcat(done(s), s.x[sel], real(y), imag(y)))
end

struct StaticSampler{XT<:Real, YT, X<:AbstractVector{XT}, Y<:AbstractVector{YT}, Done, Filename} <: AbstractSampler{XT,YT}
    x::X
    y::Y
    done::Done
    filename::Filename
end

function StaticSampler(::Type{YT}, x::AbstractVector, filename) where YT
    y = Vector{YT}(undef, length(x))
    done = falses(length(x))
    StaticSampler(x, y, done, filename)
end

StaticSampler(::Type{YT}, x::Interval, samples::Integer, filename) where YT =
    StaticSampler(YT, range(x, samples), filename)

Base.length(s::StaticSampler) = length(s.x)
done(s::StaticSampler) = findall(s.done)
not_done(s::StaticSampler) = findall(!, s.done)
isdone(s::StaticSampler) = all(s.done)

function Base.enumerate(s::StaticSampler)
    sel = not_done(s)
    zip(sel, view(s.x, sel))
end

get_sample!(s::StaticSampler, i) = s.x[i]

function Base.setindex!(s::StaticSampler, (x,y), i)
    @assert s.x[i] == x
    s.y[i] = y
    s.done[i] = true
end

function Base.map(f::Function, s::AbstractSampler; plot_fun::Union{Function,Nothing}=nothing)
    load_samples!(s)
    nx = length(s)
    for (i,x) in enumerate(s)
        println("Sample: ", i, "/", nx)
        s[i] = (x,f(x))
        save_samples!(s)
        isnothing(plot_fun) || plot_fun(view(s.x, 1:i), s.y)
    end
    isnothing(plot_fun) || plot_fun(s.x, s.y)
end

export StaticSampler


text:
"""
    struct DEposoPapeState{T, U}

An [`AbstractPathState`](@ref) designed for D`Esopo-Pape shortest-path calculations.
"""
struct DEsopoPapeState{T <:Real, U <: Integer} <: AbstractPathState
    parents::Vector{U}
    dists::Vector{T}
end

"""
    desopo_pape_shortest_paths(g, src, distmx=weights(g))

Compute shortest paths between a source `src` and all
other nodes in graph `g` using the [D'Esopo-Pape algorithm](http://web.mit.edu/dimitrib/www/SLF.pdf).
Return a [`LightGraphs.DEsopoPapeState`](@ref) with relevant traversal information.

# Examples
```jldoctest
julia> using LightGraphs

julia> ds = desopo_pape_shortest_paths(CycleGraph(5), 2);

julia> ds.dists
5-element Array{Int64,1}:
 1
 0
 1
 2
 2

julia> ds = desopo_pape_shortest_paths(PathGraph(5), 2);

julia> ds.dists
5-element Array{Int64,1}:
 1
 0
 1
 2
 3
```
"""
function desopo_pape_shortest_paths(g::AbstractGraph,
    src::Integer,
    distmx::AbstractMatrix{T} = weights(g)) where T <: Real
    U = eltype(g)
    nvg = nv(g)
    (src in 1:nvg) || throw(DomainError(src, "src should be in between 1 and $nvg"))
    dists = fill(typemax(T), nvg)
    parents = zeros(U, nvg)
    state = Vector{Int8}()
    state = fill(Int8(2), nvg)
    q = U[src]
    @inbounds dists[src] = 0

    @inbounds while !isempty(q)
        u = popfirst!(q)
        state[u] = 0

        for v in outneighbors(g, u)
            alt = dists[u] + distmx[u, v]
            if (dists[v] > alt)
                dists[v] = alt
                parents[v] = u

                if state[v] == 2
                    state[v] = 1
                    push!(q, v)
                elseif state[v] == 0
                    state[v] = 1
                    pushfirst!(q, v)
                end
            end
        end
    end

    return DEsopoPapeState{T, U}(parents, dists)
end


text:
#Copyright Nathan Parrish 2018
#University of Colorado, Boulder
#Colorado Center for Astrodynamics Research
#Celestial and Spaceflight Mechanics Laboratory
#
#DEMONSTRATION OF TRAJECTORY OPTIMIZATION AND NEURAL NETWORK FITTING:
#   - Initial guess generation
#   - Optimization with direct method (robust)
#   - Optimization with indirect method (sensitive but fast)
#   - Then a neural network is fit to map (δx₀, t) -> (δλ_v(t)).
#   - Neural network is tested, with results plotted at the end of the script.
#
#Note: It is normal for this to take a long time (a few minutes) to run the
#first time. This is because the "DifferentialEquations" and "ForwardDiff"
#packages are very slow to compile. Running subsequent solutions will be many
#times faster.
#
#Note: there are a handful of paths which must be set manually in this script --
#see notes in the code drawing attention to them

#On the first run, import all the packages
if !isdefined(:DifferentialEquations)
    println("Loading packages...")

    #Correct this path for your system:
    push!(LOAD_PATH, "/Users/nathan/Documents/GitHub")

    #Packages by Nathan, from https://github.com/travelingspaceman
    using GeneralCode
    using LowThrustOpt

    #Packages from other sources:
    using MATLAB #used for plotting
    using Plots
    using DifferentialEquations
    using Ipopt
    using JuMP
    using Interpolations
    using ForwardDiff
    plotlyjs()

    println("  Done loading packages.")
end

#Earth-Moon system:
const global MU = 0.012150585609624037 # MU = mu_moon/(mu_moon + mu_planet)
const global DU = 384747.96285603708 #Earth-Moon distance (km) (km per distance unit)
const global TU = 375699.81732246041 #time units (seconds per time unit)
const global r_moon = 1737. #km
const global r_earth = 6378. #km
const global day = 86400. #seconds
##Calculate mu_moon based on MU for consistency with other code:
const global mu_planet = 398600.4415
const global mu_moon = (MU * mu_planet) / (1 - MU)

#Read in initial & final states:
print("Reading endpoints files...")
    path_j = "/Users/nathan/Dropbox/CU/Research/Julia/multipleShooting/"
    X0_states = readdlm(string(path_j, "L2_Anderson_1.txt"))
    X0_times = linspace(0, 1, size(X0_states,2))
    Xf_states = readdlm(string(path_j, "L2_Anderson_2.txt"))
    Xf_times = linspace(0, 1, size(Xf_states,2))
    print("  Endpoints loaded.\n")

########################## initial guess generation ############################
# tof = 30 * day / TU #TU
#     p = 1.
#     rho = 1.
#     mass = 1e3 #kg
#     time_direction = -1
#
#
#     #-- Propagate a single trajectory with thrust:
#     #set up the ODE problem for numerical integration
#     time_direction = 1.
#     println("----------------------------")
#     thrustLimit = 0.1
#     params = (MU, DU, TU, thrustLimit, mass, time_direction, p, rho)
#
#     #Run numerical integration, with optimal control policy:
#     tspan = (0., tof)
#     τ1 = 0.75
#     state_0 = vcat(interpInitialStates(τ1, X0_times, X0_states, MU), zeros(6))
#     prob = ODEProblem(CRTBP_stateCostate_deriv!, state_0, tspan, params)
#     sol_forward = DifferentialEquations.solve(prob, Vern8(), reltol=1e-13, abstol=1e-13)
#
#
#     #### propagate backwards
#     time_direction = -1
#
#     τ2 = 0.5
#     state_f = interpInitialStates(τ2, Xf_times, Xf_states, MU)
#
#     xf = vcat(state_f, zeros(6))
#
#     xf[4:9] *= -1
#     params = (MU, DU, TU, thrustLimit, mass, time_direction, p, rho)
#     prob = ODEProblem(CRTBP_stateCostate_deriv!, xf, tspan, params)
#     sol_backward = DifferentialEquations.solve(prob, Vern8(), reltol=1e-13, abstol=1e-13)

################################################################################

## TRAJECTORY STACKING INITIAL GUESS ###########################################
tof1 = 10 * day / TU #TU (time in first orbit)
    tof2 = 10 * day / TU #TU (time in second orbit)
    tof = tof1 + tof2 #TU
    n_nodes = 30
    t_TU = collect(linspace(0, tof, n_nodes))
    t_TU1 = t_TU[t_TU .<  tof1]
    t_TU2 = t_TU[t_TU .>= tof1]

    p = 1.
    rho = 1.
    mass = 1e3 #kg
    time_direction = -1
    thrustLimit = 0.
    time_direction = 1.

    params = (MU, DU, TU, thrustLimit, mass, time_direction, p, rho)

    #Ballistically propagate first orbit:
    tspan = (0., tof1)
    τ1 = 0.75
    state_0 = vcat(interpInitialStates(τ1, X0_times, X0_states, MU), zeros(6))
    prob = ODEProblem(CRTBP_stateCostate_deriv!, state_0, tspan, params)
    sol_orbit1 = DifferentialEquations.solve(prob, Vern8(), reltol=1e-13, abstol=1e-13)

    #Find closest point on the second orbit:
    τ2₀ = find_τ(Xf_times, Xf_states, sol_orbit1[1:6, end], MU)
    state_f₀ = vcat(interpInitialStates(τ2₀, Xf_times, Xf_states, MU), zeros(6))

    #Ballistically propagate second orbit:
    tspan = (tof1, tof1 + tof2)
    prob = ODEProblem(CRTBP_stateCostate_deriv!, state_f₀, tspan, params)
    sol_orbit2 = DifferentialEquations.solve(prob, Vern8(), reltol=1e-13, abstol=1e-13)

    #Fix final state:
    τ2 = find_τ(Xf_times, Xf_states, sol_orbit2[1:6, end], MU)
    state_f = vcat(interpInitialStates(τ2, Xf_times, Xf_states, MU), zeros(6))

    #Concatenate states
    XC_trajectoryStack = hcat(sol_orbit1(t_TU1)[:,:], sol_orbit2(t_TU2)[:,:] )
    XC_trajectoryStack[:,1] = state_0
    XC_trajectoryStack[:,end] = state_f

    plotlyjs()
    plot(X0_states[1,:], X0_states[2,:], X0_states[3,:], label="Initial orbit")
    plot!(Xf_states[1,:], Xf_states[2,:], Xf_states[3,:], label="Final orbit")
    plot!(XC_trajectoryStack[1,:], XC_trajectoryStack[2,:], XC_trajectoryStack[3,:],
        color=:black, lw=2, label="Guess solution")
    scatter!([state_0[1],], [state_0[2],], [state_0[3],], label="X0", color=:green)
    scatter!([state_f[1],], [state_f[2],], [state_f[3],], label="Xf", color=:red)
    scatter!([1-MU,], [0.,], [0.,], color=:black, label="Moon")

################################################################################

## direct multiple shooting:
t_TU = collect(linspace(0, tof, n_nodes))

    #Random initial guess:
    # X_all = 0.5*randn(6, n_nodes)
    # X_all[1,:] += (1-MU)

    #Trajectory stacking initial guess:
    X_all = XC_trajectoryStack[1:6, :]

    u_all = zeros(3, n_nodes)
    dV1 = zeros(3)
    dV2 = zeros(3)
    Isp = 2000.
    flagEnd = false
    β = 0.
    allowImpulsive = false
    nsteps = 10
    t0_TU = t_TU[1]
    tf_TU = t_TU[end]
    # thrustLimit = 10.
    plot_yn = false
    plotlyjs(legend=false)
    maxIter = 5
    @time (X_all, u_all, τ1, τ2, t_TU, dV1, dV2, defect) =
        multiShoot_CRTBP_direct(X_all, u_all, τ1, τ2, t_TU, dV1, dV2, MU, DU, TU, n_nodes,
        nsteps, mass, Isp, X0_times, X0_states, Xf_times, Xf_states,
        plot_yn, flagEnd, β, allowImpulsive, maxIter)

    state_0 = interpInitialStates(τ1, X0_times, X0_states, MU)
    state_f = interpInitialStates(τ2, Xf_times, Xf_states, MU)

    plotTrajPlotly_direct(X_all, u_all, X0_states, Xf_states, 0.2)
    gui()

## Indirect multiple shooting, with p = 2
#Initial guess from direct method:
XC_all = vcat(X_all, 0.1*randn(6, n_nodes))

    #Fix the final state exactly:
    XC_all[1:6, 1] = state_0[1:6]
    XC_all[1:6, end] = state_f[1:6]

    #With exact derivatives, we get issues when initializing with an exact
    #numerically-integrated initial guess. Adding a tiny amount of noise makes
    #it numerically stable.
    XC_all[:, 2:end-1] += 1e-10 * randn(12, n_nodes-2)

    p = 2.
    thrustLimit = 10. #N ---> For p=2, leave thrust unconstrained. Constrain it later.

    plot_yn = false
    flag_adjointsOnly = true
    plotlyjs(legend=false)
    maxIter = 5
    @time (XC_all, defect, status_flag) = multiShoot_CRTBP_indirect(XC_all, t_TU, MU, DU, TU, n_nodes,
        mass, thrustLimit, plot_yn, flag_adjointsOnly, maxIter, p, rho)

    flag_adjointsOnly = false
    maxIter = 50
    @time (XC_all, defect, status_flag) = multiShoot_CRTBP_indirect(XC_all, t_TU, MU, DU, TU, n_nodes,
        mass, thrustLimit, plot_yn, flag_adjointsOnly, maxIter, p, rho)

    u_all = zeros(3, n_nodes)
    for ind = 1:n_nodes
        u_all[:,ind] = controlLaw_cart(XC_all[10:12,ind], thrustLimit, p, rho, mass)
    end

# "densify" (add more points) for plotting
n_desired = 300
    time_direction = 1
    params = (MU, DU, TU, thrustLimit, mass, time_direction, p, rho)
    (XC_dense, t_dense) = densify(XC_all, t_TU, params, n_desired)
    u_all = zeros(3, size(XC_dense, 2))
    for ind = 1:size(XC_dense, 2)
        u_all[:,ind] = controlLaw_cart(XC_dense[10:12,ind], thrustLimit, p, rho, mass)
    end

#plot control profile
plotlyjs(legend=true, aspect_ratio=:none)
    plot(t_dense*TU/day, norm_many(u_all), xlabel="Time [days]", ylabel="Control [N]",
        lw=2, color=:black, label="Control Magnitude", size=(800, 400))
    plot!(t_dense*TU/day, u_all[1,:], label="x", color=RGB(0,.5, .9))
    plot!(t_dense*TU/day, u_all[2,:], label="y", color=RGB(.9, .4, .4))
    plot!(t_dense*TU/day, u_all[3,:], label="z", color=RGB(0.4, .7, .4))

# Plot nicely
plotlyjs(legend=false)
    plotTrajPlotly_indirect(XC_dense, u_all, X0_states, Xf_states, 0.2)

## Indirect multiple shooting, with p = 1
p = 1.
    t_TU = collect(linspace(0, tof, n_nodes))

    #Fix the final state exactly:
    XC_all[1:6, 1] = state_0[1:6]
    XC_all[1:6, end] = state_f[1:6]

    #With exact derivatives, we get issues when initializing with an exact
    #numerically-integrated initial guess. Adding a tiny amount of noise makes
    #it numerically stable.
    XC_all[:, 2:end-1] += 1e-10 * randn(12, n_nodes-2)

    thrustLimit = 0.05
    println("thrustLimit = $thrustLimit N")
    rho = 1
    maxIter = 30
    flag_adjointsOnly = false
    plot_yn = false
    @time (XC_all, defect, status_flag) = multiShoot_CRTBP_indirect(XC_all, t_TU, MU, DU, TU, n_nodes,
        mass, thrustLimit, plot_yn, flag_adjointsOnly, maxIter, p, rho)

# "densify" (add more points) for plotting
n_desired = 300
    time_direction = 1
    params = (MU, DU, TU, thrustLimit, mass, time_direction, p, rho)
    (XC_dense, t_dense) = densify(XC_all, t_TU, params, n_desired)
    u_all = zeros(3, size(XC_dense, 2))
    for ind = 1:size(XC_dense, 2)
        u_all[:,ind] = controlLaw_cart(XC_dense[10:12,ind], thrustLimit, p, rho, mass)
    end

#plot control profile
plotlyjs(legend=true, aspect_ratio=:none)
    plot(t_dense*TU/day, norm_many(u_all), xlabel="Time [days]", ylabel="Control [N]",
        lw=2, color=:black, label="Control Magnitude", size=(800, 400))
    plot!(t_dense*TU/day, u_all[1,:], label="x", color=RGB(0,.5, .9))
    plot!(t_dense*TU/day, u_all[2,:], label="y", color=RGB(.9, .4, .4))
    plot!(t_dense*TU/day, u_all[3,:], label="z", color=RGB(0.4, .7, .4))
    plot!(title="Control profile")
    gui()

# Plot nicely
plotlyjs(legend=false)
    plotTrajPlotly_indirect(XC_dense, u_all, X0_states, Xf_states, 0.2)

## Reduce rho for sharper thrust on/off:
rho_current = copy(rho)
    rho_target = 1e-4
    plot_yn = false
    @time (XC_all_new, defect, status_flag) = reduceFuel_indirect(XC_all, t_TU, MU, DU, TU, n_nodes, mass,
        thrustLimit, rho_current, rho_target)
    if status_flag == 0
        println("Success!")
        XC_all = copy(XC_all_new)
        rho = copy(rho_target)
    else
        warn("Failed to reduce rho... :/")
    end

# "densify" (add more points) for plotting
n_desired = 300
    time_direction = 1
    params = (MU, DU, TU, thrustLimit, mass, time_direction, p, rho)
    (XC_dense, t_dense) = densify(XC_all, t_TU, params, n_desired)
    u_all = zeros(3, size(XC_dense, 2))
    for ind = 1:size(XC_dense, 2)
        u_all[:,ind] = controlLaw_cart(XC_dense[10:12,ind], thrustLimit, p, rho, mass)
    end

#plot control profile
plotlyjs(legend=true, aspect_ratio=:none)
    plot(t_dense*TU/day, norm_many(u_all), xlabel="Time [days]", ylabel="Control [N]",
        lw=2, color=:black, label="Control Magnitude", size=(800, 400))
    plot!(t_dense*TU/day, u_all[1,:], label="x", color=RGB(0,.5, .9))
    plot!(t_dense*TU/day, u_all[2,:], label="y", color=RGB(.9, .4, .4))
    plot!(t_dense*TU/day, u_all[3,:], label="z", color=RGB(0.4, .7, .4))
    plot!(title="Control profile")
    gui()

# Plot nicely
plotlyjs(legend=false)
    plotTrajPlotly_indirect(XC_dense, u_all, X0_states, Xf_states, 0.2)

println("Forcing exit.")
error("Run the next part manually, if you want to continue...")

## Generate Neural Net training samples ########################################

σ_r = 1e2/DU #Position 1-σ errors (DU)
σ_v = 1e-3 * TU/DU #Velocity 1-σ errors (DU/TU)

N_train = 200
    m = size(XC_all,2)

    #Normally distributed:
    δX0_all = randn(6, N_train)
    δX0_all[1:3, :] *= σ_r #Position, km
    δX0_all[4:6, :] *= σ_v #velocity, km/s

    δλ_all = zeros(6, N_train * m )

    t_all = zeros(N_train * m);

    NN_input = zeros(7, N_train * m)
    NN_output = zeros(6, N_train * m);

for ind = 1:N_train
    println("ind = $ind")

    #perturb the initial state:
    XC_all_mod = copy(XC_all)
    XC_all_mod[1:6,1] += δX0_all[:,ind]

    #add a tiny amount of noise, which can help convergence in some cases
    XC_all_mod[:, 2:end-1] += 1e-10 * randn(12, n_nodes-2)

    #Optimize transfer with perturbed initial condition:
    plot_yn = false
    (XC_all_mod, defect, status_flag) = reduceFuel_indirect(XC_all_mod, t_TU,
        MU, DU, TU, n_nodes, mass, thrustLimit, rho, rho)

    if status_flag == 0 #accept the results
        NN_input[:, (ind-1)* m + 1 : ind * m] = vcat(repmat(δX0_all[:,ind], 1, m), t_TU' )
        NN_output[:, (ind-1)* m + 1 : ind * m] = XC_all_mod[7:12, :] - XC_all[7:12, :]

    else #don't accept the results
        NN_input[:, (ind-1)* m + 1 : ind * m] = NaN * ones(7, m)
        NN_output[:, (ind-1)* m + 1 : ind * m] = NaN * ones(6, m)
    end
end

## TRAIN NEURAL NET ############################################################

#Remove any NaN's (un-converged cases)
temp = .!isnan.(NN_input[1,:])
    NN_input = NN_input[:, temp]
    NN_output = NN_output[:, temp]

    #use δλ_v only:
    NN_output = NN_output[4:6,:]

println("---------------------------------------------------------")
    println("Training NN in MATLAB...")
    warn("Make sure to correct the paths in MATLAB function and in this block of code!")

    NN_size = [10, 10, 10]

    @mput NN_size #pass into MATLAB this way to avoid weird issue with the Julia/MATLAB interface

    #Correct this path for your system:
    mat"addpath('/Users/nathan/Dropbox/CU/Research/ai/')"

    identifier = "L2L2_t0"
    fcn_name = string("NN_fcn_", identifier)

    #Call MATLAB to train network:
    scaleShift = mat"TrainNN_CreateFcn($NN_input, $NN_output, NN_size, $fcn_name)"

    #The training function creates the NN function, but we need to save the scale/shift parameters separately:
    #Correct this path for your system:
    writedlm(string("/Users/nathan/Dropbox/CU/Research/ai/NN_scaleShift_", identifier, ".txt"), scaleShift)

## CALL TRAINED NEURAL NETWORK WITH MATLAB #####################################

#first, read in the scaling parameters
println("---------------------------------------------------------")
    println("Evaluating NN on new test cases...")

    #Correct this path for your system:
    scaleShift = readdlm(string("/Users/nathan/Dropbox/CU/Research/ai/NN_scaleShift_", identifier, ".txt"))
    scale_in = scaleShift[:,1]
    shift_in = scaleShift[:,2]
    scale_out = scaleShift[:,3]
    shift_out = scaleShift[:,4]

    #remove the NaN padding:
    scale_in = scale_in[.!isnan.(scale_in)]
    shift_in = shift_in[.!isnan.(shift_in)]
    scale_out = scale_out[.!isnan.(scale_out)]
    shift_out = shift_out[.!isnan.(shift_out)]

N_val = 100
    δX0_val = randn(6, N_val)
    # δX0_val = zeros(6, N_val)
    δX0_val[1:3, :] *= σ_r #Position, DU
    δX0_val[4:6, :] *= σ_v #velocity, DU/TU

    err_xf_NoCorrect = zeros(6, N_val);
    err_xf_YesCorrect = zeros(6, N_val);

    err_r_NoCorrect = zeros(n_nodes, N_val)
    err_v_NoCorrect = zeros(n_nodes, N_val)

    err_r_YesCorrect = zeros(n_nodes, N_val)
    err_v_YesCorrect = zeros(n_nodes, N_val);

#Need to evaluate the control law at 1,000 times (nom + NN update) to build an
#interpolation object for control. So first, evaluate the nominal at 1,000 times
N_dense = 1000
    time_direction = 1
    params = (MU, DU, TU, thrustLimit, mass, time_direction, p, rho)
    (XC_dense_nom, t_dense) = densify(XC_all, t_TU, params, N_dense)

    nstate = 6
    λv_mag_nom = norm_many(XC_dense_nom[10:12,:]);

#Set up plot for trajectories compare:
plotlyjs(legend=false, aspect_ratio=:1)
    plotTrajPlotly_indirect(XC_dense, zeros(3, size(XC_dense,2)), X0_states, Xf_states, 0.1)
    scatter!([XC_dense[1,1],], [XC_dense[2,1],], [XC_dense[3,1],], marker=2, color=:green)
    scatter!([XC_dense[1,end],], [XC_dense[2,end],], [XC_dense[3,end],], marker=2, color=:red, camera=(45, 20))

#propagate validation samples forward in time to check error at final time.
for ind = 1:1:N_val
    #Evaluate neural network for update to λ(t)
    NN_input_val = vcat( repmat(δX0_val[:,ind], 1, N_dense), t_dense[:]') #position and velocity errors
    NN_input_val_scaled = (NN_input_val - repmat(shift_in, 1, N_dense)) ./ repmat(scale_in, 1, N_dense)

    NN_output_scaled = mat"NN_fcn_L2L2_t0($NN_input_val_scaled)"

    NN_output_val = NN_output_scaled .* repmat(scale_out, 1, N_dense) + repmat(shift_out, 1, N_dense)

    X0_val = XC_all[1:6, 1] + δX0_val[:,ind]

    ## ----- WITHOUT NN CORRECTION -----
    λv_all = XC_dense_nom[10:12,:] # --- NO UPDATE
    λv_itp = scale(interpolate(λv_all', (BSpline(Cubic(Natural())), NoInterp()), OnGrid()), t_dense, 1:3)

    params = (MU, DU, TU, Isp, time_direction, λv_itp, thrustLimit, p, rho)
    tspan = (t_TU[1], t_TU[end])
    prob = ODEProblem(CRTBP_prop_EP_NNControl_deriv!, X0_val, tspan, params)
    sol = DifferentialEquations.solve(prob, Tsit5(), reltol=1e-13, abstol=1e-13)
    plot!(sol[1,:], sol[2,:], sol[3,:], color=:red, label="No correction")

    #calculate final state error
    err_xf_NoCorrect[:, ind] = sol[1:6, end] - XC_all[1:6, end]

    #calculate error along path
    err_r_NoCorrect[:, ind] = norm_many(sol(t_TU)[1:3,:] - XC_all[1:3,:])
    err_v_NoCorrect[:, ind] = norm_many(sol(t_TU)[4:6,:] - XC_all[4:6,:])

    ## ----- WITH NN CORRECTION -----
    λv_all = XC_dense_nom[10:12,:] + NN_output_val # -- YES UPDATE

    ## Truth update (for debugging)
    # XC_all_mod = copy(XC_all)
    # XC_all_mod[1:6,1] += δX0_val[:,ind]
    # (XC_all_mod, defect, status_flag) = reduceFuel_indirect(XC_all_mod, t_TU, MU, DU, TU, n_nodes, mass,
    #     thrustLimit, rho, rho)
    # #densify truth:
    # params = (MU, DU, TU, thrustLimit, mass, time_direction, p, rho)
    # (XC_dense_truth, t_dense_truth) = densify(XC_all_mod, t_TU, params, N_dense)
    # λv_all = XC_dense_truth[10:12,:] #truth

    λv_itp = scale(interpolate(λv_all', (BSpline(Cubic(Natural())), NoInterp()), OnGrid()), t_dense, 1:3)

    params = (MU, DU, TU, Isp, time_direction, λv_itp, thrustLimit, p, rho)
    tspan = (t_TU[1], t_TU[end])
    prob = ODEProblem(CRTBP_prop_EP_NNControl_deriv!, X0_val, tspan, params)
    sol = DifferentialEquations.solve(prob, Tsit5(), reltol=1e-13, abstol=1e-13)
    plot!(sol[1,:], sol[2,:], sol[3,:], color=:green, label="Corrected")

    #calculate final state error
    err_xf_YesCorrect[:, ind] = sol[1:6, end] - XC_all[1:6, end]

    #calculate error along path
    err_r_YesCorrect[:, ind] = norm_many(sol(t_TU)[1:3,:] - XC_all[1:3,:])
    err_v_YesCorrect[:, ind] = norm_many(sol(t_TU)[4:6,:] - XC_all[4:6,:]);

    println("ind = $ind,  err = ", norm(err_xf_YesCorrect[1:3, ind]) * DU, " km")
end
gui()

max_stars_repo_name: travelingspaceman/LowThrustOpt

max_issues_repo_name: travelingspaceman/LowThrustOpt

max_forks_repo_name: travelingspaceman/LowThrustOpt


text:
function max_digit_sum(limit)
    max_sum = 0

    for a in 1:limit - 1
        for b in 1:limit - 1
            digit_sum = sum(digits(BigInt(a) ^ b))
            if digit_sum > max_sum
                max_sum = digit_sum
            end
        end
    end

    return max_sum
end

println(max_digit_sum(100))

max_stars_repo_name: HarrisonGreen/Project-Euler-Solutions

max_issues_repo_name: HarrisonGreen/Project-Euler-Solutions

max_forks_repo_name: HarrisonGreen/Project-Euler-Solutions


text:
module NNPACK

using Libdl, Requires

include("libnnpack_types.jl")
include("error.jl")
include("libnnpack.jl")
include("utils.jl")

export conv, ∇conv_data, ∇conv_filter, softmax, maxpool

const depsjl_path = joinpath(dirname(@__FILE__), "..", "deps", "deps.jl")
if !isfile(depsjl_path)
    error("NNPACK not installed properly, run Pkg.build(\"NNPACK\"), restart Julia and try again")
end
include(depsjl_path)

const interface_path = joinpath(dirname(@__FILE__), "interface.jl")
const full_interface_path = joinpath(dirname(@__FILE__), "interface_full.jl")
global shared_threadpool = Ref(C_NULL)

@init begin
    check_deps()
    status = nnp_initialize()
    status == nnp_status_unsupported_hardware && error("HARDWARE is unsupported by NNPACK")
    try
        global NNPACK_CPU_THREADS = parse(UInt64, ENV["NNPACK_CPU_THREADS"])
    catch
        global NNPACK_CPU_THREADS = 4
    end
    include(interface_path)
    try
        global NNPACK_FAST_OPS = parse(UInt64, ENV["NNPACK_FAST_OPS"])
    catch
        global NNPACK_FAST_OPS = 1
    end
    NNPACK_FAST_OPS == 1 && include(full_interface_path)
    global shared_threadpool = Ref(pthreadpool_create(NNPACK_CPU_THREADS))
end

end


text:
# Show methods
show(io::IO, d::AbstractDocument) = print(io, "A $(typeof(d))")

show(io::IO, crps::Corpus) = print(io, "A Corpus with $(length(crps.documents)) documents")

show(io::IO, dtm::DocumentTermMatrix{T}) where T = begin
    p, n = size(dtm.dtm)
    println(io, "A $(p)x$(n) DocumentTermMatrix{$T}")
end

show(io::IO, coom::CooMatrix{T}) where T = begin
    n, p = size(coom.coom)
    println(io, "A $(n)x$(p) CooMatrix{$T}")
end

titleize(str::AbstractString) = begin
    join(map(uppercasefirst, strip.(split(str, "."))), ". ","")
end

show(io::IO, md::DocumentMetadata) = begin
    _id = ifelse(isempty(md.id), "<no ID>", md.id)
    _name = ifelse(isempty(md.name), "<no name>", "\"$(titleize(md.name))\"")
    _author = ifelse(isempty(md.author), "<unknown author>", "by $(titlecase(md.author))")
    _edition_year = ifelse(isempty(md.edition_year), "?", md.edition_year)
    _published_year = ifelse(isempty(md.published_year), "?", md.published_year)
    printstyled(io, "$(_id)", bold=true)
    printstyled(io, " $(_name)",
                    " $(_author)",
                    " $(_edition_year)",
                    " ($(_published_year))")
end

# Summary methods
"""
    summary(doc)

Shows information about the document `doc`.
"""
function summary(d::AbstractDocument)
    o = ""
    o *= "A $(typeof(d))\n"
    o *= " * Language: $(language(d))\n"
    o *= " * Name: $(name(d))\n"
    o *= " * Author: $(author(d))\n"
    o *= " * Timestamp: $(timestamp(d))\n"
    td = typeof(d)
    if td<:TokenDocument|| td<:NGramDocument
        o *= " * Snippet: ***SAMPLE TEXT NOT AVAILABLE***"
    else
        l = length(text(d))
        sl = clamp(50, 0, l)
        dots = ifelse(sl<l, "...","")
        sample_text = replace(text(d)[1:sl], r"\s+"=>" ")
        o *= " * Snippet: \"$(sample_text)$(dots)\""
    end
    return o
end

"""
    summary(crps)

Shows information about the corpus `crps`.
"""
function summary(crps::Corpus)
    n = length(crps.documents)
    n_s = sum(map(d -> typeof(d)<:StringDocument, crps.documents))
    n_f = sum(map(d -> typeof(d)<:FileDocument, crps.documents))
    n_t = sum(map(d -> typeof(d)<:TokenDocument, crps.documents))
    n_ng = sum(map(d -> typeof(d)<:NGramDocument, crps.documents))
    o = ""
    o *= "A Corpus with $n documents:\n"
    o *= " * $n_s StringDocument's\n"
    o *= " * $n_f FileDocument's\n"
    o *= " * $n_t TokenDocument's\n"
    o *= " * $n_ng NGramDocument's\n\n"
    o *= "Corpus's lexicon contains $(lexicon_size(crps)) tokens\n"
    o *= "Corpus's index contains $(index_size(crps)) tokens"
    return o
end


text:
immutable TestStruct
    a::cl.CL_int
    b::cl.CL_float
end

@testset "OpenCL.Buffer" begin
    @testset "OpenCL.Buffer constructors" begin
        for device in cl.devices()

            ctx = cl.Context(device)
            testarray = zeros(Float32, 1000)

            @test cl.Buffer(Float32, ctx, cl.CL_MEM_ALLOC_HOST_PTR | cl.CL_MEM_READ_ONLY,
                                         length(testarray)) != nothing

            @test cl.Buffer(Float32, ctx, cl.CL_MEM_ALLOC_HOST_PTR | cl.CL_MEM_WRITE_ONLY,
                                         length(testarray)) != nothing

            @test cl.Buffer(Float32, ctx, cl.CL_MEM_ALLOC_HOST_PTR | cl.CL_MEM_READ_WRITE,
                                         length(testarray)) != nothing

            buf = cl.Buffer(Float32, ctx, cl.CL_MEM_ALLOC_HOST_PTR | cl.CL_MEM_READ_WRITE, length(testarray))
            @test length(buf) == length(testarray)

            @test cl.Buffer(Float32, ctx, cl.CL_MEM_COPY_HOST_PTR | cl.CL_MEM_READ_ONLY,
                                         hostbuf=testarray) != nothing

            @test cl.Buffer(Float32, ctx, cl.CL_MEM_COPY_HOST_PTR | cl.CL_MEM_WRITE_ONLY,
                                         hostbuf=testarray) != nothing

            @test cl.Buffer(Float32, ctx, cl.CL_MEM_COPY_HOST_PTR | cl.CL_MEM_READ_WRITE,
                                         hostbuf=testarray) != nothing

            buf = cl.Buffer(Float32, ctx, cl.CL_MEM_COPY_HOST_PTR | cl.CL_MEM_READ_WRITE, hostbuf=testarray)
            @test length(buf) == length(testarray)

            @test cl.Buffer(Float32, ctx, cl.CL_MEM_USE_HOST_PTR | cl.CL_MEM_READ_ONLY,
                                         hostbuf=testarray) != nothing

            @test cl.Buffer(Float32, ctx, cl.CL_MEM_USE_HOST_PTR | cl.CL_MEM_WRITE_ONLY,
                                         hostbuf=testarray) != nothing

            @test cl.Buffer(Float32, ctx, cl.CL_MEM_USE_HOST_PTR | cl.CL_MEM_READ_WRITE,
                                         hostbuf=testarray) != nothing

            buf = cl.Buffer(Float32, ctx, cl.CL_MEM_USE_HOST_PTR | cl.CL_MEM_READ_WRITE, hostbuf=testarray)
            @test length(buf) == length(testarray)

            # invalid buffer size should throw error
            @test_throws cl.CLError cl.Buffer(Float32, ctx, cl.CL_MEM_ALLOC_HOST_PTR, +0)
            @test_throws InexactError cl.Buffer(Float32, ctx, cl.CL_MEM_ALLOC_HOST_PTR, -1)

            # invalid flag combinations should throw error
            @test_throws cl.CLError cl.Buffer(Float32, ctx, cl.CL_MEM_USE_HOST_PTR | cl.CL_MEM_ALLOC_HOST_PTR,
                                             hostbuf=testarray)

            # invalid host pointer should throw error
            @test_throws TypeError cl.Buffer(Float32, ctx, cl.CL_MEM_COPY_HOST_PTR,
                                         hostbuf=C_NULL)

            @test_throws TypeError cl.Buffer(Float32, ctx, cl.CL_MEM_USE_HOST_PTR,
                                         hostbuf=C_NULL)
        end
     end

     @testset "OpenCL.Buffer constructors symbols" begin
         for device in cl.devices()
             ctx = cl.Context(device)

             for mf1 in [:rw, :r, :w]
                 for mf2 in [:copy, :use, :alloc, :null]
                     for mtype in [cl.CL_char,
                                   cl.CL_uchar,
                                   cl.CL_short,
                                   cl.CL_ushort,
                                   cl.CL_int,
                                   cl.CL_uint,
                                   cl.CL_long,
                                   cl.CL_ulong,
                                   cl.CL_half,
                                   cl.CL_float,
                                   cl.CL_double,
                                   #TODO: bool, vector_types, struct_types...
                                   ]
                         testarray = zeros(mtype, 100)
                         if mf2 == :copy || mf2 == :use
                             @test cl.Buffer(mtype, ctx, (mf1, mf2), hostbuf=testarray) != nothing
                             buf = cl.Buffer(mtype, ctx, (mf1, mf2), hostbuf=testarray)
                             @test length(buf) == length(testarray)
                         elseif mf2 == :alloc
                             @test cl.Buffer(mtype, ctx, (mf1, mf2),
                                                          length(testarray)) != nothing
                             buf = cl.Buffer(mtype, ctx, (mf1, mf2), length(testarray))
                             @test length(buf) == length(testarray)
                         end
                     end
                 end
             end

             test_array = Array(TestStruct, 100)
             @test cl.Buffer(TestStruct, ctx, :alloc, length(test_array)) != nothing
             @test cl.Buffer(TestStruct, ctx, :copy, hostbuf=test_array) != nothing

             # invalid buffer size should throw error
             @test_throws cl.CLError cl.Buffer(Float32, ctx, :alloc, +0)
             @test_throws InexactError cl.Buffer(Float32, ctx, :alloc, -1)

             # invalid flag combinations should throw error
             @test_throws ArgumentError cl.Buffer(Float32, ctx, (:use, :alloc), hostbuf=test_array)

             # invalid host pointer should throw error
             @test_throws TypeError cl.Buffer(Float32, ctx, :copy, hostbuf=C_NULL)

             @test_throws TypeError cl.Buffer(Float32, ctx, :use, hostbuf=C_NULL)

         end
     end

     @testset "OpenCL.Buffer fill" begin
        for device in cl.devices()
             if contains(device[:platform][:name], "Portable")
                 # the pocl platform claims to implement v1.2 of the spec, but does not
                 warn("Skipping test OpenCL.Buffer fill for POCL Platform")
                 continue
             end
             ctx = cl.Context(device)
             queue = cl.CmdQueue(ctx)
             testarray = zeros(Float32, 1000)
             buf = cl.Buffer(Float32, ctx, (:rw, :copy), hostbuf=testarray)
             @test length(buf) == length(testarray)

             v = cl.opencl_version(device)
             if v.major == 1 && v.minor < 2
                 platform_name = device[:platform][:name]
                 info("Skipping OpenCL.Buffer fill for $platform_name: fill is a v1.2 command")
                 continue
             end
             cl.fill!(queue, buf, 1f0)
             readback = cl.read(queue, buf)
             @test all(x -> x == 1.0, readback)
             @test all(x -> x == 0.0, testarray)
             @test buf.valid == true
        end
    end

    @testset "OpenCL.Buffer write!" begin
        for device in cl.devices()
            ctx = cl.Context(device)
            queue = cl.CmdQueue(ctx)
            testarray = zeros(Float32, 1000)
            buf = cl.Buffer(Float32, ctx, (:rw, :copy), hostbuf=testarray)
            @test length(buf) == length(testarray)
            cl.write!(queue, buf, ones(Float32, length(testarray)))
            readback = cl.read(queue, buf)
            @test all(x -> x == 1.0, readback) == true
            @test buf.valid == true
        end
    end

    @testset "OpenCL.Buffer empty_like" begin
        for device in cl.devices()
            ctx = cl.Context(device)
            queue = cl.CmdQueue(ctx)
            testarray = zeros(Float32, 1000)
            buf = cl.Buffer(Float32, ctx, (:rw, :copy), hostbuf=testarray)

            @test sizeof(cl.empty_like(ctx, buf)) == sizeof(testarray)
        end
    end

    @testset "OpenCL.Buffer copy!" begin
        for device in cl.devices()
            ctx = cl.Context(device)
            queue = cl.CmdQueue(ctx)
            test_array = fill(2f0, 1000)
            a_buf = cl.Buffer(Float32, ctx, length(test_array))
            b_buf = cl.Buffer(Float32, ctx, length(test_array))
            c_arr = Array(Float32, length(test_array))
            # host to device buffer
            cl.copy!(queue, a_buf, test_array)
            # device buffer to device buffer
            cl.copy!(queue, b_buf, a_buf)
            # device buffer to host
            cl.copy!(queue, c_arr, b_buf)
            @test all(x -> isapprox(x, 2.0), c_arr) == true
        end
    end

    @testset "OpenCL.Buffer map/unmap" begin
        for device in cl.devices()
            ctx = cl.Context(device)
            queue = cl.CmdQueue(ctx)
            b = cl.Buffer(Float32, ctx, :rw, 100)
            for f in (:r, :w, :rw)
                a, evt = cl.enqueue_map_mem(queue, b, f, 0, (10,10))
                cl.wait(evt)
                @test size(a) == (10,10)
                @test typeof(a) == Array{Float32,2}

                # cannot unmap a buffer without same host array
                bad = similar(a)
                @test_throws ArgumentError cl.unmap!(queue, b, bad)

                @test cl.ismapped(b) == true
                cl.unmap!(queue, b, a)
                @test cl.ismapped(b) == false

                # cannot unmap an unmapped buffer
                @test_throws ArgumentError cl.unmap!(queue, b, a)

                # gc here quickly force any memory errors
                Base.gc()
            end
            @test cl.ismapped(b) == false
            a, evt = cl.enqueue_map_mem(queue, b, :rw, 0, (10,10))
            @test cl.ismapped(b) == true
            evt = cl.enqueue_unmap_mem(queue, b, a, wait_for=evt)
            cl.wait(evt)
            @test cl.ismapped(b) == false
        end
    end
end