text:
using Clustering
include("datasets.jl")

function classify(n::Network, X::Array{Float64}, Y::Array{Int64}, pcfg::Dict,
                  train::Bool=true)
    labels = zeros(Int64, length(Y))
    oneh = onehot(Y)
    ma_reward = -Inf
    if ~pcfg["supervised"]
        if train
            n.da[1] = 1.0
        else
            n.da[1] = 0.0
        end
    end
    for sample in eachindex(Y)
        data_in = X[:, sample] * pcfg["fr"]
        ins = rand(n.nin, pcfg["tinput"])
        inputs = falses(n.nin, pcfg["tinput"])
        for i in 1:n.nin
            inputs[i, :] = ins[i, :] .< data_in[mod(i, length(data_in))+1]
        end
        for t in 1:pcfg["tinput"]
            outs = step!(n, inputs[:, t], train)
        end
        ocount = zeros(pcfg["n_classes"])
        for t in 1:pcfg["toutput"]
            outs = step!(n, train)
            for o in 1:length(outs)
                if outs[o]
                    ocount[mod(o, length(ocount))+1] += 1.0
                end
            end
        end
        for t in 1:pcfg["trest"]
            step!(n, train)
        end
        labels[sample] = indmax(ocount)
        ocount /= sum(ocount)
        if pcfg["supervised"] && train
            error = mean(sum((oneh[sample, :] .- ocount).^2))
            reward = 1.0 - error
            if ma_reward == -Inf
                ma_reward = reward
            else
                ma_reward = (1.0 - pcfg["ma_rate"]) * ma_reward + pcfg["ma_rate"] * reward
            end
            dop = 0.0
            if ma_reward != 0.0
                dop = reward / ma_reward
            end
            if dop > 1.0
                n.da[1] += (dop - 1.0)
            end
        end
    end

    mean(Y .== labels), randindex(Y, labels)[1]
end

function run_classify(n::Network, X::Array{Float64}, Y::Array{Int64}, pcfg::Dict)
    init_weights = copy(n.weights)
    trainx, trainy, testx, testy = train_test_split(X, Y)
    static_count = 0

    for i in 1:pcfg["n_epochs"]
        weights = copy(n.weights)
        acc, rindex = classify(n, trainx, trainy, pcfg, true)
        dweight = sum(abs.(n.weights - weights)) / length(n.weights)
        Logging.info(@sprintf("T: %s %d %d %d %d %0.6f %0.6f %e %e",
                              pcfg["data"], pcfg["seed"], pcfg["supervised"],
                              i, static_count, acc, rindex,
                              sum(abs.(n.weights - init_weights)) / length(n.weights),
                              dweight))
        if dweight < (0.001 / length(n.weights))
            static_count += 1
        else
            static_count = 0
        end
        if static_count > 10
            break
        end
    end

    acc, rindex = classify(n, testx, testy, pcfg, false)
    fit = rindex
    if pcfg["supervised"]
        fit = acc
    end
    fit
end


text:
using StokesDiffEq, Parameters
Δx = 1//2^3 # Make this much smaller (1^2-5) for your own tests
mesh = FDMMesh(Δx,mins=[-1;-1],maxs=[1;1])
prob = prob_stokes_dirichletzero
prob = prob_stokes_homogenous

@unpack Δxs,grids,dims,gridSize,square,mins,maxs = mesh
u = zeros(gridSize[1]-1,gridSize[2])
v = zeros(gridSize[1],gridSize[2]-1)
p = zeros(gridSize[1]-1,gridSize[2]-1)
rp = zeros(p)
δq = zeros(p)

#Calculate grids for u,v, and p
ux = grids[1][1:end-1,:]
uy = (grids[2]+Δxs[2]/2)[1:end-1,:]
vx = (grids[1]+Δxs[1]/2)[:,1:end-1]
vy = grids[2][:,1:end-1]
px = grids[1][1:end-1,1:end-1]+Δxs[1]/2
py = grids[2][1:end-1,1:end-1]+Δxs[2]/2
#prob = dirichletzeroStokesExample()
@unpack f₁,f₂,ugD,vgD,uanalytic,vanalytic,panalytic,g,trueknown = prob

u_analytic = float(uanalytic(ux,uy))
vTrue = float(vanalytic(vx,vy))
pTrue = float(panalytic(px,py))

# Impose boundary conditions, cut out ghost points
u[:,1]   = ugD(ux[:,1]  ,uy[:,1])
u[:,end] = ugD(ux[:,end],uy[:,end])
v[1,:]   = vgD(vx[1,:]  ,vy[1,:])
v[end,:] = vgD(vx[end,:],vy[end,:])

#Gauss-Seidel of Velocity Test
err1 = Vector{Float64}(0)
err2 = Vector{Float64}(0)
StokesDiffEq.GSu!(u,f₁,Δxs,pTrue,ugD,grids,ux,uy) # Inplace u => uhalf
StokesDiffEq.GSv!(v,f₂,Δxs,pTrue,vgD,grids,vx,vy) # Inplace v => vhalf
for j = 1:100
  StokesDiffEq.GSu!(u,f₁,Δxs,pTrue,ugD,grids,ux,uy) # Inplace u => uhalf
  StokesDiffEq.GSv!(v,f₂,Δxs,pTrue,vgD,grids,vx,vy) # Inplace v => vhalf
  StokesDiffEq.push!(err1,maximum(u-u_analytic))
  StokesDiffEq.push!(err2,maximum(v-vTrue))
end

#Should Converge
#TEST_PLOT && plot(1:100,[err1 err2],yscale=:log10)

#Now test rp
err1 = Vector{Float64}(0)
val  = Vector{Float64}(0)

StokesDiffEq.calc_rp!(rp,u_analytic,vTrue,Δxs,g,px,py)

#Should be close to zero with true values.
maximum(rp) < 1e-2

ext: jl


text:
module fier

using MKL
using ArgParse
using Base.Threads
using CSV
using DataFrames
using Dates
# using Distributions
using EmpiricalOrthogonalFunctions
using HistogramThresholding
using ImageFiltering
using LinearAlgebra
using Logging
using NCDatasets
using Statistics
using StatsBase
using Polynomials
using Random
using TOML

include("io.jl")
include("reof_fitting.jl")
include("synthesis.jl")
include("fier_interface.jl")

function julia_main()::Cint
    try
        real_main()
    catch
        Base.invokelatest(Base.display_error, Base.catch_stack())
        return 1
    end
    return 0
end

function real_main()
    args = parse_commandline()

    # unpack the command line arguments
    inconfig = args["config"]

    if !isfile(inconfig)
        @error("could not find file $inconfig")
        return
    else
        config = read_config(inconfig)
    end

    # set the BLAS backend to use as many threads as julia process
    BLAS.set_num_threads(nthreads())

    if !args["verbose"]
        Logging.disable_logging(Logging.Info)
    else
        # Run information
        @info "Logging set to info"
        @info "Using $(Threads.nthreads()) threads"
    end

    if args["subcommand"] == "synthesis"
        @info "Running subcommand $(args["subcommand"])"

        synthesize(config)

    elseif args["subcommand"] == "buildregressions"
        @info "Running subcommand $(args["subcommand"])"

        buildregressions(config)

    else
        @error(
            "Could not understand subcommand argument, " *
            "options are 'buildregressions', 'synthesis'"
        )
        return
    end

    return
end

function parse_commandline()
    s = ArgParseSettings(
        description = "Forecasting of Inundation Extents using REOF (FIER) Command Line Application",
        add_version = true,
        version = @project_version
    )

    @add_arg_table s begin
        "subcommand"
        help = "which operation to run, options are 'buildregressions', 'synthesis'"
        required = true

        "--config", "-c"
        arg_type = String
        help = "toml file specifying the configuration for command"
        required = true

        "--verbose"
        help = "flag to set the process to print info"
        action = :store_true

    end

    return parse_args(s)
end

if abspath(PROGRAM_FILE) == @__FILE__
    real_main()
end

end # module

max_stars_repo_name: SERVIR/fier-cli

max_issues_repo_name: SERVIR/fier-cli

max_forks_repo_name: SERVIR/fier-cli


text:
using Blend2D

image = Blend2D.BLImageCore()
@show Blend2D.blImageInit(image)
@show image

codec = Blend2D.BLImageCodecCore()
@show Blend2D.blImageCodecInit(codec)
@show codec

codecs = Blend2D.blImageCodecBuiltInCodecs()
@show codecs
@show Blend2D.blImageCodecFindByName(codec, codecs, "PNG")
@show codec

@show Blend2D.blImageReadFromFile(image, "julia.png", codecs)
@show image

@show Blend2D.blImageWriteToFile(image, "julia_out.png", codec)

data = Blend2D.BLImageData()
@show Blend2D.blImageGetData(image, data)
@show data


text:
isdefined(Base, :__precompile__) && __precompile__()

"""
Module for visualizing julia type hierarchies
"""
module D3TypeTrees

#################################################################################
# Exports

export TypeTree
export DisplayTypeTree

#################################################################################
# Imports

using JSON
using D3Trees
using InteractiveUtils

function wrap(t::Type)
	s = InteractiveUtils.subtypes(t)
	if isempty(s)
		d = try
			[string(fn) for fn in fieldnames(t)]
		catch
			[]
		end
	else
		d = Dict()
		for t_inner in s
			if t_inner != t
				d[string(t_inner)] =  wrap(t_inner)
			end
		end
	end
	return d
end

function get_children(d::Dict, c = Array{Array{Any,1},1}(),n = Array{String,1}(),f = Array{String,1}(); scopesep = ".")
	cc = Array{Any,1}()
	for (k,v) in d
		K = join(split(String(k),"."),scopesep)
		if isa(v,Dict)
			get_children(v,c,n,f,scopesep = scopesep)
			push!(cc,K)
		elseif isa(v,Array)
			push!(c,[])
			push!(cc,K)
		end
		if !(K in n)
			push!(n,K)
			push!(f,isa(v,Array) ? join(v,"\n") : "Abstract")
		end
	end
	push!(c,reverse(cc))
	c = [[findall(x->x==i, reverse(n))[1] for i in j] for j in reverse(c)]
	#c = reverse(c)
	return c[2:end], reverse(n), reverse(f)
end

function TypeTree(t::Type; init_expand=1, scopesep = ".")
	t = Dict(string(t)=>wrap(t))
	children, names, fields = get_children(t,scopesep = scopesep)
	tree = D3Tree(children,text = names,tooltip = fields,init_expand=init_expand)
	return tree
end

function DisplayTypeTree(t::Type; browser="google chrome", init_expand=1, scopesep = ".")
    tree = TypeTree(t,init_expand = init_expand, scopesep = scopesep)
    DisplayTypeTree(tree,browser = browser)
end

function DisplayTypeTree(tree::D3Trees.D3Tree; browser="google chrome")
    inbrowser(tree,browser)
end

end #module


text:
module TestSubtypeMatching
    abstract type GPU end
    struct Context <: GPU end

    struct T{A} end
end

@testset "subtype matching" begin
    function foobar()
        foo()
    end

    @noinline foo() = "foo"
    @noinline bar() = "bar"

    # regular execution, without a context
    @test foobar() == "foo"
    @test TinyCassette.execute(nothing, foobar) == "foo"

    # define a context, as a subtype of some other context
    TinyCassette.execute(ctx::TestSubtypeMatching.GPU, f::typeof(foo)) = bar()
    @test TinyCassette.execute(TestSubtypeMatching.Context(), foobar) == "bar"
end

@testset "subtype matching directly" begin
    foo() = 0
    TinyCassette.execute(ctx::C, f::typeof(foo)) where {C} = 1
    @test TinyCassette.execute(TestSubtypeMatching.Context(), foo) == 1
end

@testset "access type var" begin
    foo() = 0
    TinyCassette.execute(ctx::TestSubtypeMatching.T{A}, f::typeof(foo)) where {A} = A

    @test TinyCassette.execute(TestSubtypeMatching.T{3}(), foo) == 3
    @test TinyCassette.execute(TestSubtypeMatching.T{5}(), foo) == 5
end


text:
"""

The sum of the primes below 10 is 2 + 3 + 5 + 7 = 17.

Find the sum of all the primes below two million.
"""

function primesum(n::Int)
    isprime = trues(n)
    isprime[1] = false

    for i in 2:Int(trunc(√n))
        if isprime[i]
            for j in i^2:i:n
                isprime[j] = false
            end
        end
    end
    return sum(filter(x->isprime[x],1:n))
end

@assert primesum(10) == 17
@show primesum(2000000)


text:
# Estimate derivatives numerically
# These use centered-differencing for higher accuracy
function derivative_numer{T<:Number}(func::Function, x::T, h::T)
    xp = x + h
    vp = func(xp)
    xm = x - h
    vm = func(xm)
    return (vp-vm)/(xp-xm)
end
derivative_numer{T<:Number}(func::Function, x::T, index::Int, h::T) = derivative_numer(func, x, h)  # compatibility with partial derivatives
function derivative_numer{T<:Number}(func::Function, x::Array{T}, index::Int, h::T)
    xsave = x[index]
    xp = xsave + h
    xm = xsave - h
    x[index] = xp
    vp = func(x)
    x[index] = xm
    vm = func(x)
    x[index] = xsave
    return (vp-vm)/(xp-xm)
end
function derivative_numer{T<:Number}(func::Function, c::(T...,), index::Int, h::T)
    x = [c...]
    xsave = x[index]
    xp = xsave + h
    xm = xsave - h
    x[index] = xp
    vp = func(x...)
    x[index] = xm
    vm = func(x...)
    x[index] = xsave
    return (vp-vm)/(xp-xm)
end
function derivative_numer{T<:Number}(func::Function, x::T, h::Vector{T})
    d = zeros(T, length(h))
    for i = 1:length(h)
        d[i] = derivative_numer(func, x, h[i])
    end
    return d
end
function derivative_numer{T<:Number}(func::Function, x, index::Int, h::Vector{T})
    d = zeros(T, length(h))
    for i = 1:length(h)
        d[i] = derivative_numer(func, x, index, h[i])
    end
    return d
end
derivative_numer{T<:Number}(func::Function, x::T) = derivative_numer(func, x, (eps(max(abs(x),one(T))))^convert(T, 1/3))
derivative_numer{T<:Number}(func::Function, x::Array{T}, index::Int) = derivative_numer(func, x, index, (eps(max(abs(x[index]),one(T))))^convert(T, 1/3))
derivative_numer{T<:Number}(func::Function, c::(T...,), index::Int) = derivative_numer(func, c, index, (eps(max(abs(c[index]),one(T))))^convert(T, 1/3))


text:
# MIT License
#
# Copyright (c) 2020 Joshua E Gentry

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.
module DataFiles

using Dates

const INDEX_VERSION = 2
const LANGUAGES_FILE = ".languages"
const PDF_FONT_INFO_FILE = ".pdffont"

const DATA_FILE_EXT = ".traineddata"
const DATE_FORMAT = @dateformat_str("HH:MM:SS dd-mm-yyyy")

const GITHUB_CONTENTS_BASE = "https://api.github.com/repos"

const DATA_URL  = "https://raw.githubusercontent.com/tesseract-ocr/tessdata_best/main/"
const PDF_FONT_URL  = "https://raw.githubusercontent.com/tesseract-ocr/tessconfigs/main/"

include("errors.jl")
include("common.jl")
include("download.jl")
include("github_file.jl")
include("github_project.jl")
include("languages.jl")
include("pdf_font_info.jl")

const DATA_REPO     = GitHubProject("tesseract-ocr", "tessdata_best", "main", "")
const PDF_FONT_REPO = GitHubProject("tesseract-ocr", "tessconfigs", "main", "")

export DATA_FILE_EXT, DATA_URL, PDF_FONT_URL, DATA_REPO, PDF_FONT_REPO
export Languages, GitHubProject, FileError, NetworkError
export load_languages, save_languages, refresh_languages
export has_language, has_all_languages, get_language_file
export file_out_of_date
export load_pdf_font_info, save_pdf_font_info, refresh_pdf_font_info

export download_file

end


text:
using CairoMakie, Random
Random.seed!(145)
n = 50
x = 1:n
y = sin.(x) .+ randn(n)/10
yerr = 0.2*abs.(randn(n))
xerr = abs.(randn(n))
set_theme!(theme_dark())

fig, ax, =  errorbars(x, y, yerr; color = 1:n, colormap = (:white, :red),
    whiskerwidth = 10, linewidth = 1.5,
    figure = (;resolution = (600,400)))
errorbars!(x, y, xerr; color = 1:n, colormap = (:dodgerblue, :white),
    direction = :x, linewidth = 1.5, whiskerwidth = 10)
scatter!(x, y; markersize = 10, color = (:white, 0.25),
    strokewidth = 1, strokecolor = :white)
lines!(x, y; linewidth = 0.5, linestyle = :dashdot)
hidedecorations!(ax; grid = false)
save("errorbars.png", fig, px_per_unit = 2)
display(fig)


text:
# Concrete materials

struct Lambertian{T} <: Material{T}
	albedo::Vec3{T}
end

"""Create a scattered ray emitted by `mat` from incident Ray `r`.

	Args:
		rec: the HitRecord of the surface from which to scatter the ray.

	Return `nothing`` if it's fully absorbed. """
@inline @fastmath function scatter(mat::Lambertian{T}, r::Ray{T}, rec::HitRecord{T})::Scatter{T} where T
	scatter_dir = rec.n⃗ + random_vec3_on_sphere(T)
	if near_zero(scatter_dir) # Catch degenerate scatter direction
		scatter_dir = rec.n⃗
	else
		scatter_dir = normalize(scatter_dir)
	end
	scattered_r = Ray{T}(rec.p, scatter_dir)
	attenuation = mat.albedo
	return Scatter(scattered_r, attenuation)
end

struct Metal{T} <: Material{T}
	albedo::Vec3{T}
	fuzz::T # how big the sphere used to generate fuzzy reflection rays. 0=none
	@inline Metal(a::Vec3{T}, f::T=0.0) where T = new{T}(a,f)
end

@inline @fastmath function scatter(mat::Metal{T}, r_in::Ray{T}, rec::HitRecord)::Scatter{T} where T
	reflected = normalize(reflect(r_in.dir, rec.n⃗) + mat.fuzz*random_vec3_on_sphere(T))
	Scatter(Ray(rec.p, reflected), mat.albedo)
end

"Dielectric, i.e. transparent/refractive material"
struct Dielectric{T} <: Material{T}
	ir::T # index of refraction, i.e. η.
end

@inline @fastmath function scatter(mat::Dielectric{T}, r_in::Ray{T}, rec::HitRecord{T}) where T
	attenuation = SA{T}[1,1,1]
	refraction_ratio = rec.front_face ? (one(T)/mat.ir) : mat.ir # i.e. ηᵢ/ηₜ
	cosθ = min(-r_in.dir⋅rec.n⃗, one(T))
	sinθ = √(one(T) - cosθ^2)
	cannot_refract = refraction_ratio * sinθ > one(T)
	if cannot_refract || reflectance(cosθ, refraction_ratio) > trand(T)
		dir = reflect(r_in.dir, rec.n⃗)
	else
		dir = refract(r_in.dir, rec.n⃗, refraction_ratio)
	end
	Scatter(Ray{T}(rec.p, dir), attenuation) # TODO: rename reflected -> !absorbed?
end


text:
using SensorFeatureTracking, FixedPointNumbers, Colors
using ProfileView

projdir = joinpath(dirname(@__FILE__), "..")
# srcdir = joinpath(projdir,"src")
datadir = joinpath(projdir,"Data")
##
#read images
im1 = load(joinpath(datadir,"testSequence/image_$(70).jpg"))
im2 = load(joinpath(datadir,"testSequence/image_$(71).jpg"))

# feats = getapproxbestharris(im1, 20)
feats = getApproxBestShiTomasi(im1,nfeatures=1000, stepguess=0.95)

flow = BlockTracker(feats)

#Profile---
# flwtest = deepcopy(flow)
# @profile block_tracker!(flwtest, im1, im2) # JIT compiling
# Profile.clear()
# @profile block_tracker!(flwtest, im1, im2)
# ProfileView.view()

oldfeats = deepcopy(flow.features)

@time block_tracker!(flow, im1, im2)

# draw flow lines
image2 = deepcopy(im2)
map((ft1,ft2)->drawfeatureLine!(image2, Feature(ft1),Feature(ft2)),flow.features, oldfeats)
image2

##
# draw and mark features
image1 = deepcopy(im1)
map(ft->drawfeatureX!(image1, Feature(ft), length=3),oldfeats)

image2 = deepcopy(im2)
map(ft->drawfeatureX!(image2, Feature(ft), length=3),flow.features)
map(ft->drawfeatureX!(image2, Feature(ft), length=1),oldfeats)
grid = hcat(image1,image2)
grid

## ---------- Run in loop -------------------------------------------------------
startFrame = 50
numFrames = 50

im1 = load(joinpath(datadir,"testSequence/image_$(startFrame).jpg"))
im2 = load(joinpath(datadir,"testSequence/image_$(startFrame+1).jpg"))

(ro,co) = size(im1)
feats = getApproxBestShiTomasi(im1,nfeatures=200, stepguess=0.95)

flow = BlockTracker(feats, matchFunction = compute_ssd)
# grid_features!(flow,co,ro)

image1 = deepcopy(im1)
oldfeats = deepcopy(flow.features)
# @time block_tracker!(flow, im1, im2)
@time block_tracker!(flow, im1, im2)

blankImg = zeros(Gray{N0f8},ro,co)
map((ft1,ft2)->drawfeatureLine!(blankImg, Feature(ft1),Feature(ft2)),oldfeats,flow.features)

#
tic()
for frame_idx = startFrame+2:startFrame+numFrames -1
    println("Processing Frame: ", frame_idx)
    im1 = deepcopy(im2)
    im2 = load(joinpath(datadir,"testSequence/image_$(frame_idx).jpg"))

    # im2 = im2[range...]
    block_tracker!(flow, im1, im2)

    map((ft1,ft2)->drawfeatureLine!(blankImg, Feature(ft1),Feature(ft2)),flow.features, oldfeats)
    oldfeats = deepcopy(flow.features)
end
toc()

# overlay tracking traces with image
# o_img = (reinterpret(UInt8,blankImg)) .| (reinterpret(UInt8,im2))
# b_img = reinterpret(Gray{N0f8},o_img)
o_img = (reinterpret(UInt8,blankImg)) .| (reinterpret(UInt8,im2))
r_img = reinterpret(Gray{N0f8},o_img)
n_img = .~(reinterpret(UInt8,blankImg)) .& (reinterpret(UInt8,im2))
n_img = reinterpret(Gray{N0f8},n_img)
#
# colorview(RGB, n_img, n_img, b_img)
# colorview(RGB, n_img, g_img, n_img)
colorview(RGB, r_img, n_img, n_img)

max_stars_repo_name: dehann/SensorFeatureTracking

max_issues_repo_name: dehann/SensorFeatureTracking

max_forks_repo_name: dehann/SensorFeatureTracking


text:
include("4.7.1 Sine normalization rules.jl")
include("4.7.2 Tangent normalization rules.jl")
include("4.7.3 Secant normalization rules.jl")
include("4.7.4 (c trig)^m (d trig)^n.jl")
include("4.7.5 Inert trig functions.jl")
include("4.7.6 (c+d x)^m trig(a+b x)^n trig(a+b x)^p.jl")
include("4.7.7 F^(c (a+b x)) trig(d+e x)^n.jl")
include("4.7.8 u trig(a+b log(c x^n))^p.jl")
include("4.7.9 Active trig functions.jl")

max_stars_repo_path: src/rules/4 Trig functions/4.7 Miscellaneous/.jl


text:
###
# Distributed broadcast implementation
##

using Base.Broadcast
import Base.Broadcast: BroadcastStyle, Broadcasted, _max

# We define a custom ArrayStyle here since we need to keep track of
# the fact that it is Distributed and what kind of underlying broadcast behaviour
# we will encounter.
struct DArrayStyle{Style <: BroadcastStyle} <: Broadcast.AbstractArrayStyle{Any} end
DArrayStyle(::S) where {S} = DArrayStyle{S}()
DArrayStyle(::S, ::Val{N}) where {S,N} = DArrayStyle(S(Val(N)))
DArrayStyle(::Val{N}) where N = DArrayStyle{Broadcast.DefaultArrayStyle{N}}()

BroadcastStyle(::Type{<:DArray{<:Any, N, A}}) where {N, A} = DArrayStyle(BroadcastStyle(A), Val(N))

# promotion rules
function BroadcastStyle(::DArrayStyle{AStyle}, ::DArrayStyle{BStyle}) where {AStyle, BStyle}
    DArrayStyle{BroadcastStyle(AStyle, BStyle)}()
end

# # deal with one layer deep lazy arrays
# BroadcastStyle(::Type{<:LinearAlgebra.Transpose{<:Any,T}}) where T <: DArray = BroadcastStyle(T)
# BroadcastStyle(::Type{<:LinearAlgebra.Adjoint{<:Any,T}}) where T <: DArray = BroadcastStyle(T)
# BroadcastStyle(::Type{<:SubArray{<:Any,<:Any,<:T}}) where T <: DArray = BroadcastStyle(T)

# # This Union is a hack. Ideally Base would have a Transpose <: WrappedArray <: AbstractArray
# # and we could define our methods in terms of Union{DArray, WrappedArray{<:Any, <:DArray}}
# const DDestArray = Union{DArray,
#                          LinearAlgebra.Transpose{<:Any,<:DArray},
#                          LinearAlgebra.Adjoint{<:Any,<:DArray},
#                          SubArray{<:Any, <:Any, <:DArray}}
const DDestArray = DArray

# This method is responsible for selection the output type of broadcast
function Base.similar(bc::Broadcasted{<:DArrayStyle{Style}}, ::Type{ElType}) where {Style, ElType}
    DArray(map(length, axes(bc))) do I
        # create fake Broadcasted for underlying ArrayStyle
        bc′ = Broadcasted{Style}(identity, (), map(length, I))
        similar(bc′, ElType)
    end
end

##
# We purposefully only specialise `copyto!`,
# Broadcast implementation that defers to the underlying BroadcastStyle. We can't
# assume that `getindex` is fast, furthermore  we can't assume that the distribution of
# DArray accross workers is equal or that the underlying array type is consistent.
#
# Implementation:
#   - first distribute all arguments
#     - Q: How do decide on the cuts
#   - then localise arguments on each node
##
@inline function Base.copyto!(dest::DDestArray, bc::Broadcasted)
    axes(dest) == axes(bc) || Broadcast.throwdm(axes(dest), axes(bc))

    # Distribute Broadcasted
    # This will turn local AbstractArrays into DArrays
    dbc = bcdistribute(bc)

    asyncmap(procs(dest)) do p
        remotecall_fetch(p) do
            # get the indices for the localpart
            lpidx = localpartindex(dest)
            @assert lpidx != 0
            # create a local version of the broadcast, by constructing views
            # Note: creates copies of the argument
            lbc = bclocal(dbc, dest.indices[lpidx])
            Base.copyto!(localpart(dest), lbc)
            return nothing
        end
    end
    return dest
end

@inline function Base.copy(bc::Broadcasted{<:DArrayStyle})
    dbc = bcdistribute(bc)
    # TODO: teach DArray about axes since this is wrong for OffsetArrays
    DArray(map(length, axes(bc))) do I
        lbc = bclocal(dbc, I)
        copy(lbc)
    end
end

# _bcview creates takes the shapes of a view and the shape of a broadcasted argument,
# and produces the view over that argument that constitutes part of the broadcast
# it is in a sense the inverse of _bcs in Base.Broadcast
_bcview(::Tuple{}, ::Tuple{}) = ()
_bcview(::Tuple{}, view::Tuple) = ()
_bcview(shape::Tuple, ::Tuple{}) = (shape[1], _bcview(tail(shape), ())...)
function _bcview(shape::Tuple, view::Tuple)
    return (_bcview1(shape[1], view[1]), _bcview(tail(shape), tail(view))...)
end

# _bcview1 handles the logic for a single dimension
function _bcview1(a, b)
    if a == 1 || a == 1:1
        return 1:1
    elseif first(a) <= first(b) <= last(a) &&
           first(a) <= last(b)  <= last(b)
        return b
    else
        throw(DimensionMismatch("broadcast view could not be constructed"))
    end
end

# Distribute broadcast
# TODO: How to decide on cuts
@inline bcdistribute(bc::Broadcasted{Style}) where Style = Broadcasted{DArrayStyle{Style}}(bc.f, bcdistribute_args(bc.args), bc.axes)
@inline bcdistribute(bc::Broadcasted{Style}) where Style<:DArrayStyle = Broadcasted{Style}(bc.f, bcdistribute_args(bc.args), bc.axes)

# ask BroadcastStyle to decide if argument is in need of being distributed
bcdistribute(x::T) where T = _bcdistribute(BroadcastStyle(T), x)
_bcdistribute(::DArrayStyle, x) = x
# Don't bother distributing singletons
_bcdistribute(::Broadcast.AbstractArrayStyle{0}, x) = x
_bcdistribute(::Broadcast.AbstractArrayStyle, x) = distribute(x)
_bcdistribute(::Any, x) = x

@inline bcdistribute_args(args::Tuple) = (bcdistribute(args[1]), bcdistribute_args(tail(args))...)
bcdistribute_args(args::Tuple{Any}) = (bcdistribute(args[1]),)
bcdistribute_args(args::Tuple{}) = ()

# dropping axes here since recomputing is easier
@inline bclocal(bc::Broadcasted{DArrayStyle{Style}}, idxs) where Style = Broadcasted{Style}(bc.f, bclocal_args(_bcview(axes(bc), idxs), bc.args))

# bclocal will do a view of the data and the copy it over
# except when the data already is local
function bclocal(x::DArray{T, N, AT}, idxs) where {T, N, AT}
    bcidxs = _bcview(axes(x), idxs)
    makelocal(x, bcidxs...)
end
bclocal(x, idxs) = x

@inline bclocal_args(idxs, args::Tuple) = (bclocal(args[1], idxs), bclocal_args(idxs, tail(args))...)
bclocal_args(idxs, args::Tuple{Any}) = (bclocal(args[1], idxs),)
bclocal_args(idxs, args::Tuple{}) = ()


text:
@static if VERSION < v"0.7-"
    srand(12345)
else
    Random.seed!(12345)
end

@testset "swap!" begin
    for elty in eltypes
        x = rand(elty, n_L1)
        x_cl = cl.CLArray(queue, x)
        y = rand(elty, n_L1)
        y_cl = cl.CLArray(queue, y)
        CLBlast.swap!(length(x_cl), x_cl, 1, y_cl, 1, queue=queue)
        @test cl.to_host(y_cl, queue=queue) ≈ x
        @test cl.to_host(x_cl, queue=queue) ≈ y
    end
end

@testset "scal!" begin
    for elty in eltypes
        x = rand(elty, n_L1)
        x_cl = cl.CLArray(queue, x)
        y_cl = cl.CLArray(queue, x)
        α = rand(elty)
        CLBlast.scal!(length(x_cl), α, x_cl, 1, queue=queue)
        LinearAlgebra.BLAS.scal!(length(x), α, x, 1)
        @test cl.to_host(x_cl, queue=queue) ≈ x

        CLBlast.scal!(α, y_cl, queue=queue)
        @test cl.to_host(y_cl, queue=queue) ≈ x

        for α in (2, 2.f0, 2.0, 2+0im)
            CLBlast.scal!(length(x_cl), α, x_cl, 1, queue=queue)
            x .= α .* x
            @test cl.to_host(x_cl, queue=queue) ≈ x
        end
    end
end

@testset "copy!" begin
    for elty in eltypes
        x = rand(elty, n_L1)
        x_cl = cl.CLArray(queue, x)
        y = rand(elty, n_L1)
        y_cl = cl.CLArray(queue, y)
        CLBlast.copy!(length(x_cl), x_cl, 1, y_cl, 1, queue=queue)
        @test cl.to_host(y_cl, queue=queue) ≈ x
        @test cl.to_host(x_cl, queue=queue) ≈ x
    end
end

@testset "axpy!" begin
    for elty in eltypes
        x = rand(elty, n_L1)
        x_cl = cl.CLArray(queue, x)
        y = rand(elty, n_L1)
        y_cl = cl.CLArray(queue, y)
        α = rand(elty)

        CLBlast.axpy!(length(x_cl), α, x_cl, 1, y_cl, 1, queue=queue)
        y_true = α .* x .+ y
        @test cl.to_host(x_cl, queue=queue) ≈ x
        @test cl.to_host(y_cl, queue=queue) ≈ y_true

        x_cl = cl.CLArray(queue, x)
        y_cl = cl.CLArray(queue, y)
        CLBlast.axpy!(α, x_cl, y_cl, queue=queue)
        @test cl.to_host(x_cl, queue=queue) ≈ x
        @test cl.to_host(y_cl, queue=queue) ≈ y_true

        x_cl = cl.CLArray(queue, x[1:n_L1-1])
        @test_throws DimensionMismatch CLBlast.axpy!(α, x_cl, y_cl, queue=queue)
    end
end

@testset "dot" begin
    for elty in eltypes
        elty <: Real || continue
        x = rand(elty, n_L1)
        x_cl = cl.CLArray(queue, x)
        y = rand(elty, n_L1)
        y_cl = cl.CLArray(queue, y)
        @test LinearAlgebra.BLAS.dot(length(x), x, 1, y, 1) ≈ CLBlast.dot(length(x_cl), x_cl, 1, y_cl, 1, queue=queue)
    end
end

@testset "dotu" begin
    for elty in eltypes
        elty <: Complex || continue
        x = rand(elty, n_L1)
        x_cl = cl.CLArray(queue, x)
        y = rand(elty, n_L1)
        y_cl = cl.CLArray(queue, y)
        @test LinearAlgebra.BLAS.dotu(length(x), x, 1, y, 1) ≈ CLBlast.dotu(length(x_cl), x_cl, 1, y_cl, 1, queue=queue)
    end
end

@testset "dotc" begin
    for elty in eltypes
        elty <: Complex || continue
        x = rand(elty, n_L1)
        x_cl = cl.CLArray(queue, x)
        y = rand(elty, n_L1)
        y_cl = cl.CLArray(queue, y)
        @test LinearAlgebra.BLAS.dotc(length(x), x, 1, y, 1) ≈ CLBlast.dotc(length(x_cl), x_cl, 1, y_cl, 1, queue=queue)
    end
end

@testset "nrm2" begin
    for elty in eltypes
        x = rand(elty, n_L1)
        x_cl = cl.CLArray(queue, x)
        @test LinearAlgebra.BLAS.nrm2(length(x), x, 1) ≈ CLBlast.nrm2(length(x_cl), x_cl, 1, queue=queue)
    end
end

@testset "asum" begin
    for elty in eltypes
        x = rand(elty, n_L1)
        x_cl = cl.CLArray(queue, x)
        @test LinearAlgebra.BLAS.asum(length(x), x, 1) ≈ CLBlast.asum(length(x_cl), x_cl, 1, queue=queue)
    end
end

@testset "sum" begin
    for elty in eltypes
        x = rand(elty, n_L1)
        x_cl = cl.CLArray(queue, x)
        sumx = sum(x)
        @test real(sumx)+imag(sumx) ≈ CLBlast.sum(length(x_cl), x_cl, 1, queue=queue)
    end
end

_internalnorm(z) = abs(real(z)) + abs(imag(z))
@testset "iamax" begin
    for elty in eltypes
        x = rand(elty, n_L1)
        x_cl = cl.CLArray(queue, x)
        #NOTE: +1 due to zero based indexing in OpenCL vs. 1 based indexing in Julia
        idx_CLBlast = CLBlast.iamax(length(x_cl), x_cl, 1, queue=queue) + 1
        idx_BLAS = LinearAlgebra.BLAS.iamax(length(x), x, 1)
        if elty <: Real
            @test _internalnorm(x[idx_BLAS]) ≈ _internalnorm(x[idx_CLBlast])
        else
            @test_broken _internalnorm(x[idx_BLAS]) ≈ _internalnorm(x[idx_CLBlast])
        end
    end
end

@testset "iamin" begin
    for elty in eltypes
        x = rand(elty, n_L1)
        x_cl = cl.CLArray(queue, x)
        #NOTE: +1 due to zero based indexing in OpenCL vs. 1 based indexing in Julia
        idx_CLBlast = CLBlast.iamin(length(x_cl), x_cl, 1, queue=queue) + 1
        @test_broken minimum(_internalnorm, x) ≈ _internalnorm(x[idx_CLBlast])
    end
end

@testset "had!" begin
    for elty in eltypes
        x = rand(elty, n_L1)
        x_cl = cl.CLArray(queue, x)
        y = rand(elty, n_L1)
        y_cl = cl.CLArray(queue, y)
        z = rand(elty, n_L1)
        z_cl = cl.CLArray(queue, z)
        α = rand(elty)
        β = rand(elty)

        CLBlast.had!(length(x_cl), α, x_cl, 1, y_cl, 1, β, z_cl, 1, queue=queue)
        z .= α .* x .* y .+ β .* z
        @test cl.to_host(x_cl, queue=queue) ≈ x
        @test cl.to_host(y_cl, queue=queue) ≈ y
        @test cl.to_host(z_cl, queue=queue) ≈ z
    end
end


text:
module TestMockedTests

using Test
using TestFunctionRunner.Internal: test_functions, test_modules

using ..Samples: SampleMockedTests
using ..Utils: runtests_process, ⊜

const SCRIPT = joinpath(@__DIR__, "../../samples/mockedtests/runtests.jl")

function cleanenv()
    env = Dict{String,String}()
    for (k, v) in ENV
        if startswith(k, "SampleMockedTests.")
            continue
        end
        env[k] = v
    end
    return env
end

function test_runtests_pass()
    if VERSION < v"1.1"
        @test_broken runtests_process(SCRIPT; env = cleanenv()) ⊜ true
        return
    end
    @test runtests_process(SCRIPT; env = cleanenv()) ⊜ true
end

function test_runtests_fail()
    env = cleanenv()
    env["SampleMockedTests.TestNested2.TestNested3.test"] = "false"
    @test runtests_process(SCRIPT; env = env) ⊜ false
end

function test_should_test_module()
    env = cleanenv()
    env["SampleMockedTests.TestShouldTest.test"] = "false"
    @test runtests_process(SCRIPT; env = env) ⊜ false
    env["SampleMockedTests.TestShouldTest.should_test_module"] = "false"
    if VERSION < v"1.1"
        @test_broken runtests_process(SCRIPT; env = env) ⊜ true
        return
    end
    @test runtests_process(SCRIPT; env = env) ⊜ true
end

function test_module_context()
    env = cleanenv()
    env["SampleMockedTests.TestModuleContext.use_module_context"] = "false"
    @test runtests_process(SCRIPT; env = env) ⊜ false
end

function test_discovery()
    @test test_functions(SampleMockedTests) ==
          [SampleMockedTests.test, SampleMockedTests.test_2]
    @test test_modules(SampleMockedTests) == [
        SampleMockedTests.TestModuleContext,
        SampleMockedTests.TestNested,
        SampleMockedTests.TestNested2,
        SampleMockedTests.TestShouldTest,
    ]
end

end  # module


text:
export EndgameTracker, EndgameOptions, track

"""
    AbstractPathTracker

Supertype for path trackers.
"""
abstract type AbstractPathTracker end

Base.broadcastable(T::AbstractPathTracker) = Ref(T)

"""
    EndgameOptions(; options...)

Options controlling the behaviour of a [`EndgameTracker`](@ref).

## Options

* `at_infinity_check = true`: Whether divering paths should be truncated.
* `endgame_start = 0.1`: The point `t` in time where the endgame starts. Set it to `0.0`
  to disable the endgame.
* `only_nonsingular = false`: If `true` don't run the Cauchy endgame to handle singular
  solutions.
* `zero_is_at_infinity = false`: Whether paths going to a solution where at least one
  coordinates is zero should also be considered diverging.

### Parameters
These parameters control the behaviour during the endgame. See [^BT20] for details.

* `max_endgame_steps = 2000`: The maximal number of steps performed during the endgame.
* `max_winding_number = 6`: The maximal winding number which is attempted in the
 Cauchy endgame.
* `min_cond = 1e6`: The minimal condition number after which an endgame strategy is
  considered to be applied.
* `min_cond_growth = 1e4`: The minimal condition number growth after which an
   endgame strategy is considered to be applied.
* `min_coord_growth = 100`: The minimal relative growth of a coordinate necessary to
  to be considered going to infininity (resp. zero).
* `val_at_infinity_tol = 1e-3`: Tolerance on the valuation which has to be
  satisfied before a path is considered to diverge / go to infinity.
* `val_finite_tol = 1e-3`: Tolerance on the valuation which has to be satisfied
  before the Cauchy endgame is started.

[^BT20]: Breiding, P. and Timme, S. "Tropical Endgame", In preparation (2020)
"""
Base.@kwdef mutable struct EndgameOptions
    endgame_start::Float64 = 0.1
    max_endgame_steps::Int = 2000
    max_endgame_extended_steps::Int = 400
    # eg parameter
    min_cond::Float64 = 1e6
    min_cond_growth::Float64 = 1e4
    min_coord_growth::Float64 = 100.0
    zero_is_at_infinity::Bool = false
    at_infinity_check::Bool = true
    only_nonsingular::Bool = false
    # valuation etc
    val_finite_tol::Float64 = 0.05
    val_at_infinity_tol::Float64 = 0.01

    # singular solutions parameters
    max_winding_number::Int = 6
    singular_min_accuracy::Float64 = 1e-6
end

Base.show(io::IO, opts::EndgameOptions) = print_fieldnames(io, opts)
Base.show(io::IO, ::MIME"application/prs.juno.inline", opts::EndgameOptions) = opts

"""
    EndgameTrackerCode

The possible states an `EndgameTracker` can be in:

* `EndgameTrackerCode.tracking`
* `EndgameTrackerCode.success`
* `EndgameTrackerCode.at_infinity`
* `EndgameTrackerCode.at_zero`
* `EndgameTrackerCode.excess_solution`
* `EndgameTrackerCode.post_check_failed`
* `EndgameTrackerCode.polyhedral_failed`
* `EndgameTrackerCode.terminated_accuracy_limit`
* `EndgameTrackerCode.terminated_ill_conditioned`
* `EndgameTrackerCode.terminated_invalid_startvalue`
* `EndgameTrackerCode.terminated_invalid_startvalue_singular_jacobian`
* `EndgameTrackerCode.terminated_max_winding_number`
* `EndgameTrackerCode.terminated_max_steps`
* `EndgameTrackerCode.terminated_step_size_too_small`
"""
module EndgameTrackerCode
using ..TrackerCode: TrackerCode

@enum codes begin
    tracking
    success
    at_infinity
    at_zero
    terminated_accuracy_limit
    terminated_invalid_startvalue
    terminated_invalid_startvalue_singular_jacobian
    terminated_ill_conditioned
    terminated_max_steps
    terminated_max_extended_steps
    terminated_max_winding_number
    terminated_step_size_too_small
    terminated_unknown
    post_check_failed
    excess_solution
    polyhedral_failed
end

function Base.convert(::Type{codes}, code::TrackerCode.codes)
    if code == TrackerCode.success
        return success
    elseif code == TrackerCode.terminated_max_steps
        return terminated_max_steps
    elseif code == TrackerCode.terminated_accuracy_limit
        return terminated_accuracy_limit
    elseif code == TrackerCode.terminated_ill_conditioned
        return terminated_ill_conditioned
    elseif code == TrackerCode.terminated_invalid_startvalue
        return terminated_invalid_startvalue
    elseif code == TrackerCode.terminated_invalid_startvalue_singular_jacobian
        return terminated_invalid_startvalue_singular_jacobian
    elseif code == TrackerCode.terminated_step_size_too_small
        return terminated_step_size_too_small
    elseif code == TrackerCode.terminated_unknown
        return terminated_unknown
    else
        return tracking
    end
end

end

"""
    is_success(code::EndgameTrackerCode.codes)

Returns `true` if `status` indicates a success in tracking.
"""
is_success(code::EndgameTrackerCode.codes) = code == EndgameTrackerCode.success

"""
    is_success(code::EndgameTrackerCode.codes)

Returns `true` if `status` indicates that a path diverged towards infinity.
"""
is_at_infinity(code::EndgameTrackerCode.codes) = code == EndgameTrackerCode.at_infinity

"""
    is_tracking(code::EndgameTrackerCode.codes)

Returns `true` if `status` indicates the tracking is not going on.
"""
is_tracking(code::EndgameTrackerCode.codes) = code == EndgameTrackerCode.tracking

"""
    is_invalid_startvalue(code::EndgameTrackerCode.codes)

Returns `true` if the provided start value was not valid.
"""
is_invalid_startvalue(code::EndgameTrackerCode.codes) =
    code == EndgameTrackerCode.terminated_invalid_startvalue

###########
## State ##
###########

Base.@kwdef mutable struct EndgameTrackerState
    code::EndgameTrackerCode.codes = EndgameTrackerCode.tracking
    # modi
    singular_endgame::Bool = false

    val::Valuation
    winding_number::Union{Nothing,Int} = nothing

    solution::Vector{ComplexF64}
    accuracy::Float64 = NaN
    cond::Float64 = 1.0
    singular::Bool = false
    steps_eg::Int = 0
    ext_steps_eg_start::Int = 0

    jump_to_zero_failed::Tuple{Bool,Bool} = (false, false)
    last_point::Vector{ComplexF64} = copy(solution)
    last_t::Float64 = NaN

    # condition number
    row_scaling::Vector{Float64}
    col_scaling::Vector{Float64} = zeros(length(solution))
    cond_base::Vector{Float64} = copy(row_scaling)

    # at_infinity
    at_infinity_starts::Vector{Float64} = copy(col_scaling)
    at_infinity_tols::Vector{Float64} = copy(col_scaling)
    at_infinity_abs_coords::Vector{Float64} = copy(col_scaling)
    at_infinity_conds::Vector{Float64} = copy(col_scaling)

    # singular endgame
    samples::Vector{TaylorVector{2,ComplexF64}}
    sample_times::Vector{Float64} = zeros(3)
    sample_conds::Vector{Float64} = zeros(3)
    singular_start::Float64 = NaN
    singular_steps::Int = 0
    prediction::Vector{ComplexF64} = copy(solution)
    prev_prediction::Vector{ComplexF64} = copy(solution)
end

EndgameTrackerState(npolynomials::Integer, x::AbstractVector) = EndgameTrackerState(;
    val = Valuation(length(x)),
    solution = copy(x),
    row_scaling = zeros(npolynomials),
    samples = [TaylorVector{2}(ComplexF64, length(x)) for _ = 1:3],
)

"""
    EndgameTracker(tracker::Tracker; options = EndgameOptions())
    EndgameTracker(H::AbstractHomotopy; options = EndgameOptions())

A `EndgameTracker` combines a [`Tracker`](@ref) with an endgame. That is,
while a [`Tracker`](@ref) assumes that the solution path is non-singular and convergent, the endgame
allows to handle singular endpoints as well as diverging paths.
To compute singular solutions the *Cauchy endgame* used, for divering paths a strategy
based on the valuation of local Puiseux series expansion of the path is used.
See [^BT20] for a detailed description.
By convention, a `EndgameTracker` always tracks from ``t=1`` to ``t = 0``.
See [`EndgameOptions`](@ref) for the possible options.

[^BT20]: Breiding, P. and Timme, S. "Tropical Endgame", In preparation (2020)
"""
struct EndgameTracker{H<:AbstractHomotopy,M<:AbstractMatrix{ComplexF64}} <:
       AbstractPathTracker
    tracker::Tracker{H,M}
    state::EndgameTrackerState
    options::EndgameOptions
end

EndgameTracker(H::Homotopy; compile::Union{Bool,Symbol} = COMPILE_DEFAULT[], kwargs...) =
    EndgameTracker(fixed(H; compile = compile); kwargs...)
function EndgameTracker(H::AbstractHomotopy; tracker_options = TrackerOptions(), kwargs...)
    EndgameTracker(Tracker(H; options = tracker_options); kwargs...)
end
function EndgameTracker(tracker::Tracker; options = EndgameOptions())
    if !(options isa EndgameOptions)
        options = EndgameOptions(; options...)
    end
    state = EndgameTrackerState(size(tracker.homotopy, 1), tracker.state.x)
    EndgameTracker(tracker, state, options)
end

Base.show(io::IO, x::EndgameTracker) = print(io, typeof(x), "()")
Base.show(io::IO, ::MIME"application/prs.juno.inline", x::EndgameTracker) = x
Base.broadcastable(T::EndgameTracker) = Ref(T)

function init!(
    endgame_tracker::EndgameTracker,
    x,
    t₁::Real;
    ω::Float64 = NaN,
    μ::Float64 = NaN,
    extended_precision::Bool = false,
)
    @unpack tracker, state, options = endgame_tracker

    tracker.options.min_rel_step_size = 0.0
    init!(tracker, x, t₁, 0.0; ω = ω, μ = μ, extended_precision = extended_precision)
    state.code = status(tracker)
    state.singular_endgame = false
    state.jump_to_zero_failed = (false, false)

    init!(state.val)
    state.winding_number = nothing

    state.solution .= NaN
    state.accuracy = NaN
    state.cond = NaN
    state.singular = false
    state.steps_eg = 0
    state.ext_steps_eg_start = typemax(Int)

    state.row_scaling .= 1
    state.col_scaling .= 1

    state.at_infinity_starts .= NaN

    state.singular_steps = 0

    endgame_tracker
end

function track!(
    endgame_tracker::EndgameTracker,
    x,
    t₁::Real;
    ω::Float64 = NaN,
    μ::Float64 = NaN,
    extended_precision::Bool = false,
    debug::Bool = false,
)
    init!(endgame_tracker, x, t₁; ω = ω, μ = μ, extended_precision = extended_precision)

    while is_tracking(endgame_tracker.state.code)
        step!(endgame_tracker, debug)
    end

    endgame_tracker.state.code
end

# The endgame tracker can be in three different states
# 1) Pre Endgame
# Just forward to a tracker step
#
# 2) Endgame
# Monitor valuation, and decide to possibly switch to singular eg, or cut off path.
# Each step! still corresponds to a single tracker step
#
# 3) Singular Endgame
#
# Here we expect a singular solution at the end. For this the tracker starts a geometric
# progression in which zero is approached.

function step!(endgame_tracker::EndgameTracker, debug::Bool = false)
    @unpack tracker, state, options = endgame_tracker
    # For performance reasons (and to catch degenerate situations) we only allow a maximal
    # number of steps in the endgame. Check if we reached this and if so terminate.
    if state.steps_eg ≥ options.max_endgame_steps
        return (state.code = EndgameTrackerCode.terminated_max_steps)
    elseif ext_steps(tracker.state) - state.ext_steps_eg_start >
           options.max_endgame_extended_steps
        # check if we had previously a singular solution attempt and return this
        if all(!isnan, state.solution) &&
           !isnothing(state.winding_number) &&
           state.accuracy < 1e-5
            state.cond =
                LA.cond(tracker, state.solution, 0.0, state.row_scaling, state.col_scaling)
            state.singular = true
            return (state.code = EndgameTrackerCode.success)
        end
        return (state.code = EndgameTrackerCode.terminated_max_extended_steps)
    end

    state.last_point .= tracker.state.x
    state.last_t = tracker.state.t

    # if we are in the singular_endgame actually, move there
    if state.singular_endgame
        return singular_endgame_step!(endgame_tracker, debug)
    end

    # keep track on whether we wanted to finish tracking / jump to zero
    is_jump_to_zero = iszero(tracker.state.t′)

    # perform a simple tracker step
    step_success = step!(tracker, debug)
    state.code = status(tracker)

    if !is_tracking(state.code)
        tracking_stopped!(endgame_tracker)
        return state.code
    end

    # Keep track on whether we attempted to jump to 0 and it failed.
    # This is indication that we have a solution on a reduced positive dimensional component
    # at the end of the path
    state.jump_to_zero_failed = (last(state.jump_to_zero_failed), is_jump_to_zero)

    # check if we can start the endgame
    t = real(tracker.state.t)
    t ≤ options.endgame_start || return state.code

    if state.steps_eg == 0
        state.ext_steps_eg_start = ext_steps(tracker.state)
    end
    state.steps_eg += 1

    # continue if we didn't make progress
    if !step_success
        return state.code
    end

    # Update the valuation to get more information about the path
    update!(state.val, tracker.predictor, t)
    debug && println(state.val)
    if check_finite!(state, options)
        switch_to_singular!(state, tracker; debug = debug)
        return state.code
    elseif check_at_infinity!(state, tracker, options; debug = debug)
        return state.code
    end

    return state.code
end

function check_finite!(state, options)
    is_finite(
        state.val;
        finite_tol = options.val_finite_tol,
        zero_is_finite = !options.zero_is_at_infinity,
        max_winding_number = options.max_winding_number,
    ) || return false

    m, merr =
        estimate_winding_number(state.val; max_winding_number = options.max_winding_number)
    if merr < options.val_finite_tol
        if m == 1 && !first(state.jump_to_zero_failed)
            return false
        else
            state.winding_number = m
            return true
        end
    else
        return false
    end
end

function check_at_infinity!(state, tracker, options; debug::Bool = false)
    options.at_infinity_check || return false
    # compute new tolerance
    at_infinity_tol!(
        state.at_infinity_tols,
        state.val;
        finite_tol = options.val_finite_tol,
        zero_is_finite = !options.zero_is_at_infinity,
    )
    κ = NaN
    t = real(tracker.state.t)
    debug && @show state.at_infinity_tols
    for (i, tolᵢ) in enumerate(state.at_infinity_tols)
        if tolᵢ < options.val_at_infinity_tol
            # Check if this is the first time the coordinate indicates going to infinity
            if isnan(state.at_infinity_starts[i])
                # check if we need to compute a row and column scaling
                # this is the case if we have no coordinate marked
                if all(isnan, state.at_infinity_starts)
                    state.col_scaling .= weights(tracker.state.norm)
                    row_scaling!(
                        state.row_scaling,
                        workspace(tracker.state.jacobian),
                        state.col_scaling,
                    )
                end
                κ = LA.cond(tracker.state.jacobian, state.row_scaling, state.col_scaling)
                state.at_infinity_conds[i] = κ
                state.at_infinity_abs_coords[i] = fast_abs(tracker.state.x[i])
                state.at_infinity_starts[i] = t

                # The coordinate is already marked as going to infinity, check
                # whether other criteria are also sufficient
            else
                if isnan(κ)
                    κ = LA.cond(
                        tracker.state.jacobian,
                        state.row_scaling,
                        state.col_scaling,
                    )
                end
                vᵢ = state.val.val_x[i]
                at_zero = state.val.val_x[i] > 0
                cond_growth = κ / state.at_infinity_conds[i]
                if at_zero
                    coord_growth =
                        state.at_infinity_abs_coords[i] / fast_abs(tracker.state.x[i])
                else
                    coord_growth =
                        fast_abs(tracker.state.x[i]) / state.at_infinity_abs_coords[i]
                end
                if debug
                    δt = t / state.at_infinity_starts[i]
                    @show i
                    @show coord_growth, 0.25^(4vᵢ)
                    @show κ, cond_growth, δt, log(cond_growth / δt)
                end

                if coord_growth > clamp(0.25^(4vᵢ), 20, options.min_coord_growth) &&
                   (cond_growth > options.min_cond_growth || κ > max(1e8, options.min_cond))
                    if at_zero
                        state.code = EndgameTrackerCode.at_zero
                    else
                        state.code = EndgameTrackerCode.at_infinity
                    end
                    return true
                end
            end
            # Coordinate was marked as going to infinity but it is no more
        elseif !isnan(state.at_infinity_starts[i])
            state.at_infinity_starts[i] = NaN
        end
    end
    false
end

function switch_to_singular!(state, tracker; debug::Bool)
    state.singular_endgame = true
    t = real(tracker.state.t)
    state.singular_start = t
    debug && printstyled("SWITCH TO SINGULAR"; color = :red, bold = true)

    if all(isone, state.row_scaling)
        state.col_scaling .= weights(tracker.state.norm)
        row_scaling!(
            state.row_scaling,
            workspace(tracker.state.jacobian),
            state.col_scaling,
        )
    end
    add_sample!(state, tracker, t)
    state.singular_steps = 0

    winding_number!(tracker.predictor, state.winding_number)
    state.at_infinity_conds[1] = state.sample_conds[1]

    tracker.state.keep_extended_prec = true
    return state
end
function switch_to_regular!(state, tracker)
    state.singular_endgame = false
    winding_number!(tracker.predictor, 1)
    init!(tracker, 0.0)
    return state
end

function singular_endgame_step!(endgame_tracker::EndgameTracker, debug::Bool = false)
    @unpack tracker, state, options = endgame_tracker

    λ = 0.25
    # move in geometric series forward
    t = real(tracker.state.t)
    init!(tracker, λ * t)
    max_steps = false
    while is_tracking(tracker.state.code)
        step!(tracker, debug)

        if (state.steps_eg += 1) ≥ options.max_endgame_steps
            state.code = EndgameTrackerCode.terminated_max_steps
            max_steps = true
            break
        elseif ext_steps(tracker.state) - state.ext_steps_eg_start >
               options.max_endgame_extended_steps
            state.code = EndgameTrackerCode.terminated_max_extended_steps
            max_steps = true
            break
        end
    end
    state.singular_steps += 1
    debug && @show λ * t
    debug && @show max_steps
    if max_steps
        @goto prediction
    end
    if !is_success(status(tracker))
        state.code = status(tracker)
        tracking_stopped!(endgame_tracker)
        return state.code
    end

    # Update the valuation to get more information about the path
    update!(state.val, tracker.predictor, λ * t)
    debug && println(state.val)
    # check that the valuation still is finite and that the winding number estimate
    # matches
    m̂, m̂_err =
        estimate_winding_number(state.val; max_winding_number = options.max_winding_number)
    if m̂ != state.winding_number || m̂_err > 0.1
        # if !check_finite!(state, options)
        debug && printstyled("SWITCH TO REGULAR"; color = :yellow, bold = true)
        switch_to_regular!(state, tracker)
        return state.code
    end

    add_sample!(state, tracker, λ * t)

    state.singular_steps ≥ 2 || return state.code

    # estimate solution at 0
    acc = predict_endpoint!(state, tracker.state.norm)
    if state.singular_steps == 2
        state.accuracy = acc
        state.solution .= state.prediction
        return state.code
    elseif acc < state.accuracy && state.accuracy > 1e-12
        state.accuracy = acc
        state.solution .= state.prediction
        return state.code
    end

    @label prediction
    m = state.winding_number
    κ = state.sample_conds[3]
    zero_cond = 1 / (m + 1)
    for i = 1:length(state.prediction)
        state.solution[i] = (state.val.val_x[i] < zero_cond) * state.prediction[i]
    end
    κ₀ = LA.cond(tracker, state.solution, 0.0, state.row_scaling, state.col_scaling)
    J₀_norm = inf_norm(workspace(tracker.state.jacobian), state.row_scaling)
    if debug
        @show κ κ₀ state.accuracy nanmax(κ₀, inv(J₀_norm))
    end
    if state.accuracy < options.singular_min_accuracy && (
        (
            (m > 1 && κ > options.min_cond && nanmax(κ₀, inv(J₀_norm)) > κ) ||
            (m == 1 && κ₀ > 1e12)
        ) ||
        # # some solution is better than none
        max_steps ||
        # handle univariate case
        (length(state.solution) == 1 && inv(J₀_norm) < options.min_cond)
    )
        state.cond = max(κ₀, inv(J₀_norm))
        state.singular = true
        return (state.code = EndgameTrackerCode.success)
    elseif !max_steps
        switch_to_regular!(state, tracker)
        return state.code
    end

    return state.code
end

function add_sample!(state, tracker, t)
    tx¹ = tracker.predictor.tx¹
    m = state.winding_number
    s = nthroot(t, m)
    # will transform the samples to s-plane
    μ = m * s^(m - 1)
    κ = LA.cond(tracker.state.jacobian, state.row_scaling, state.col_scaling)
    if state.singular_steps ≤ 2
        tyk = state.samples[state.singular_steps+1]
        for i = 1:length(tx¹)
            tyk[i, 1] = tx¹[i, 1]
            tyk[i, 2] = μ * tx¹[i, 2]
        end
        state.sample_times[state.singular_steps+1] = s
        state.sample_conds[state.singular_steps+1] = κ
    else
        # need to shift previous solutions to the left
        ty3 = state.samples[1]
        state.samples[1] = state.samples[2]
        state.sample_times[1] = state.sample_times[2]
        state.sample_conds[1] = state.sample_conds[2]
        state.samples[2] = state.samples[3]
        state.sample_times[2] = state.sample_times[3]
        state.sample_conds[2] = state.sample_conds[3]
        state.samples[3] = ty3
        for i = 1:length(tx¹)
            ty3[i, 1] = tx¹[i, 1]
            ty3[i, 2] = μ * tx¹[i, 2]
        end
        state.sample_times[3] = s
        state.sample_conds[3] = κ
    end
end

function predict_endpoint!(state, norm)
    state.singular_steps ≥ 2 || return Inf
    if state.singular_steps == 2
        cubic_hermite!(
            state.prediction,
            state.samples[1],
            state.sample_times[1],
            state.samples[2],
            state.sample_times[2],
            0.0,
        )
    end
    state.prev_prediction .= state.prediction
    cubic_hermite!(
        state.prediction,
        state.samples[2],
        state.sample_times[2],
        state.samples[3],
        state.sample_times[3],
        0.0,
    )
    p = state.sample_times[3] / state.sample_times[2]
    # Use error estimate from MSW92 (7)
    err = InfNorm()(state.prediction, state.prev_prediction) / abs((p^2)^2 - 1)
    norm_s = InfNorm()(state.prediction)
    if norm_s > 1e-8
        err /= norm_s
    end
    return err
end

function tracking_stopped!(endgame_tracker::EndgameTracker)
    @unpack tracker, state, options = endgame_tracker

    state.accuracy = tracker.state.accuracy
    if is_success(state.code) && state.accuracy > 1e-14
        refine_current_solution!(tracker; min_tol = 1e-14)
    end
    state.solution .= tracker.state.x
    state.winding_number = nothing
    # only update condition number for successfull paths
    if is_success(state.code)
        state.col_scaling .= weights(tracker.state.norm)
        row_scaling!(
            state.row_scaling,
            workspace(tracker.state.jacobian),
            state.col_scaling,
        )
        state.cond =
            LA.cond(tracker, state.solution, 0.0, state.row_scaling, state.col_scaling)
        state.singular = state.cond > 1e14 || state.accuracy > 1e-12
    end
end

#
# """
#     CauchyEndgameResult
#
# An enum indicating the result of the [`cauchy!`](@ref) computation.
#
# # Cases
# * `CAUCHY_SUCCESS`: The endgame was successfull.
# * `CAUCHY_TERMINATED_MAX_WINDING_NUMBER`: The endgame was terminated since the winding
#   number is larger than the provided threshold.
# * `CAUCHY_TERMINATED`: The endgame was terminated due to some other error in the path
#   tracking.
# """
# @enum CauchyEndgameResult begin
#     CAUCHY_SUCCESS
#     CAUCHY_TERMINATED_MAX_WINDING_NUMBER
#     CAUCHY_TERMINATED
# end
#
# """
#     cauchy!(state::EndgameTrackerState, tracker::Tracker, options::EndgameOptions)
#
# Try to predict the value of `x(0)` using the [`CauchyEndgame`](@ref).
# For this we track the polygon defined by ``te^{i2πk/n}`` until we end again at ``x``.
# Here ``n`` is the number of samples we take per loop, `samples_per_loop`.
# The computation gives up if we have a winding number larger than `max_winding_number`.
# It returns a tuple denoting the success ([`CauchyEndgameResult`](@ref)) the computed
# winding number `m::Int` and th expected accuracy of the solution.
#
# [Cauchy's integral formula]: https://en.wikipedia.org/wiki/Cauchy%27s_integral_formula
# """
# function cauchy!(state::EndgameTrackerState, tracker::Tracker, options::EndgameOptions)
#     @unpack last_point, prediction = state
#
#     # winding_number!(tracker.predictor, 1)
#
#     t = real(tracker.state.t)
#     # Mathemtically, we only need `n₀ = ceil(Int, log(eps()) / log(t))` many sample points
#     # to achieve the maximal accuracy since the error is ≈ t^n₀.
#     # However, we have to be careful that we are not getting too close to the singularity
#     # during tracking. E.g. with `n₀ = 3` we track fairly close to the origin for the
#     # first first and third path. So we require at least 8 sample points.
#     n₀ = max(ceil(Int, log(eps()) / log(t)), 8)
#     @unpack x, μ, ω = tracker.state
#
#     # # always use extended precision for cauchy endgame
#     prediction_acc = refine_current_solution!(tracker)
#     # fix tracker to not flip between extended precision and and mach. precision
#     tracker.state.keep_extended_prec = true
#     # disallow hermite predictor
#     tracker.predictor.use_hermite = false
#     tracker.predictor.branch = 0
#
#     state.last_point .= tracker.state.x
#     state.last_t = tracker.state.t
#     prediction .= 0.0
#     sample_point_acc = Inf
#     m = 1
#     Δθ = 2π / n₀
#     result = CAUCHY_TERMINATED_MAX_WINDING_NUMBER
#     while m ≤ options.max_winding_number
#         θⱼ = 0.0
#         tⱼ = complex(t, 0.0)
#         for j = 1:n₀
#             θⱼ += Δθ
#             if j == n₀
#                 # tracker.predictor.branch += 1
#                 tⱼ = complex(t, 0.0)
#             else
#                 tⱼ = t * cis(θⱼ)
#             end
#             # @show tⱼ
#             # @show t_to_s_plane(tⱼ, 3; branch = tracker.predictor.branch)
#             res = track!(tracker, tⱼ; debug = true)
#             sample_point_acc = tracker.state.accuracy
#             prediction_acc = max(prediction_acc, sample_point_acc)
#
#             if !is_success(res)
#                 result = CAUCHY_TERMINATED
#                 @goto _return
#             end
#             prediction .+= x
#         end
#         # Check that loop is closed
#         d = tracker.state.norm(last_point, x)
#         if d < 100 * max(prediction_acc, sample_point_acc)
#             n = n₀ * m
#             prediction .= prediction ./ n
#
#             result = CAUCHY_SUCCESS
#             break
#         end
#         m += 1
#     end
#
#
#     @label _return
#
#     init!(tracker, last_point, t, 0.0; ω = ω, μ = μ, keep_steps = true)
#
#     result, m, 100prediction_acc
# end

"""
    start_parameters!(tracker::EndgameTracker, p)

Set the start parameters of the homotopy of the tracker.
"""
start_parameters!(T::EndgameTracker, p) = (start_parameters!(T.tracker, p); T)

"""
    target_parameters!(tracker::EndgameTracker, p)

Set the target parameters of the homotopy of the tracker.
"""
target_parameters!(T::EndgameTracker, p) = (target_parameters!(T.tracker, p); T)

parameters!(T::EndgameTracker, p, q) = (parameters!(T.tracker, p, q); T)

function solution(endgame_tracker::EndgameTracker)
    get_solution(endgame_tracker.tracker.homotopy, endgame_tracker.state.solution, 0.0)
end

function PathResult(
    endgame_tracker::EndgameTracker,
    start_solution = nothing,
    path_number = nothing,
)
    @unpack tracker, state, options = endgame_tracker
    H = tracker.homotopy

    if is_success(state.code)
        t = 0.0
        solution = get_solution(H, state.solution, 0.0)
        evaluate!(tracker.corrector.r, H, state.solution, complex(0.0))
        residual = LA.norm(tracker.corrector.r, InfNorm())
    else
        t = real(tracker.state.t)
        solution = get_solution(H, tracker.state.x, t)
        evaluate!(tracker.corrector.r, H, tracker.state.x, complex(t))
        residual = LA.norm(tracker.corrector.r, InfNorm())
    end

    PathResult(
        return_code = Symbol(state.code),
        solution = solution,
        t = t,
        singular = state.singular,
        accuracy = state.accuracy,
        residual = residual,
        condition_jacobian = state.cond,
        winding_number = state.winding_number,
        last_path_point = (get_solution(H, state.last_point, state.last_t), state.last_t),
        valuation = t > options.endgame_start ? nothing : copy(state.val.val_x),
        start_solution = start_solution,
        path_number = path_number,
        ω = tracker.state.ω,
        μ = tracker.state.μ,
        extended_precision = tracker.state.extended_prec,
        accepted_steps = tracker.state.accepted_steps,
        rejected_steps = tracker.state.rejected_steps,
        steps_eg = state.steps_eg,
        extended_precision_used = tracker.state.used_extended_prec,
    )
end

"""
    track(endgame_tracker::EndgameTracker, x::AbstractVector, t::Real = 1.0;
          path_number = nothing, debug = false)

Track the given start solution `x` from `t` towards `0` using the given `endgame_tracker`.
Returns a [`PathResult`](@ref).

    track(endgame_tracker::EndgameTracker, r::PathResult, t::Real = 1.0;
          path_number = nothing, debug = false)

Track `solution(r)` from `t` towards `0` using the given `endgame_tracker`.
"""
function track(
    endgame_tracker::EndgameTracker,
    x,
    t₁::Real = 1.0;
    path_number::Union{Nothing,Int} = nothing,
    kwargs...,
)
    track!(endgame_tracker, x, t₁; kwargs...)
    PathResult(endgame_tracker, x, path_number)
end
function track(endgame_tracker::EndgameTracker, r::PathResult, t₁::Real; kwargs...)
    track(endgame_tracker, solution(r), t₁; ω = r.ω, μ = r.μ, kwargs...)
end


text:
using Turing
using Distributions

import Turing: Trace, TraceR, TraceC, current_trace, fork, fork2, randr

function f()
  t = TArray(Int, 1);
  t[1] = 0;
  while true
    rand(current_trace(), Normal(0,1))
    produce(t[1])
    rand(current_trace(), Normal(0,1))
    t[1] = 1 + t[1]
  end
end

# test task copy version of trace
t = TraceC(f)

consume(t); consume(t)
a = fork(t);
consume(a); consume(a)

Base.@assert consume(t) == 2
Base.@assert consume(a) == 4

# test replaying version of trace
t = TraceR(f)

consume(t); consume(t)
a = fork(t);
consume(a); consume(a)

Base.@assert consume(t) == 2
Base.@assert consume(a) == 4

a2 = fork(t)
Base.@assert length(a2.randomness) == 5
Base.@assert t.randomness == a2.randomness
Base.@assert t.index == a2.index


text:
module NonconvexBayesian

export BayesOptOptions, BayesOptAlg

using Reexport, AbstractGPs, IntervalArithmetic, Parameters
using Sobol, ChainRulesCore
@reexport using NonconvexCore
using NonconvexCore: VecModel, AbstractResult, AbstractOptimizer, getinit
using NonconvexCore: VectorOfFunctions, IneqConstraint, EqConstraint
using NonconvexCore: Objective, flatten, @params
import NonconvexCore: optimize!, Workspace
using NonconvexIpopt

include("bayesian.jl")

end


text:
using JuMP, VariationalInequality
using Test

# https://supernet.isenberg.umass.edu/articles/SPE_Model_Information_Asymmetry_in_Quality.pdf
# Problem (15), Data in Table 1, Example 1
m = 2; n = 1

model = VIPModel()

@variable(model, s[i=1:m] >=0)
@variable(model, d[j=1:n] >=0)
@variable(model, Q[i=1:m, j=1:n] >= 0)
@variable(model, q[i=1:m] >= 0)

@constraint(model, supply[i=1:m], s[i] == sum(Q[i,j] for j in 1:n))
@constraint(model, demand[j=1:n], d[j] == sum(Q[i,j] for i in 1:m))

as = [5; 2]
bs = [5; 10]
@mapping(model, pi[i=1:m], as[i] * s[i] + q[i] + bs[i])

ac = [1; 2]
bc = [15; 20]
@mapping(model, c[i=1:m, j=1:n], ac[i,j] * Q[i,j] + bc[i,j] )

ad = [2]
bd = [100]
@NLexpression(model, qhat[j=1:n], sum(q[i]*Q[i,j] for i in 1:m) / ( sum(Q[i,j] for i in 1:m) + 1e-6 ) )
@mapping(model, nrho[j=1:n], ad[j] * d[j] - qhat[j] - bd[j] )

aq = [5; 10]
@mapping(model, Fq[i=1:m], aq[i] * q[i] - pi[i] )

@innerproduct(model, pi, s)
@innerproduct(model, c, Q)
@innerproduct(model, nrho, d)
@innerproduct(model, Fq, q)

for i=1:m, j=1:n
    setvalue(Q[i,j], 1.0)
end

sol1, Fval1, gap1 = solveVIP(model, algorithm=:fixed_point, max_iter=10000, step_size=0.1, tolerance=1e-10)
@assert 0<= gap1 < 1e-6

@show sol1
@show gap1

@show sol1[Q[1,1]]
@show sol1[Q[2,1]]
@show sol1[q[1]]
@show sol1[q[2]]


text:
using PyCall,DataFrames,RCall
unshift!(PyVector(pyimport("sys")["path"]), "")
@pyimport ncdata
len = length

function sorteqn(x,delim)
    r,p = split(string(x),delim)
    r = join(sort(split(r,'+')),'+')
    p = join(sort(split(p,'+')),'+')
    return join([r,p],delim)
end

organics = readtable("./mcm_3_3_1_reactions.csv", separator = ',', header = true)

reactions = map((x) -> vcat([join(matchall(r"[^\h]*",organics[x,4]),""),organics[x,2]]), 1:len(organics[1]))

organics[:reaction]=[sorteqn(i[1],"=") for i in reactions]

groups  =  [match(r"RO2NO\d?|RO2HO2|RO2|NO\d?|HO2|PAN|J",i[2]) for i in reactions]
allocations = [(r"\b(OH)\b","OH"),(r"\b(CH3O|HOCH2CH2O)\b","ALKOXYL"),(r"\b(HCOCO)\b","GLYOXL"),(r"\b(PPN)\b","PAN")]#[[regex,group]]

for i in 1:length(groups)
    groups[i] == nothing ? nothing : continue;
    rxn = split(reactions[i][1],'=')[1]
    groups[i] = match(r"\b(NO\d?|HO2)\b",rxn)

    if groups[i] == nothing

        for aloc in allocations
            if (ismatch( aloc[1] , rxn))
                groups[i] = aloc[2]
                @goto skip
            end
        end

        for j in split(rxn,"+")
            try
                if (ismatch(r"(\[O\-\]\[O\+\]\=C|\[C\]O\[O\]|\[O\]O\[C\])",mcm3[mcm3[:name].==j,:][:smiles][1]) )
                    groups[i] = "CRI"
                    @goto skip
                end
            catch
                groups[i] = "INORGANIC"
                @goto skip
            end
        end
        #print(reactions[i][1],groups[i],"\n")
        groups[i]="DISSOCIATION"
        @label skip
    end
end

function rmrx(i)
    try
        return i.match
    end
    return i
end

organics[Symbol("groups")]= [rmrx(i) for i in groups] #remove RegexMatch objects

allgroups = Set(organics[:groups])

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

eqnmatch = map((x)->ismatch(r"RO2|RO2NO2|RO2NO|RO2HO2|HO2",x),organics[:groups])

organics[Symbol("keep")]=eqnmatch

selection = organics[organics[:keep].==true,:]

matchingspecies = Set(matchall(r"\b\w+\b",join(selection[:reaction]," ")))

filename = "./edwards11.nc"

data = ncdata.get(filename)
specs = names!(DataFrame(data["spec"]),[Symbol(i)for i in data["sc"]])
rates = names!(DataFrame(data["rate"])[9:len(data["rc"])],[Symbol(i)for i in data["rc"][9:len(data["rc"])]])
#rates = rates[:,[Symbol(i) for i in filter(x->!ismatch(r"EMISS|DUMMY",x),data["rc"])]]
specs = specs[:,[Symbol(i) for i in filter(x->!ismatch(r"EMISS|DUMMY",x),data["sc"])]]
rates = rates[:,[Symbol(i) for i in filter(x->sum(rates[x])>0,names(rates))]]
joint = "\n"*(join(names(rates),"\n"))*"\n"

selects = intersect(Set(names(specs)),Set([Symbol(x) for x in matchingspecies]))

rtype = map((x)-> organics[organics[:reaction].== replace(sorteqn(string(x),"-->"),"-->","="),:groups],names(rates))

products= [split(i[4:len(i)],"+") for i in matchall(r"-->([A-z0-9+]*)",joint)]
reactants = [split(i[2:len(i)-1],"+") for i in matchall(r"\n([A-z0-9+]{1,60})([-->]{0,1})",joint)]

node_set = Set([ i[2:len(i)-1] for i in matchall(r"[\n+>-]{1}[A-z0-9]+[\n+>-]{1}",joint) ])
nodes = [i for i in node_set]
dict = Dict( [[i,x] for (x,i) in enumerate(nodes)])

fluxtype = []
flux = []
edges = []
for i in 1:length(reactants)
  rcol=[]
  #println(i)
  for j in reactants[i]
    coeff = match( r"([\d\s\.]*)(\D[\d\D]*)", j)
    dummy = specs[Symbol(coeff[2])]
    try
      push!(rcol,parse(Float32,coeff[1]))
    catch
      push!(rcol,1.0f0)
    end

    prod = 1
    for k in rcol
      prod*=k
    end

    push!(fluxtype, rtype[i])
    push!(flux, prod*rates[i])
    # make edges array
    for l in products[i]
        push!(edges,[match( r"([\d\s\.]*)(\D[\d\D]*)", j)[2],match( r"([\d\s\.]*)(\D[\d\D]*)", l)[2],i, fluxtype[i]])
    end

  end
end

  function smrxn(x,t,spcsel)
    loss = [flux[i[3]][t] for i in x if i[1]==spcsel]
    prod = [flux[i[3]][t] for i in x if i[2]==spcsel]

    #println(prod,loss)

    if len(loss)!=0
      loss = sum(loss)
    else
      loss = 0
    end

    if len(prod)!=0
      prod = sum(prod)
    else
      prod = 0
    end
    return(abs(prod-loss))
  end

  spcsel = "CH3COCH2O2"#"CH3O2"

newedges =  filter((x)-> ((Symbol(x[1]) in selects)&(  Symbol(x[2]) in selects)) , edges)
newedges =  filter((x)-> (((x[1]) == spcsel)|(  (x[2]) == spcsel)) , newedges)

RO2HO2 = filter((x)-> x[4]==AbstractString["RO2HO2"] , newedges)
RO2NO2 = filter((x)-> x[4]==AbstractString["RO2NO2"] , newedges)
RO2NO = filter((x)-> x[4]==AbstractString["RO2NO"] , newedges)
RO2 = filter((x)-> x[4]==AbstractString["RO2"] , newedges)
HO2 = filter((x)-> x[4]==AbstractString["RO2HO2"] , newedges)
NO2 =  filter((x)-> x[4]==AbstractString["NO2"] , newedges)
OH= filter((x)-> x[4]==AbstractString["OH"] , newedges)
NO3= filter((x)-> x[4]==AbstractString["NO3"] , newedges)
J= filter((x)-> x[4]==AbstractString["J"] , newedges)
NO = filter((x)-> x[4]==AbstractString["NO"] , newedges)
RO2NO3= filter((x)-> x[4]==AbstractString["RO2NO3"] , newedges)

file = open("results.csv", "w")
var = Float64(89)
typevar = string(typeof(var), "\n")

write(file,"t,ts,RO2HO2,RO2NO,RO2NO2,RO2,HO2,OH,NO3,J,NO,RO2NO3,concs,\n")
#names = "    ,ts,RO2HO2,RO2NO,RO2NO2,RO2,HO2,OH,NO3,J,NO,RO2NO3,conc        "
for t in 1:144#linspace(2,len(specs[2]),100)
         t = trunc(Int,t)
         write(file, string(t)*",");write(file, string(t)*",")
       newarr = [];push!(newarr,t);push!(newarr,t)
       for i in [RO2HO2,RO2NO,RO2NO2,RO2,HO2,OH,NO3,J,NO,RO2NO3,]
         val=smrxn(i,t,spcsel)
         write(file, string(smrxn(i,t,spcsel))*",")
         println("s2wwww ",t,spcsel)
         #push!(newarr,val)
       end
         #push!(newarr,specs[t,Symbol(specsel)])
          write(file, string(specs[t,Symbol(spcsel)])*",")
         write(file,"\n")
end
#writedlm("results.csv", dmy, ",")
close(file)

print("fini")

max_stars_repo_name: wolfiex/RopaJL

max_issues_repo_name: wolfiex/RopaJL

max_forks_repo_name: wolfiex/RopaJL


text:
#

if VERSION >= v"0.7.0-DEV.1660"
    _maxind(A) = indmax(A).I
else
    _maxind(A) = ind2sub(size(A), indmax(A))
end

function perf_lucompletepivCopy!(A)
    n = size(A, 1)
    rowpiv=zeros(Int, n-1)
    colpiv=zeros(Int, n-1)
    for k = 1:n-1
        As = @compat abs.(A[k:n, k:n])
        μ, λ = _maxind(As)
        μ += k-1; λ += k-1
        rowpiv[k] = μ
        A[[k,μ], 1:n] = A[[μ,k], 1:n]
        colpiv[k] = λ
        A[1:n, [k,λ]] = A[1:n, [λ,k]]
        if A[k,k] ≠ 0
            ρ = k+1:n
            A[ρ, k] = A[ρ, k]/A[k, k]
            A[ρ, ρ] = A[ρ, ρ] - A[ρ, k:k] * A[k:k, ρ]
        end
    end
    return (A, rowpiv, colpiv)
end

function perf_lucompletepivSub!(A)
    n = size(A, 1)
    rowpiv=zeros(Int, n-1)
    colpiv=zeros(Int, n-1)
    for k = 1:n-1
        As = @compat abs.(view(A, k:n, k:n))
        μ, λ = _maxind(As)
        μ += k-1; λ += k-1
        rowpiv[k] = μ
        A[[k,μ], 1:n] = view(A, [μ,k], 1:n)
        colpiv[k] = λ
        A[1:n, [k,λ]] = view(A, 1:n, [λ,k])
        if A[k,k] ≠ 0
            ρ = k+1:n
            A[ρ, k] = view(A, ρ, k)/A[k, k]
            A[ρ, ρ] = view(A, ρ, ρ) - view(A, ρ, k:k) * view(A, k:k, ρ)
        end
    end
    return (A, rowpiv, colpiv)
end


text:
using Manifolds, Manopt, Test
import Random: seed!

include("trust_region_model.jl")

@testset "Riemannian Trust-Region" begin
    seed!(141)
    n = size(A, 1)
    p = 2
    N = Grassmann(n, p)
    M = PowerManifold(N, ArrayPowerRepresentation(), 2)
    x = random_point(M)

    @test_throws ErrorException trust_regions(M, cost, rgrad, rhess, x; ρ_prime=0.3)
    @test_throws ErrorException trust_regions(
        M, cost, rgrad, rhess, x; max_trust_region_radius=-0.1
    )
    @test_throws ErrorException trust_regions(
        M, cost, rgrad, rhess, x; trust_region_radius=-0.1
    )
    @test_throws ErrorException trust_regions(
        M, cost, rgrad, rhess, x; max_trust_region_radius=0.1, trust_region_radius=0.11
    )

    @testset "Allocating Variant" begin
        X = trust_regions(
            M, cost, rgrad, rhess, x; max_trust_region_radius=8.0, debug=[:Stop]
        )
        opt = trust_regions(
            M, cost, rgrad, rhess, x; max_trust_region_radius=8.0, return_options=true
        )
        @test isapprox(M, X, get_solver_result(opt))

        X2 = deepcopy(x)
        trust_regions!(M, cost, rgrad, rhess, X2; max_trust_region_radius=8.0)
        @test isapprox(M, X, X2)
        XuR = trust_regions(
            M, cost, rgrad, rhess, x; max_trust_region_radius=8.0, randomize=true
        )

        @test cost(M, XuR) ≈ cost(M, X)

        XaH = trust_regions(
            M,
            cost,
            rgrad,
            ApproxHessianFiniteDifference(
                M,
                x,
                rgrad;
                steplength=2^(-9),
                vector_transport_method=ProjectionTransport(),
            ),
            x;
            max_trust_region_radius=8.0,
            stopping_criterion=StopAfterIteration(2000) | StopWhenGradientNormLess(1e-6),
        )
        XaH2 = deepcopy(x)
        trust_regions!(
            M,
            cost,
            rgrad,
            ApproxHessianFiniteDifference(
                M,
                x,
                rgrad;
                steplength=2^(-9),
                vector_transport_method=ProjectionTransport(),
            ),
            XaH2;
            stopping_criterion=StopAfterIteration(2000) | StopWhenGradientNormLess(1e-6),
            max_trust_region_radius=8.0,
        )
        @test isapprox(M, XaH, XaH2; atol=1e-6)
        @test cost(M, XaH) ≈ cost(M, X)

        ξ = random_tangent(M, x)
        @test_throws MethodError get_hessian(SubGradientProblem(M, cost, rgrad), x, ξ)

        η = truncated_conjugate_gradient_descent(
            M, cost, rgrad, x, ξ, rhess; trust_region_radius=0.5
        )
        ηOpt = truncated_conjugate_gradient_descent(
            M, cost, rgrad, x, ξ, rhess; trust_region_radius=0.5, return_options=true
        )
        @test get_solver_result(ηOpt) == η
    end
    @testset "Mutating" begin
        g = RGrad(M, A)
        h = RHess(M, A, p)
        x3 = deepcopy(x)
        trust_regions!(
            M,
            cost,
            g,
            h,
            x3;
            max_trust_region_radius=8.0,
            evaluation=MutatingEvaluation(),
            debug=[:Stop],
        )
        x4 = deepcopy(x)
        opt = trust_regions!(
            M,
            cost,
            g,
            h,
            x4;
            max_trust_region_radius=8.0,
            evaluation=MutatingEvaluation(),
            return_options=true,
        )
        @test isapprox(M, x3, x4)
        XaH = deepcopy(x)
        trust_regions!(
            M,
            cost,
            g,
            ApproxHessianFiniteDifference(
                M,
                deepcopy(x),
                g;
                steplength=2^(-9),
                vector_transport_method=ProjectionTransport(),
                evaluation=MutatingEvaluation(),
            ),
            XaH;
            stopping_criterion=StopAfterIteration(2000) | StopWhenGradientNormLess(1e-6),
            max_trust_region_radius=8.0,
            evaluation=MutatingEvaluation(),
        )
        @test cost(M, XaH) ≈ cost(M, x3)
    end
end


text:
"""
    NamedMatch{T,M<:QueryMatch}

This object has two fields,

* `match::M`, which holds a `QueryMatch` object corresponding to the match
* and `metadata::T`, which holds a `NamedTuple` of metadata.

and is created by the method [`match(query::NamedQuery, obj)`](@ref).

`NamedMatch` satisfies the Tables.jl `AbstractRow` interface. This means that
a vector of `NamedMatch` objects is a valid Tables.jl-compatible table.

## Example

```jldoctest
julia> document_1 = Document("one", (; document_name = "a"))
Document with text "one ". Metadata: (document_name = "a",)

julia> document_2 = Document("Two but there's also a one here.", (; document_name = "b"))
Document starting with "Two but there's…". Metadata: (document_name = "b",)

julia> query = NamedQuery(Query("one"), "find one")
NamedQuery
├─ (query_name = "find one",)
└─ Query("one")

julia> matches = match_all(query, Corpus([document_1, document_2], (;corpus_name="corpus")));

julia> using Tables

julia> Tables.istable(matches)
true

julia> Tables.schema(Tables.rowtable(matches))
Tables.Schema:
 :document       Document{NamedTuple{(:document_name,),Tuple{String}}}
 :distance       Int64
 :indices        UnitRange{Int64}
 :query          Query
 :query_name     String
 :corpus_name    String
 :document_name  String

```
"""
struct NamedMatch{T,M<:QueryMatch}
    match::M
    metadata::T
    function NamedMatch(match::M, metadata::T) where {M,T}
        check_keys(T)
        return new{T,M}(match, metadata)
    end
end

function Tables.getcolumn(m::NamedMatch, i::Int)
    if i === 1
        return m.match.query
    elseif i === 2
        return m.match.document
    elseif i === 3
        return m.match.distance
    elseif i === 4
        return m.match.indices
    else
        return m.metadata[i - 4]
    end
end

function Tables.getcolumn(m::NamedMatch, s::Symbol)
    if s === :query
        return m.match.query
    elseif s === :document
        return m.match.document
    elseif s === :distance
        return m.match.distance
    elseif s === :indices
        return m.match.indices
    else
        return m.metadata[s]
    end
end

function Tables.columnnames(m::NamedMatch)
    return (:document, :distance, :indices, :query, keys(m.metadata)...)
end
Tables.isrowtable(::Type{<:AbstractVector{<:NamedMatch}}) = true

explain(io::IO, m::NamedMatch; context=40) = explain(io, m.match; context=context)

struct NamedQuery{T<:NamedTuple,Q<:AbstractQuery} <: AbstractQuery
    query::Q
    metadata::T

    function NamedQuery(query::Q, metadata::T) where {T,Q}
        check_keys(T)
        return new{T,Q}(query, metadata)
    end
end

function NamedQuery(query::AbstractQuery, name::AbstractString)
    return NamedQuery(query, (; query_name=name))
end

# Here, it seems like we could just dispach on `obj::Union{Corpus, Document}`.
# However, that leads to an ambiguity with `match(query::KeywordSearch.AbstractQuery, p::Corpus)`
# in `src/corpus.jl`. We solve that by creating a separate methods for `Document`s and `Corpus`s
# here.
for OT in (:Corpus, :Document)
    @eval function Base.match(Q::NamedQuery, obj::$(OT))
        disjoint_keys_check(Q, obj)
        m = match(Q.query, obj)
        m === nothing && return nothing
        return NamedMatch(m, (; Q.metadata..., obj.metadata..., m.document.metadata...))
    end

    @eval function match_all(Q::NamedQuery, obj::$(OT))
        disjoint_keys_check(Q, obj)
        matches = match_all(Q.query, obj)
        return [NamedMatch(m, (; Q.metadata..., obj.metadata..., m.document.metadata...))
                for m in matches]
    end
end

"""
    NamedQuery(metadata::Union{String, NamedTuple}, query::AbstractQuery)

Creates a `NamedQuery` that stores a metadata field holding information
about the query. When used with `match`, returns a [`NamedMatch`](@ref), which
carries the metadata of the `NamedQuery` as well as the metadata of the
object which was matched.

## Example

```jldoctest
julia> document_1 = Document("One", (; document_name = "a"))
Document with text "One ". Metadata: (document_name = "a",)

julia> document_2 = Document("Two", (; document_name = "b"))
Document with text "Two ". Metadata: (document_name = "b",)

julia> corpus = Corpus([document_1, document_2], (; corpus_name = "Numbers"))
Corpus with 2 documents, each with metadata keys: (:document_name,)
Corpus metadata: (corpus_name = "Numbers",)

julia> query = NamedQuery(FuzzyQuery("one"), "find one")
NamedQuery
├─ (query_name = "find one",)
└─ FuzzyQuery("one", DamerauLevenshtein{Nothing}(nothing), 2)

julia> m = match(query, corpus)
NamedMatch
├─ (query_name = "find one", corpus_name = "Numbers", document_name = "a")
└─ QueryMatch with distance 1 at indices 1:3.
   ├─ FuzzyQuery("one", DamerauLevenshtein{Nothing}(nothing), 2)
   └─ Document with text "One ". Metadata: (document_name = "a",)

julia> m.metadata
(query_name = "find one", corpus_name = "Numbers", document_name = "a")
```
"""
NamedQuery

max_stars_repo_name: beacon-biosignals/KeywordSearch

max_issues_repo_name: beacon-biosignals/KeywordSearch

max_forks_repo_name: beacon-biosignals/KeywordSearch


text:
using Plots

timings = #Data from julia.mit.edu
Dict(
     1 => ( 15.585,  15.721),
     2 => ( 16.325,  15.893),
     3 => ( 32.079,  31.678),
     4 => ( 46.037,  33.257),
     5 => ( 61.947,  51.412),
     6 => ( 77.247,  50.769),
     7 => ( 92.588,  67.796),
     8 => (108.057,  65.185),
     9 => (123.556,  68.197),
    10 => (138.408,  66.658),
    11 => (153.392,  82.580),
    12 => (168.874,  83.373),
    13 => (184.227,  82.875),
    14 => (199.741,  83.098),
    15 => (215.881, 101.271),
    16 => (230.531, 104.860),
    17 => (246.576, 103.229),
    18 => (261.412, 102.497),
    19 => (276.926, 103.367),
    20 => (294.404, 104.162),
    21 => (308.995, 104.944),
    22 => (323.287, 104.838),
    23 => (340.173, 121.495),
    24 => (353.717, 119.729),
    25 => (369.113, 120.281),
    26 => (384.638, 118.710),
    27 => (402.484, 119.237),
    28 => (417.036, 119.980),
    29 => (431.273, 120.463),
    30 => (446.288, 120.560),
    31 => (461.756, 137.209),
    32 => (476.653, 135.873),
    33 => (493.989, 138.150),
    34 => (508.089, 137.998),
    35 => (523.703, 139.059),
    36 => (539.305, 135.201),
    37 => (554.201, 138.242),
    38 => (569.276, 137.470),
    39 => (585.268, 137.261),
    40 => (599.921, 137.871),
    41 => (615.718, 137.887),
    42 => (632.463, 137.579),
    43 => (647.355, 138.181),
    44 => (664.344, 139.917),
    45 => (678.899, 138.382),
    46 => (695.815, 138.354),
    47 => (717.068, 154.826),
    48 => (730.902, 155.603),
    49 => (745.608, 154.549),
    50 => (756.012, 154.244),
    51 => (769.611, 154.000),
    52 => (793.350, 152.893),
    53 => (801.568, 155.276),
    54 => (818.174, 153.436),
    55 => (832.994, 155.923),
    56 => (849.526, 155.721),
    57 => (863.337, 155.289),
    58 => (881.577, 156.188),
    59 => (893.735, 157.889),
    60 => (911.597, 155.557),
    61 => (925.117, 157.048),
    62 => (940.897, 153.338),
    63 => (954.825, 173.942),
    64 => (990.757, 173.236),
    65 => (987.457, 172.005),
    66 => (1009.688, 173.869),
    67 => (1022.357, 173.303),
    68 => (1043.536, 175.044),
    69 => (1052.156, 173.361),
    70 => (1067.174, 172.230),
    71 => (1079.543, 173.491),
    72 => (1101.516, 174.311),
    73 => (1119.993, 174.057),
    74 => (1136.092, 175.804),
    75 => (1154.758, 175.075),
    76 => (1170.732, 173.453),
    77 => (1184.598, 171.936),
    78 => (1190.922, 172.809),
    79 => (1204.645, 172.878),
)

gflops = zeros(79)
gflopsser = zeros(79)
N = 4096
for (k, (t1, t2)) in timings
    gflops[k] = k*(2*N^3)/t2*1e-9
    gflopsser[k] = k*(2*N^3)/t1*1e-9
end

plot(1:79, gflops, xlabel="#cores", ylabel="GFlops", label="parallel prefix")
plot!(1:79, gflopsser, label="serial")
png("../plots/flops-prefix")


text:
# --------------------------------------------------------------------------
# ACE.jl and SHIPs.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# All rights reserved.
# --------------------------------------------------------------------------

##

using ACE, JuLIP, BenchmarkTools, LinearAlgebra, Test, Random, StaticArrays
using JuLIP
using JuLIP: evaluate, evaluate_d
using JuLIP.Testing

function randR()
   R = rand(JVecF) .- 0.5
   return (0.9 + 2 * rand()) * R/norm(R)
end
randR(N) = [ randR() for n=1:N ], zeros(Int16, N), 0
randcoeffs(B) = rand(length(B)) .* (1:length(B)).^(-2)

##
@info("Generate a basis")
basis = ACE.Utils.TestBasis(6,13)
@show length(basis)
@info("Convert to the SHIP")
ship = SHIP(basis, randcoeffs(basis))
@info("Compress the coefficients")
shipc = ACE.compressA(ship)
@info("check the lengths")
length(ship.coeffs[1])
length(shipc.coeffs[1])

@info("Evaluation test")
Rs = ACE.rand_vec(ship.J, 30)
Zs = zeros(Int16, length(Rs))
z0 = 0
tmp = ACE.alloc_temp(ship, length(Rs))
tmp_d = ACE.alloc_temp_d(ship, length(Rs))
dEs = zeros(JVecF, 30)
@info("  evaluate!:")
@btime ACE.evaluate!($tmp, $ship, $Rs, $Zs, $z0)
@btime ACE.evaluate!($tmp, $shipc, $Rs, $Zs, $z0)
@info("  evaluate_d!:")
@btime ACE.evaluate_d!($dEs, $tmp_d, $ship, $Rs, $Zs, $z0)
@btime ACE.evaluate_d!($dEs, $tmp_d, $shipc, $Rs, $Zs, $z0)


text:
# Variables in Gurobi model

# variable kinds

const GRB_CONTINUOUS = convert(Cchar, 'C')
const GRB_BINARY     = convert(Cchar, 'B')
const GRB_INTEGER    = convert(Cchar, 'I')

# add_var!
function add_var!(model::Model, numnz::Integer, vind::Vector, vval::Vector{Float64}, obj::Float64, lb::Float64, ub::Float64, vtype::Cchar)
    ret = @grb_ccall(addvar, Cint, (
        Ptr{Void},    # model
        Cint,         # numnz
        Ptr{Cint},    # vind
        Ptr{Float64}, # vval
        Float64,      # obj
        Float64,      # lb
        Float64,      # ub
        UInt8,        # vtype
        Ptr{UInt8}    # name
        ),
        model, numnz, ivec(vind.-1), vval, obj, lb, ub, vtype, C_NULL)

    if ret != 0
        throw(GurobiError(model.env, ret))
    end
    nothing
end

function add_var!(model::Model, vtype::Cchar, c::Float64, lb::Float64, ub::Float64)
    ret = @grb_ccall(addvar, Cint, (
        Ptr{Void},    # model
        Cint,         # numnz
        Ptr{Cint},    # vind
        Ptr{Float64}, # vval
        Float64,      # obj
        Float64,      # lb
        Float64,      # ub
        UInt8,        # vtype
        Ptr{UInt8}    # name
        ),
        model, 0, C_NULL, C_NULL, c, lb, ub, vtype, C_NULL)

    if ret != 0
        throw(GurobiError(model.env, ret))
    end
    nothing
end

function add_var!(model::Model, vtype::GChars, c::Real, lb::Real, ub::Real)
    add_var!(model, cchar(vtype), Float64(c), Float64(lb), Float64(ub))
end
add_var!(model::Model, vtype::GChars, c::Real) = add_var!(model, vtype, c, -Inf, Inf)

add_cvar!(model::Model, c::Real, lb::Real, ub::Real) = add_var!(model, GRB_CONTINUOUS, c, lb, ub)
add_cvar!(model::Model, c::Real) = add_cvar!(model, c, -Inf, Inf)

add_bvar!(model::Model, c::Real) = add_var!(model, GRB_BINARY, c, 0., 1.)

add_ivar!(model::Model, c::Real, lb::Real, ub::Real) = add_var!(model, GRB_INTEGER, c, lb, ub)
add_ivar!(model::Model, c::Real) = add_ivar!(model, c, -Inf, Inf)

# add_vars!

function add_vars!(model::Model, vtypes::CVec, c::FVec, lb::FVec, ub::FVec)

    # check dimensions
    n = length(vtypes)
    _chklen(c, n)
    _chklen(lb, n)
    _chklen(ub, n)

    # main call
    ret = @grb_ccall(addvars, Cint, (
        Ptr{Void},  # model
        Cint,       # numvars
        Cint,       # numnz
        Ptr{Cint},  # vbeg
        Ptr{Cint},  # vind
        Ptr{Float64}, # vval
        Ptr{Float64}, # obj
        Ptr{Float64}, # lb
        Ptr{Float64}, # ub
        Ptr{Cchar},   # vtypes
        Ptr{Ptr{UInt8}}, # varnames
        ),
        model, n, 0, C_NULL, C_NULL, C_NULL, c, lb, ub, vtypes, C_NULL)

    if ret != 0
        throw(GurobiError(model.env, ret))
    end
    nothing
end

function add_vars!(model::Model, vtypes::GCharOrVec, c::Vector, lb::Bounds, ub::Bounds)
    n = length(c)
    add_vars!(model, cvecx(vtypes, n), fvec(c), fvecx(lb, n), fvecx(ub, n))
end

add_cvars!(model::Model, c::Vector, lb::Bounds, ub::Bounds) = add_vars!(model, GRB_CONTINUOUS, c, lb, ub)
add_cvars!(model::Model, c::Vector) = add_cvars!(model, c, -Inf, Inf)

add_bvars!(model::Model, c::Vector) = add_vars!(model, GRB_BINARY, c, 0, 1)

add_ivars!(model::Model, c::Vector, lb::Bounds, ub::Bounds) = add_vars!(model, GRB_INTEGER, c, lb, ub)
add_ivars!(model::Model, c::Vector) = add_ivars!(model, GRB_INTEGER, c, -Inf, Inf)


text:
# Exports

export Shape, PostOpObj, create, combine

abstract type Shape end

mutable struct PostOpObj
    objs::Vector{Any}
    operations::Array{Function}
    function PostOpObj(obj, func)
        if ~isa(obj, Vector)
            obj = [obj]
        end
        if ~isa(func, Vector)
                func = [func]
        end
        new(obj, func)
    end
end

function Base.copy(x::T) where T<:Union{Shape, PostOpObj}
    T(copy.([getfield(x, fn) for fn in fieldnames(T)])...)
end

"""
    create(shape::T; resolution=nothing, rand_=0.0, type::Int64=1) where T <: Shape

Abstact function for creating **Shape** objects.

## Returns
    - X : Initial reference position
    - V : Initial velocity
    - Y : Initial position
    - volume : Volume per particle point
    - type: Type of particle point

"""
function create(shape::T; resolution=nothing, rand_=0.0, type::Int64=1) where T <: Shape
    error("Not implemented for type **$(typeof(shape))** yet.")
end

function mycat(a, b)
    if length(size(a))==1
        return vcat(a, b)
    else
        return hcat(a, b)
    end
end

function create(pobj::PostOpObj, args...; kwargs...)
    temp = nothing
    out = nothing
    for obj in pobj.objs
        temp = create(obj, args...; kwargs...)
        if isa(out, Nothing)
            out = temp
        else
            out = mycat.(out, temp)
        end
    end
    for func in pobj.operations
        out = func(out)
    end
    return out
end

function combine(obj1::T1, obj2::T2) where {T1<:Union{Shape, PostOpObj}, T2<:Union{Shape, PostOpObj}}
    PostOpObj([obj1, obj2], [])
end

include("./cuboid.jl")
include("./sphere.jl")
include("./disk.jl")
include("./cylinder.jl")
include("./cone.jl")

max_stars_repo_name: ravinderbhattoo/PDMesh

max_issues_repo_name: ravinderbhattoo/PDMesh

max_forks_repo_name: ravinderbhattoo/PDMesh


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:4] <= 1)

                     @variable(m, -6.019960436071491 <= q <= 11.696817216382863)

                     add_NL_constraint(m, :(softplus(0.9042935565441375 + 0.31236426273908524*softplus(0.5868252904954097 + 0.04712499510197299*$(x[1]) + 0.4005456947892885*$(x[2]) + 0.34544434364772814*$(x[3]) + 0.48957180414623247*$(x[4])) + -0.4128166127064836*softplus(0.5797949667252156 + -0.7245574143823843*$(x[1]) + -0.08207866517990237*$(x[2]) + 0.4631961331947827*$(x[3]) + 0.9179781125579716*$(x[4])) + -0.2476169775927386*softplus(0.8773151656273224 + 0.003673243009740812*$(x[1]) + 0.698013040735662*$(x[2]) + 0.6787695040522581*$(x[3]) + -0.9128403322093082*$(x[4]))) + softplus(0.8600467319105785 + -0.4006275013725942*softplus(0.5868252904954097 + 0.04712499510197299*$(x[1]) + 0.4005456947892885*$(x[2]) + 0.34544434364772814*$(x[3]) + 0.48957180414623247*$(x[4])) + 0.22106887260345687*softplus(0.5797949667252156 + -0.7245574143823843*$(x[1]) + -0.08207866517990237*$(x[2]) + 0.4631961331947827*$(x[3]) + 0.9179781125579716*$(x[4])) + 0.7210213948433983*softplus(0.8773151656273224 + 0.003673243009740812*$(x[1]) + 0.698013040735662*$(x[2]) + 0.6787695040522581*$(x[3]) + -0.9128403322093082*$(x[4]))) + softplus(0.982756228863964 + 0.34792310275142935*softplus(0.5868252904954097 + 0.04712499510197299*$(x[1]) + 0.4005456947892885*$(x[2]) + 0.34544434364772814*$(x[3]) + 0.48957180414623247*$(x[4])) + 0.7924515917531365*softplus(0.5797949667252156 + -0.7245574143823843*$(x[1]) + -0.08207866517990237*$(x[2]) + 0.4631961331947827*$(x[3]) + 0.9179781125579716*$(x[4])) + -0.9399735455074207*softplus(0.8773151656273224 + 0.003673243009740812*$(x[1]) + 0.698013040735662*$(x[2]) + 0.6787695040522581*$(x[3]) + -0.9128403322093082*$(x[4]))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
############################################################################################
#                              OBJECTIVES                                                  #
############################################################################################

abstract type AbstractObjective end
Base.length(obj::AbstractObjective) = length(obj.cost)

"""```
cost(obj::Objective, Z::Traj)::Float64
cost(obj::Objective, dyn_con::DynamicsConstraint{Q}, Z::Traj)
```
Evaluate the cost for a trajectory.
Calculate the cost gradient for an entire trajectory. If a dynamics constraint is given,
    use the appropriate integration rule, if defined.
"""
function cost(obj::AbstractObjective, Z)
    cost!(obj, Z)
    J = get_J(obj)
    return sum(J)
end

"""$(TYPEDEF)
Objective: stores stage cost(s) and terminal cost functions

Constructors:
```julia
Objective(cost, N)
Objective(cost, cost_term, N)
Objective(costs::Vector{<:CostFunction}, cost_term)
Objective(costs::Vector{<:CostFunction})
```
"""
struct Objective{C} <: AbstractObjective
    cost::Vector{C}
    J::Vector{Float64}
end

# Constructors
function Objective(cost::CostFunction,N::Int)
    Objective([cost for k = 1:N], zeros(N))
end

function Objective(cost::CostFunction,cost_terminal::CostFunction,N::Int)
    Objective([k < N ? cost : cost_terminal for k = 1:N], zeros(N))
end

function Objective(cost::Vector{<:CostFunction},cost_terminal::CostFunction)
    N = length(cost) + 1
    Objective([cost...,cost_terminal], zeros(N))
end

function Objective(cost::Vector{<:CostFunction})
    N = length(cost)
    Objective(cost, zeros(N))
end

# Methods
"Get the vector of costs at each knot point. `sum(get_J(obj))` is equal to the cost"
get_J(obj::Objective) = obj.J

Base.copy(obj::Objective) = Objective(copy(obj.cost), copy(obj.J))

Base.getindex(obj::Objective,i::Int) = obj.cost[i]

Base.iterate(obj::Objective, start=1) = Base.iterate(obj.cost, start)

Base.show(io::IO, obj::Objective{C}) where C = print(io,"Objective")

# Convenience constructors
@doc raw"""```julia
LQRObjective(Q, R, Qf, xf, N)
```
Create an objective of the form
`` (x_N - x_f)^T Q_f (x_N - x_f) + \sum_{k=0}^{N-1} (x_k-x_f)^T Q (x_k-x_f) + u_k^T R u_k``
"""
function LQRObjective(Q::AbstractArray, R::AbstractArray, Qf::AbstractArray,
        xf::AbstractVector, N::Int; checks=true, uf=zeros(size(R,1)))
    H = zeros(size(R,1),size(Q,1))
    q = -Q*xf
    r = -R*uf
    c = 0.5*xf'*Q*xf + 0.5*uf'R*uf
    qf = -Qf*xf
    cf = 0.5*xf'*Qf*xf

    ℓ = QuadraticCost(Q, R, H, q, r, c, checks=checks)
    ℓN = QuadraticCost(Qf, qf, cf, check=checks)

    Objective(ℓ, ℓN, N)
end

function LQRObjective(
        Q::Union{Diagonal{T,<:SVector{n}},SMatrix{n,n}},
        R::Union{Diagonal{T,<:SVector{m}},SMatrix{m,m}},
        Qf::AbstractArray, xf::AbstractVector, N::Int;
        uf=(@SVector zeros(m)),
        checks=true) where {T,n,m}
    H = @SMatrix zeros(m,n)
    q = -Q*xf
    r = -R*uf
    c = 0.5*xf'*Q*xf + 0.5*uf'R*uf
    qf = -Qf*xf
    cf = 0.5*xf'*Qf*xf

    ℓ = QuadraticCost(Q, R, H, q, r, c, checks=checks)

    ℓN = QuadraticCost(Qf, R, H, qf, r, cf, checks=checks)

    Objective(ℓ, ℓN, N)
end


text:
module GenericInstruments

# lib locations
include("lib_locations.jl")
# common instrument containers
include("INSTR_TYPES.jl")
# visa wrapper,constants and utils
include("visa/VISA.jl")
# instruments
include("psu/PSU.jl")
include("awg/AWG.jl")
include("scope/SCOPE.jl")
include("dmm/DMM.jl")

using .PSU
using .AWG
using .SCOPE
using .DMM

end #endmodule


text:
# Data Section
# ============

# Each section is compressed separately.
struct Section
    chromid::UInt32
    chromstart::UInt32
    chromend::UInt32
    offset::UInt64
    size::UInt64
end