text:
using SpecialFunctions
using Base.MathConstants
using DataFrames, CSV
using Printf

"""
    gamma_derivatives()

Derive a table of derivatives of the Gamma function at 1, i.e., Gamma^{(m)}(a)

This code is used to create the results, which are included in the main code via a lookup table.

"""
function gamma_derivatives(k::Integer, x::Float64)
    #    general form from: http://erikerlandson.github.io/blog/2016/06/15/computing-derivatives-of-the-gamma-function/
    return SpecialFunctions.gamma(x) * D(k,0,x)
end

function D(k::Integer, n::Integer, x::Float64)
    # general form from: http://erikerlandson.github.io/blog/2016/06/15/computing-derivatives-of-the-gamma-function/
    if n<0
        throw(DomainError(n))
    elseif k<0
        throw(DomainError(k))
    end
    if k==0
        return 1
    elseif k==1
        return SpecialFunctions.polygamma(n,x)
    else
        total = D(k-1,n+1,x)
        for j=0:n
            total += binomial(n,j) * D(1,j,x) * D(k-1,n-j,x)
        end
        return total
    end
end

function g1(t::Int) # Crandall,2012, p.17
    # t derivate of Gamma function at 1
    PG21 = -2.4041138063191902 # SpecialFunctions.polygamma(2,1)
    PG41 = -24.8862661234409   # SpecialFunctions.polygamma(4,1)
    if t==0
        return 1.0
    elseif t==1
        return -γ
    elseif t==2
        # return γ^2 - γ + pi^2/6, Crandall seems to be wrong on this
        return γ^2 + pi^2/6  # from Mathematica
    elseif t==3
        # return -2*γ^3 + 9*γ^2 - (π^2+6)*γ + 3*π^2/2 - 4*SpecialFunctions.zeta(3), Crandall seems to be wrong on this
        return -γ^3  - π^2*γ/2 + PG21 # from Mathematica
    elseif t==4
        return γ^4 + γ^2*π^2 + 3*π^4/20 - 4*γ*PG21 # from Mathematica
    elseif t==5
        return -γ^5 - (20/12)*γ^3*π^2 - (9/12)*γ*π^4 + (10*γ^2 + (20/12)*π^2)*PG21 + PG41 # from Mathematica
    else
        return NaN
    end
end

function g2(t::Int) # http://erikerlandson.github.io/blog/2016/06/15/computing-derivatives-of-the-gamma-function/
    # t derivate of Gamma function at 1
    return gamma_derivatives(t,1.0)
end

# #### output a table
# K = 14
# k = collect(0:K)
# G1 = zeros(Float64, K+1)
# G2 = zeros(Float64, K+1)
# for i=1:K+1
#     G1[i] = g1(i-1)
#     G2[i] = g2(i-1)
#     @printf("G1[i] = %.20f, G2[i] = %.20f\n", G1[i], G2[i])
# end

# df = DataFrame( k=k, G1=G1, G2=G2 )
# CSV.write("../data/gamma_derivatives.csv", df)


text:
"""Code generator which aims at automating idenfitied patterns.
"""

abstract type Declaration end

generate(cgs::Declaration) = error("No generator available for type $(typeof(cgs))")

struct Statement
    body::AbstractString
    assigned_id
end

mutable struct SDefinition <: Declaration
    name::AbstractString
    is_mutable::Bool
    fields::OrderedDict
    abstract_type
    inner_constructor
    function SDefinition(name, is_mutable, fields, abstract_type, inner_constructor)
        if !isnothing(inner_constructor)
            @assert inner_constructor.name == name
        end
        new(name, is_mutable, fields, abstract_type, inner_constructor)
    end
end

struct FDefinition <: Declaration
    name::AbstractString
    signature::Signature
    short::Bool
    body::AbstractArray{Statement}
    docstring::AbstractString
    return_type::AbstractString
end

struct CDefinition <: Declaration
    name::AbstractString
    value
end

@with_kw struct EDefinition <: Declaration
    name::AbstractString
    fields
    with_begin_block::Bool = length(fields) > 8
    type = nothing
    enum_macro::AbstractString = "@enum"
end

OrderedDict(defs::AbstractArray{T}) where {T <: Declaration} = OrderedDict{AbstractString,T}(map(x -> x.name => x, defs))

generate(cdef::CDefinition) = format_text("const $(cdef.name) = $(cdef.value)")
function generate(edef::EDefinition)
    bg = edef.with_begin_block
    format_text("$(edef.enum_macro) $(typed_field(edef.name, edef.type)) $(bg ? "begin\n" : "") $(join(edef.fields, bg ? "\n" : " ")) $(bg ? "\nend" : "")")
end

Statement(body::AbstractString) = Statement(strip(body), nothing)
SDefinition(name::AbstractString, is_mutable::Bool, fields::OrderedDict) = SDefinition(name, is_mutable, fields, nothing, nothing)
SDefinition(name::AbstractString, is_mutable::Bool; fields=()) = SDefinition(name, is_mutable, OrderedDict(fields))
Signature(sdef::SDefinition) = isnothing(sdef.inner_constructor) ? Signature(sdef.name, PositionalArgument.(keys(sdef.fields), values(sdef.fields)), KeywordArgument[]) : sdef.inner_constructor.signature

function extract_args(str)
    split_str = split(str, ";") # get kwargs first
    @assert length(split_str) <= 2 "More than one kwargs separator ';' in arguments ($str)"
    args, kwargs = first(split_str), length(split_str) == 1 ? nothing : split_str[2]
    splitstrip(args), splitstrip(kwargs)
end

function FDefinition(str::AbstractString)
    short = !startswith(str, "function")
    if short
        id, args, body = match(r"([\w\.]+)\((.*)\)\s+=(?!=)(.*)$", str).captures
        body = [Statement(body)]
    else
        id, args = match(r"function ([\w\.]+)\((.*)\)", str).captures
        body = Statement.(split(str, "\n")[2:end - 1]) # remove function... and end lines
    end
    if occursin(";", args)
        args_str, kwargs_str = split(args, ";")
        kwargs = KeywordArgument.(splitargs(kwargs_str))
    else
        args_str, kwargs_str = args, []
        kwargs = []
    end
    args = PositionalArgument.(splitargs(args_str))
    sig = Signature(id, args, kwargs)
    FDefinition(id, sig, short, body)
end

FDefinition(name::AbstractString, signature::Signature, short::Bool, body::AbstractArray{Statement}) = FDefinition(name, signature, short, body, "", "")
FDefinition(name::AbstractString, signature::Signature, short::Bool, body::AbstractArray{Statement}, docstring::AbstractString) = FDefinition(name, signature, short, body, docstring, "")

function CDefinition(str::AbstractString)
    split_str = split(str, " ")
    id = split_str[2]
    value_multiline = splitjoin(str, [1], delim="=")
    value = replace(join(splitstrip(value_multiline; delim="\n"), "\n"), "\n" => " ") # remove any whitespace for each line
    CDefinition(id, value)
end

function EDefinition(str::AbstractString)
    split_str = split(str, " ")
    enum_macro, id = split_str[1:2]
    id, type = decompose_field_decl(id)
    split_str_n = splitstrip(str, delim="\n")
    with_begin_block = last(split(split_str_n[1], " ")) == "begin"
    values = with_begin_block ? strip.(split(join(split_str_n[2:end - 1], "\n"), "\n")) : strip.(splitjoin(str, [1, 2], delim=" "))
    EDefinition(id, values, with_begin_block, type, enum_macro)
end

function SDefinition(str::AbstractString)
    split_str = strip.(split(str, "\n"))
    def = split_str[1]
    def_parts = split(def, " ")
    is_mutable = def_parts[1] == "mutable"
    name = is_mutable ? def_parts[3] : def_parts[2]
    field_decls = OrderedDict(decompose_field_decl(field) for field ∈ split_str[2:end - 1]) # skip struct... and end lines
    SDefinition(name, is_mutable, field_decls)
end

function generate(statements::AbstractArray{Statement})
    isempty(statements) && return ""
    format_text(join(getproperty.(statements, :body), "\n"))
end

function generate(s::SDefinition)
    def = (s.is_mutable ? "mutable " : "") * "struct $(s.name)" * (isnothing(s.abstract_type) ? "" : "<: $(s.abstract_type)")
    fields = join(typed_field.(keys(s.fields), values(s.fields)), "\n")
    inner_constructor = isnothing(s.inner_constructor) ? "" : "\n" * generate(s.inner_constructor) * "\n"
    format_text("$def $fields $inner_constructor end")
end

function generate(f::FDefinition)
    body = generate(f.body)
    docstring = isempty(f.docstring) ? "" : """\"\"\"
    $(f.docstring)
    \"\"\"\n"""
    if f.short
        str = docstring * generate(f.signature) * " = " * body
    else
        str = docstring * "function $(generate(f.signature)) $(body) end"
    end
    format_text(str)
end

ext: jl


text:
using Test, Pkg
using MPI

if get(ENV,"JULIA_MPI_TEST_ARRAYTYPE","") == "CuArray"
    using CuArrays
    ArrayType = CuArray
else
    ArrayType = Array
end

MPI.Init()

comm_size = MPI.Comm_size(MPI.COMM_WORLD)

if ArrayType != Array || Sys.iswindows() && Sys.WORD_SIZE == 32
    operators = [MPI.SUM, +]
else
    operators = [MPI.SUM, +, (x,y) -> 2x+y-x]
end

for T = [Int]
    for dims = [1, 2, 3]
        send_arr = ArrayType(zeros(T, Tuple(3 for i in 1:dims)))
        send_arr[:] .= 1:length(send_arr)

        for op in operators

            # Non allocating version
            recv_arr = ArrayType{T}(undef, size(send_arr))
            MPI.Allreduce!(send_arr, recv_arr, op, MPI.COMM_WORLD)
            @test recv_arr == comm_size .* send_arr

            # Assertions when output buffer too small
            recv_arr = ArrayType{T}(undef, size(send_arr).-1)
            @test_throws AssertionError MPI.Allreduce!(send_arr, recv_arr,
                                                       length(send_arr),
                                                       op, MPI.COMM_WORLD)
            # IN_PLACE
            recv_arr = copy(send_arr)
            MPI.Allreduce!(recv_arr, op, MPI.COMM_WORLD)
            @test recv_arr == comm_size .* send_arr

            # Allocating version
            val = MPI.Allreduce(2, op, MPI.COMM_WORLD)
            @test val == comm_size*2

            vals = MPI.Allreduce(send_arr, op, MPI.COMM_WORLD)
            @test vals isa ArrayType{T}
            @test vals == comm_size .* send_arr
        end
    end
end

MPI.Barrier( MPI.COMM_WORLD )

GC.gc()
MPI.Finalize()
@test MPI.Finalized()


text:
"""
	dt = DotTheta( (x,y) -> dot(x,y) / length(x) )

This parametric type allows to define a new dot product from the one saved in `dt::dot`. More precisely:

	dt(u1, u2, p1::T, p2::T, theta::T) where {T <: Real}

computes, the weigthed dot product ``\\langle (u_1,p_1), (u_2,p_2)\\rangle_\\theta = \\theta \\Re \\langle u_1,u_2\\rangle  +(1-\\theta)p_1p_2`` where ``u_i\\in\\mathbb R^N``. The ``\\Re`` factor is put to ensure a real valued result despite possible complex valued arguments.

	normtheta(u, p::T, theta::T)

Compute, the norm associated to weighted dot product ``\\langle (u_1,p_1), (u_2,p_2)\\rangle_\\theta``.

!!! info "Info"
    This is used in the pseudo-arclength constraint with the dot product ``\\frac{1}{N} \\langle u_1,u_2\\rangle,\\quad u_i\\in\\mathbb R^N``
"""
struct DotTheta{Tdot}
	dot::Tdot		# defaults to (x,y) -> dot(x,y) / length(x)
end

DotTheta() = DotTheta( (x, y) -> dot(x, y) / length(x) )

# Implementation of the dot product associated to DotTheta
# we restrict the type of the parameters because for complex problems, we still want the parameter to be real
(dt::DotTheta)(u1, u2, p1::T, p2::T, θ::T) where {T <: Real} = real(dt.dot(u1, u2) * θ + p1 * p2 * (one(T) - θ))

# Implementation of the norm associated to DotTheta
# we restrict the type of the parameters because for complex problems, we still want the parameter to be real
(dt::DotTheta)(u, p::T, θ::T) where T = sqrt(dt(u, u, p, p, θ))

(dt::DotTheta)(a::BorderedArray{vec, T}, b::BorderedArray{vec, T}, θ::T) where {vec, T} = dt(a.u, b.u, a.p, b.p, θ)
(dt::DotTheta)(a::BorderedArray{vec, T}, θ::T) where {vec, T} = dt(a.u, a.p, θ)
####################################################################################################
# equation of the arc length constraint
arcLengthEq(dt::DotTheta, u, p, du, dp, θ, ds) = dt(u, du, p, dp, θ) - ds
####################################################################################################
abstract type AbstractTangentPredictor end
abstract type AbstractSecantPredictor <: AbstractTangentPredictor end

# wrapper to use iterators and state
getPredictor!(state::AbstractContinuationState, iter::AbstractContinuationIterable, nrm = false) = getPredictor!(state.z_pred, state.z_old, state.tau, state.ds, iter.tangentAlgo, nrm)
getTangent!(state::AbstractContinuationState, it::AbstractContinuationIterable, verbosity) = getTangent!(state.tau, state.z_new, state.z_old, it, state.ds, θ, it.tangenAlgo::NaturalPred, it.verbosity)

# reset the predictor
emptypredictor!(::Union{Nothing, AbstractTangentPredictor}) = nothing

# clamp p value
clampPred(p::Number, it::AbstractContinuationIterable) = clamp(p, it.contParams.pMin, it.contParams.pMax)

# this function only mutates z_pred
# the nrm argument allows to just increment z_pred.p by ds
function getPredictor!(z_pred::M, z_old::M, tau::M, ds, algo::Talgo, nrm = false) where {T, vectype, M <: BorderedArray{vectype, T}, Talgo <: AbstractTangentPredictor}
	# we perform z_pred = z_old + ds * tau
	copyto!(z_pred, z_old) # z_pred .= z_old
	nrm ? axpy!(ds / tau.p, tau, z_pred) : axpy!(ds, tau, z_pred)
end

# generic corrector based on Bordered formulation
function corrector(it, z_old::M, τ::M, z_pred::M, ds, θ,
			algo::Talgo, linearalgo = MatrixFreeBLS();
			normC = norm, callback = cbDefault, kwargs...) where
			{T, vectype, M <: BorderedArray{vectype, T}, Talgo <: AbstractTangentPredictor}
	if z_pred.p <= it.contParams.pMin || z_pred.p >= it.contParams.pMax
		z_pred.p = clampPred(z_pred.p, it)
		return corrector(it, z_old, τ, z_pred, ds, θ, NaturalPred(), linearalgo;
						normC = normC, callback = callback, kwargs...)
	end
	return newtonPALC(it, z_old, τ, z_pred, ds, θ; linearbdalgo = linearalgo, normN = normC, callback = callback, kwargs...)
end
####################################################################################################
"""
	Natural predictor / corrector
"""
struct NaturalPred <: AbstractTangentPredictor end

function getPredictor!(z_pred::M, z_old::M, τ::M, ds, algo::NaturalPred, nrm = false) where {T, vectype, M <: BorderedArray{vectype, T}}
	# we do z_pred .= z_old
	copyto!(z_pred, z_old) # z_pred .= z_old
	z_pred.p += ds
end

# corrector based on natural formulation
function corrector(it, z_old::M, τ::M, z_pred::M, ds, θ,
			algo::NaturalPred, linearalgo = MatrixFreeBLS();
			normC = norm, callback = cbDefault, kwargs...) where
			{T, vectype, M <: BorderedArray{vectype, T}}
	res = newton(it.F, it.J, z_pred.u, setParam(it, clampPred(z_pred.p, it)), it.contParams.newtonOptions; normN = normC, callback = callback, kwargs...)
	return BorderedArray(res[1], z_pred.p), res[2:end]...
end

function getTangent!(τ::M, z_new::M, z_old::M, it::AbstractContinuationIterable, ds, θ, algo::NaturalPred, verbosity) where {T, vectype, M <: BorderedArray{vectype, T}}
	(verbosity > 0) && println("--> predictor = ", algo)
	# we do nothing here, the predictor will just copy z_old into z_pred
end
####################################################################################################
"""
	Secant tangent predictor
"""
struct SecantPred <: AbstractSecantPredictor end

# tangent computation using Secant predictor
# tau is the tangent prediction
function getTangent!(τ::M, z_new::M, z_old::M, it::AbstractContinuationIterable, ds, θ, algo::SecantPred, verbosity) where {T, vectype, M <: BorderedArray{vectype, T}}
	(verbosity > 0) && println("--> predictor = ", algo)
	# secant predictor: tau = z_new - z_old; tau *= sign(ds) / normtheta(tau)
	copyto!(τ, z_new)
	minus!(τ, z_old)
	α = algo isa SecantPred ? sign(ds) / it.dottheta(τ, θ) : sign(ds) / abs(τ.p)
	rmul!(τ, α)
end
####################################################################################################
"""
	Bordered Tangent predictor
"""
struct BorderedPred <: AbstractTangentPredictor end

# tangent computation using Bordered system
# tau is the tangent prediction
function getTangent!(τ::M, z_new::M, z_old::M, it::AbstractContinuationIterable, ds, θ, algo::BorderedPred, verbosity) where {T, vectype, M <: BorderedArray{vectype, T}}
	(verbosity > 0) && println("--> predictor = Bordered")
	# tangent predictor
	ϵ = it.contParams.finDiffEps
	# dFdl = (F(z_new.u, z_new.p + ϵ) - F(z_new.u, z_new.p)) / ϵ
	dFdl = it.F(z_new.u, setParam(it, z_new.p + ϵ))
	minus!(dFdl, it.F(z_new.u, setParam(it, z_new.p)))
	rmul!(dFdl, 1/ϵ)

	# tau = getTangent(J(z_new.u, z_new.p), dFdl, tau_old, theta, contparams.newtonOptions.linsolve)
	tau_normed = copy(τ)#copyto!(similar(tau), tau) #copy(tau_old)
	rmul!(tau_normed, θ / length(τ.u), 1 - θ)
	# extract tangent as solution of bordered linear system, using zero(z_new.u)
	tauu, taup, flag, itl = it.linearAlgo( it.J(z_new.u, setParam(it, z_new.p)), dFdl,
			tau_normed, 0*z_new.u, T(1), θ)

	# the new tangent vector must preserve the direction along the curve
	α = T(1) / it.dottheta(tauu, τ.u, taup, τ.p, θ)

	# tau_new = α * tau
	copyto!(τ.u, tauu)
	τ.p = taup
	rmul!(τ, α)
end
####################################################################################################
"""
	Multiple Tangent predictor
$(TYPEDFIELDS)

# Constructor(s)

	MultiplePred(algo, x0, α, n)

	MultiplePred(x0, α, n)

- `α` damping in Newton iterations
- `n` number of predictors
- `x0` example of vector solution to be stored
"""
@with_kw mutable struct MultiplePred{T <: Real, Tvec, Talgo} <: AbstractTangentPredictor
	"Tangent algorithm used"
	tangentalgo::Talgo

	"Save the current tangent"
	τ::Tvec

	"Damping factor"
	α::T

	"Number of predictors"
	nb::Int64

	"Index of the largest converged predictor"
	currentind::Int64 = 0

	"Index for lookup in residual history"
	pmimax::Int64 = 1

	"Maximum index for lookup in residual history"
	imax::Int64 = 4

	"Factor to increase ds upon successful step"
	dsfact::T = 1.5
end
MultiplePred(algo::AbstractTangentPredictor, x0, α::T, nb) where T = MultiplePred(tangentalgo= algo, τ = BorderedArray(x0, T(0)), α = α, nb = nb)
MultiplePred(x0, α, nb) = MultiplePred(SecantPred(), x0, α, nb)
emptypredictor!(mpd::MultiplePred) = (mpd.currentind = 1; mpd.pmimax = 1)

# callback for newton
function (mpred::MultiplePred)(x, f, J, res, iteration, itlinear, options; kwargs...)
	resHist = get(kwargs, :resHist, nothing)
	if mpred.currentind > 1
		return iteration - mpred.pmimax > 0 ?  resHist[end] <= mpred.α * resHist[end-mpred.pmimax] : true
	end
	return true
end

function getTangent!(τ::M, z_new::M, z_old::M, it::AbstractContinuationIterable, ds, θ, algo::MultiplePred{T, M, Talgo}, verbosity) where {T, vectype, M <: BorderedArray{vectype, T}, Talgo}
	# compute tangent and store it
	(verbosity > 0) && print("--> predictor = MultiplePred\n--")
	getTangent!(τ, z_new, z_old, it, ds, θ, algo.tangentalgo, verbosity)
	copyto!(algo.τ, τ)
end

function getPredictor!(z_pred::M, z_old::M, τ::M, ds, algo::MultiplePred, nrm = false) where {T, vectype, M <: BorderedArray{vectype, T}}
	# we do nothing!
	# emptypredictor!(algo)
	return nothing
end

function corrector(it, z_old::M, tau::M, z_pred::M, ds, θ,
		algo::MultiplePred, linearalgo = MatrixFreeBLS(); normC = norm,
		callback = cbDefault, kwargs...) where {T, vectype, M <: BorderedArray{vectype, T}}
	# we combine the callbacks for the newton iterations
	cb = (x, f, J, res, iteration, itlinear, options; k...) -> callback(x, f, J, res, iteration, itlinear, options; k...) & algo(x, f, J, res, iteration, itlinear, options; k...)
	# note that z_pred already contains ds * τ, hence ii=0 corresponds to this case
	for ii in algo.nb:-1:1
		# printstyled(color=:magenta, "--> ii = $ii\n")
		# record the current index
		algo.currentind = ii
		zpred = _copy(z_pred)
		axpy!(ii * ds, algo.τ, zpred)
		# we restore the original callback if it reaches the usual case ii == 0
		zold, res, flag, itnewton, itlinear = corrector(it, z_old, tau, zpred, ds, θ,
				algo.tangentalgo, linearalgo; normC = normC, callback = cb, kwargs...)
		if flag || ii == 1 # for i==1, we return the result anyway
			return zold, res, flag, itnewton, itlinear
		end
	end
	return zold, res, flag, itnewton, itlinear
end

function stepSizeControl(ds, θ, contparams::ContinuationPar, converged::Bool, it_newton_number::Int, tau::M, mpd::MultiplePred, verbosity) where {T, vectype, M<:BorderedArray{vectype, T}}
	if converged == false
		dsnew = ds
		if abs(ds) < (1 + mpd.nb) * contparams.dsmin
			if mpd.pmimax < mpd.imax
				@error "--> Increase pmimax"
				mpd.pmimax += 1
			else
				(verbosity > 0) && printstyled("*"^80*"\nFailure to converge with given tolerances\n"*"*"^80, color=:red)
				# we stop the continuation
				return ds, θ, true
			end
		else
			@error "--> Decrease ds"
			dsnew = ds / (1 + mpd.nb)
			(verbosity > 0) && printstyled("Halving continuation step, ds = $(dsnew)\n", color=:red)
		end
	else # the newton correction has converged
		dsnew = ds
		if mpd.currentind == mpd.nb && abs(ds) * mpd.dsfact <= contparams.dsmax
			(verbosity > 0) && @show dsnew
			# println("--> Increase ds")
			dsnew = ds *  mpd.dsfact
		end
	end

	# control step to stay between bounds
	dsnew = clampDs(dsnew, contparams)

	return dsnew, θ, false
end
####################################################################################################
"""
	Polynomial Tangent predictor

$(TYPEDFIELDS)

# Constructor(s)

	PolynomialPred(algo, n, k, v0)

	PolynomialPred(n, k, v0)

- `n` order of the polynomial
- `k` length of the last solutions vector used for the polynomial fit
- `v0` example of solution to be stored. It is only used to get the `eltype` of the tangent!!
"""
mutable struct PolynomialPred{T <: Real, Tvec, Talgo} <: AbstractTangentPredictor
	"Order of the polynomial"
	n::Int64

	"Length of the last solutions vector used for the polynomial fit"
	k::Int64

	"Matrix for the interpolation"
	A::Matrix{T}

	"Algo for tangent when polynomial predictor is not possible"
	tangentalgo::Talgo

	"Vector of solutions"
	solutions::CircularBuffer{Tvec}

	"Vector of parameters"
	parameters::CircularBuffer{T}

	"Vector of arclengths"
	arclengths::CircularBuffer{T}

	"Coefficients for the polynomials for the solution"
	coeffsSol::Vector{Tvec}

	"Coefficients for the polynomials for the parameter"
	coeffsPar::Vector{T}

	"Update the predictor ?"
	update::Bool
end

PolynomialPred(algo, n, k, v0) = (@assert n<k; ;PolynomialPred(n,k,zeros(eltype(v0),k,n+1),algo,
	CircularBuffer{typeof(v0)}(k),CircularBuffer{eltype(v0)}(k),
	CircularBuffer{eltype(v0)}(k),
	Vector{typeof(v0)}(undef, n+1),
	Vector{eltype(v0)}(undef, n+1),true))
PolynomialPred(n, k, v0) = PolynomialPred(SecantPred(),n,k,v0)

isready(ppd::PolynomialPred) = length(ppd.solutions) >= ppd.k

function emptypredictor!(ppd::PolynomialPred)
	empty!(ppd.solutions); empty!(ppd.parameters); empty!(ppd.arclengths);
end

function getStats(polypred::PolynomialPred)
	Sbar = sum(polypred.arclengths) / polypred.n
	σ = sqrt(sum(x->(x-Sbar)^2, polypred.arclengths ) / (polypred.n))
	# return 0,1
	return Sbar, σ
end

function (polypred::PolynomialPred)(ds)
	sbar, σ = getStats(polypred)
	s = polypred.arclengths[end] + ds
	S = [((s-sbar)/σ)^(jj-1) for jj=1:polypred.n+1]
	p = sum(S .* polypred.coeffsPar)
	x = sum(S .* polypred.coeffsSol)
	return x, p
end

function updatePred!(polypred::PolynomialPred)
	Sbar, σ = getStats(polypred)
	Ss = (polypred.arclengths .- Sbar) ./ σ
	# construction of the Vandermond Matrix
	polypred.A[:, 1] .= 1
	for jj in 1:polypred.n
		polypred.A[:, jj+1] .= polypred.A[:, jj] .* Ss
	end
	# invert linear system for least square fitting
	B = (polypred.A' * polypred.A) \ polypred.A'
	mul!(polypred.coeffsSol, B, polypred.solutions)
	mul!(polypred.coeffsPar, B, polypred.parameters)
end

function getTangent!(tau::M, z_new::M, z_old::M, it::AbstractContinuationIterable, ds, θ, polypred::PolynomialPred, verbosity) where {T, vectype, M <: BorderedArray{vectype, T}, Talgo}
	# compute tangent and store it
	(verbosity > 0) && println("--> predictor = PolynomialPred")

	if polypred.update
		# update the list of solutions
		if length(polypred.arclengths)==0
			push!(polypred.arclengths, ds)
		else
			push!(polypred.arclengths, polypred.arclengths[end]+ds)
		end
		push!(polypred.solutions, z_new.u)
		push!(polypred.parameters, z_new.p)
	end

	if ~isready(polypred) || ~polypred.update
		return getTangent!(tau, z_new, z_old, it, ds, θ, polypred.tangentalgo, verbosity)
	else
		return polypred.update ? updatePred!(polypred) : true
	end
end

function getPredictor!(z_pred::M, z_old::M, tau::M, ds, polypred::PolynomialPred, nrm = false) where {T, vectype, M <: BorderedArray{vectype, T}}
	if ~isready(polypred)
		return getPredictor!(z_pred, z_old, tau, ds, polypred.tangentalgo, nrm)
	else
		x, p = polypred(ds)
		copyto!(z_pred.u, x)
		z_pred.p = p
		return true
	end
end
####################################################################################################
function arcLengthScaling(θ, contparams, tau::M, verbosity) where {M <: BorderedArray}
	# the arclength scaling algorithm is based on Salinger, Andrew G, Nawaf M Bou-Rabee,
	# Elizabeth A Burroughs, Roger P Pawlowski, Richard B Lehoucq, Louis Romero, and Edward D
	# Wilkes. “LOCA 1.0 Library of Continuation Algorithms: Theory and Implementation Manual,
	# ” March 1, 2002. https://doi.org/10.2172/800778.
	thetanew = θ
	g = abs(tau.p * θ)
	(verbosity > 0) && print("Theta changes from $(θ) to ")
	if (g > contparams.gMax)
		thetanew = contparams.gGoal / tau.p * sqrt( abs(1.0 - g^2) / abs(1.0 - tau.p^2) )
		if (thetanew < contparams.thetaMin)
		  thetanew = contparams.thetaMin;
		end
	end
	(verbosity > 0) && print("$(thetanew)\n")
	return thetanew
end
####################################################################################################
function clampDs(dsnew, contparams::ContinuationPar)
	return sign(dsnew) * clamp(abs(dsnew), contparams.dsmin, contparams.dsmax)
end

function stepSizeControl(ds, θ, contparams::ContinuationPar, converged::Bool, it_newton_number::Int, tau::M, algo::AbstractTangentPredictor, verbosity) where {T, vectype, M<:BorderedArray{vectype, T}}
	if converged == false
		if  abs(ds) <= contparams.dsmin
			(verbosity > 0) && printstyled("*"^80*"\nFailure to converge with given tolerances\n"*"*"^80, color=:red)
			# we stop the continuation
			return ds, θ, true
		end
		dsnew = sign(ds) * max(abs(ds) / 2, contparams.dsmin);
		(verbosity > 0) && printstyled("Halving continuation step, ds=$(dsnew)\n", color=:red)
	else
		# control to have the same number of Newton iterations
		Nmax = contparams.newtonOptions.maxIter
		factor = (Nmax - it_newton_number) / Nmax
		dsnew = ds * (1 + contparams.a * factor^2)
		# (verbosity > 0) && @show 1 + contparams.a * factor^2
	end

	# control step to stay between bounds
	dsnew = clampDs(dsnew, contparams)

	thetanew = contparams.doArcLengthScaling ? arcLengthScaling(θ, contparams, tau, verbosity) : θ
	# we do not stop the continuation
	return dsnew, thetanew, false
end
####################################################################################################
"""
This is the classical Newton-Krylov solver used to solve `F(x, p) = 0` together
with the scalar condition `n(x, p) ≡ θ ⋅ <x - x0, τx> + (1-θ) ⋅ (p - p0) * τp - n0 = 0`. This makes a problem of dimension N + 1.

Here, we specify `p` as a subfield of `par` with the `paramLens::Lens`

# Arguments
- `(x, par) -> F(x, par)` where `par` is a set of parameters like `(a=1.0, b=1)`
- `(x, par) -> Jh(x, par)` the jacobian Jh = ∂xF
"""
function newtonPALC(F, Jh, par, paramlens::Lens,
					z0::BorderedArray{vectype, T},
					τ0::BorderedArray{vectype, T},
					z_pred::BorderedArray{vectype, T},
					ds::T, θ::T,
					contparams::ContinuationPar{T, S, E},
					dottheta::DotTheta;
					linearbdalgo = BorderingBLS(DefaultLS()),
					normN = norm,
					callback = cbDefault, kwargs...) where {T, S, E, vectype}
	# Extract parameters
	@unpack tol, maxIter, verbose, α, αmin, linesearch = contparams.newtonOptions
	@unpack finDiffEps, pMin, pMax = contparams

	N = (x, p) -> arcLengthEq(dottheta, minus(x, z0.u), p - z0.p, τ0.u, τ0.p, θ, ds)
	normAC = (resf, resn) -> max(normN(resf), abs(resn))

	# Initialise iterations
	x = _copy(z_pred.u) 					# copy(z_pred.u)
	p = z_pred.p
	x_pred = _copy(x) 						# copy(x)

	# Initialise residuals
	res_f = F(x, set(par, paramlens, p));  res_n = N(x, p)

	dX = _copy(res_f) # copy(res_f)
	dp = T(0)
	up = T(0)
	# dFdp = (F(x, p + finDiffEps) - res_f) / finDiffEps
	dFdp = _copy(F(x, set(par, paramlens, p + finDiffEps)))
	minus!(dFdp, res_f)						# dFdp = dFdp - res_f
	rmul!(dFdp, T(1) / finDiffEps)

	res     = normAC(res_f, res_n)
	resHist = [res]
	it = 0
	itlineartot = 0

	# Displaying results
	verbose && displayIteration(it, res)
	line_step = true

	# invoke callback before algo really starts
	compute = callback(x, res_f, nothing, res, 0, 0, contparams; p = p, resHist = resHist, fromNewton = false, kwargs...)

	# Main loop
	while (res > tol) & (it < maxIter) & line_step & compute
		# dFdp = (F(x, p + epsi) - F(x, p)) / epsi)
		copyto!(dFdp, F(x, set(par, paramlens, p + finDiffEps)))
			minus!(dFdp, res_f); rmul!(dFdp, T(1) / finDiffEps)

		J = Jh(x, set(par, paramlens, p))
		u, up, flag, itlinear = linearbdalgo(J, dFdp, τ0, res_f, res_n, θ)
		itlineartot += sum(itlinear)

		if linesearch
			line_step = false
			while !line_step & (α > αmin)
				# x_pred = x - α * u
				copyto!(x_pred, x); axpy!(-α, u, x_pred)

				p_pred = p - α * up
				copyto!(res_f, F(x_pred, set(par, paramlens, p_pred)))

				res_n  = N(x_pred, p_pred)
				res = normAC(res_f, res_n)

				if res < resHist[end]
					if (res < resHist[end] / 2) & (α < 1)
						α *= 2
					end
					line_step = true
					copyto!(x, x_pred)
					p  = clamp(p_pred, pMin, pMax)
					# p = p_pred
				else
					α /= 2
				end
			end
			α = contparams.newtonOptions.α 	# we put back the initial value
		else
			minus!(x, u) 					# x .= x .- u
			p = clamp(p - up, pMin, pMax)
			# p = p - up

			copyto!(res_f, F(x, set(par, paramlens, p)))

			res_n  = N(x, p); res = normAC(res_f, res_n)
		end

		push!(resHist, res)
		it += 1

		verbose && displayIteration(it, res, itlinear)

		# shall we break the loop?
		compute = callback(x, res_f, J, res, it, itlinear, contparams; p = p, resHist = resHist, fromNewton = false, kwargs...)
	end
	verbose && displayIteration(it, res, 0, true) # display last line of the table
	flag = (resHist[end] < tol) & callback(x, res_f, nothing, res, it, -1, contparams; p = p, resHist = resHist, fromNewton = false, kwargs...)
	return BorderedArray(x, p), resHist, flag, it, itlineartot
end

# conveniency for use in continuation
newtonPALC(it::AbstractContinuationIterable, z0::M, τ0::M, z_pred::M, ds::T, θ::T; kwargs...) where {T, vectype, M <: BorderedArray{vectype, T}} = newtonPALC(it.F, it.J, it.par, it.lens, z0, τ0, z_pred, ds, θ, it.contParams, it.dottheta; kwargs...)
####################################################################################################


text:
@testset "path_info" begin
    @polyvar x y
    f = [x^2 + y^2 + 1, x + y - 3]
    tracker, starts = coretracker_startsolutions(f)
    ctinfo = path_info(tracker, first(starts))
    @test !isempty(sprint(show, ctinfo))

    tracker = pathtracker(f)
    ptinfo = path_info(tracker, first(starts))
    @test !isempty(sprint(show, ptinfo))
end


text:
## Aircraft analysis benchmarking
using BenchmarkTools
using StaticArrays

# All subsequent analyses use no symmetry tricks for performance as AeroMDAO hasn't implemented it and apples must be compared to apples.

## AeroMDAO tests: https://github.com/GodotMisogi/AeroMDAO.jl
#=======================================================#
using AeroMDAO

println("AeroMDAO Aircraft Functional -")
@benchmark begin
    # Wing
    wing = Wing(foils     = Foil.(fill(naca4((0,0,1,2)), 2)),
                chords    = [1.0, 0.6],
                twists    = [2.0, 2.0],
                spans     = [5.0],
                dihedrals = [11.31],
                sweep_LEs = [2.29]);

    # Horizontal tail
    htail = Wing(foils     = Foil.(fill(naca4((0,0,1,2)), 2)),
                 chords    = [0.7, 0.42],
                 twists    = [0.0, 0.0],
                 spans     = [1.25],
                 dihedrals = [0.],
                 sweep_LEs = [6.39])

    # Vertical tail
    vtail = HalfWing(foils     = Foil.(fill(naca4((0,0,0,9)), 2)),
                     chords    = [0.7, 0.42],
                     twists    = [0.0, 0.0],
                     spans     = [1.0],
                     dihedrals = [0.],
                     sweep_LEs = [7.97])

    wing_panels  = panel_wing(wing, 20, 20, spacing = "cosine")

    htail_panels = panel_wing(htail, 12, 12;
                              position = [4., 0, 0],
                              angle    = deg2rad(0.),
                              axis 	   = [0., 1., 0.],
                              spacing  = "cosine"
                             )

    vtail_panels = panel_wing(vtail, 12, 10;
                              position = [4., 0, 0],
                              angle    = π/2,
                              axis 	   = [1., 0., 0.],
                              spacing  = "cosine"
                             )

    aircraft = Dict("Wing" 			  	=> wing_panels,
                    "Horizontal Tail" 	=> htail_panels,
                    "Vertical Tail"   	=> vtail_panels)

    # display(size.([ wing_panels[1], htail_panels[1], vtail_panels[1] ])) # Checking sizes

    ρ       = 1.225
    x_ref   = [0.5, 0., 0.]
    V, α, β = 1.0, 5.0, 0.0
    Ω       = [0.0, 0.0, 0.0]
    fs      = AeroMDAO.Freestream(V, α, β, Ω)
    S, b, c = 9.0, 10.0, 0.9

    data = solve_case(aircraft, fs;
                      rho_ref   = ρ,
                      r_ref     = x_ref,
                      area_ref  = S,
                      span_ref  = b,
                      chord_ref = c)

    nf, ff = data["Aircraft"][1:2]

    nf[1:3], nf[4:6], ff[1]
end

## BYU FLOW Lab tests: https://github.com/byuflowlab/VortexLattice.jl
#=======================================================#

# Example taken from: https://flow.byu.edu/VortexLattice.jl/dev/examples/

# Pulled these out to avoid compilation for proper benchmarking
fc   = fill((xc) -> 0, 2) # camberline function for each section
fc_h = fill((xc) -> 0, 2) # camberline function for each section
fc_v = fill((xc) -> 0, 2) # camberline function for each section

##
using VortexLattice

println("BYU FLOW Lab VortexLattice.jl - ")
@benchmark begin
    # wing
    xle = [0.0, 0.2]
    yle = [0.0, 5.0]
    zle = [0.0, 1.0]
    chord = [1.0, 0.6]
    theta = [2.0*pi/180, 2.0*pi/180]
    phi = [0.0, 0.0]
    ns = 20
    nc = 20
    spacing_s = Cosine()
    spacing_c = Cosine()
    mirror = true

    # horizontal stabilizer
    xle_h = [0.0, 0.14]
    yle_h = [0.0, 1.25]
    zle_h = [0.0, 0.0]
    chord_h = [0.7, 0.42]
    theta_h = [0.0, 0.0]
    phi_h = [0.0, 0.0]
    ns_h = 12
    nc_h = 12
    spacing_s_h = Cosine()
    spacing_c_h = Cosine()
    mirror_h = true

    # vertical stabilizer
    xle_v = [0.0, 0.14]
    yle_v = [0.0, 0.0]
    zle_v = [0.0, 1.0]
    chord_v = [0.7, 0.42]
    theta_v = [0.0, 0.0]
    phi_v = [0.0, 0.0]
    ns_v = 12
    nc_v = 10
    spacing_s_v = Cosine()
    spacing_c_v = Cosine()
    mirror_v = false

    Sref = 9.0
    cref = 0.9
    bref = 10.0
    rref = [0.5, 0.0, 0.0]
    Vinf = 1.0
    ref = Reference(Sref, cref, bref, rref, Vinf)

    alpha = 5.0*pi/180
    beta = 0.0
    Omega = [0.0; 0.0; 0.0]
    fs = VortexLattice.Freestream(Vinf, alpha, beta, Omega)

    symmetric = [false, false, false]

    # generate surface panels for wing
    wgrid, wing = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc;
        mirror=mirror, fc = fc, spacing_s=spacing_s, spacing_c=spacing_c)

    # generate surface panels for horizontal tail
    hgrid, htail = wing_to_surface_panels(xle_h, yle_h, zle_h, chord_h, theta_h, phi_h, ns_h, nc_h;
        mirror=mirror_h, fc=fc_h, spacing_s=spacing_s_h, spacing_c=spacing_c_h)
    translate!(hgrid, [4.0, 0.0, 0.0])
    translate!(htail, [4.0, 0.0, 0.0])

    # generate surface panels for vertical tail
    vgrid, vtail = wing_to_surface_panels(xle_v, yle_v, zle_v, chord_v, theta_v, phi_v, ns_v, nc_v;
        mirror=mirror_v, fc=fc_v, spacing_s=spacing_s_v, spacing_c=spacing_c_v)
    translate!(vgrid, [4.0, 0.0, 0.0])
    translate!(vtail, [4.0, 0.0, 0.0])

    # display(size.([ wing, htail, vtail ])) # Checking sizes

    grids = [wgrid, hgrid, vgrid]
    surfaces = [wing, htail, vtail]
    surface_id = [1, 2, 3]

    system = steady_analysis(surfaces, ref, fs; symmetric=symmetric, surface_id=surface_id)

    CF, CM = body_forces(system; frame=Wind())

    CDiff = far_field_drag(system)

    CF, CM, CDiff
end

max_forks_repo_name: GodotMisogi/AeroMDAO


text:
type Biconnected_components_output
    bcc_edge_list::Vector{Int64}  #list of bcc edges
    components::Vector{Int64} #Indices of bcc corresponding to the edges in List
    number::Int64 # Total number of Biconnected components
    A::MatrixNetwork #MatrixNetwork
end

function BiCC(A::MatrixNetwork)
    n=A.n
    rp=A.rp
    ci=A.ci
    bcc_edges=zeros(Int64,2*length(ci))
    bcc_edges_count=0
    cn=1
    low=zeros(Int64,n)
    dt=zeros(Int64,n)
    pred=zeros(Int64,n)
    components=-1 * ones(Int64,length(ci))
    t=0
    cs=zeros(Int64,2*length(ci))
    css=0
    rs=zeros(Int64,2*n)
    rss=0 #recursion stack (v,ri)
    children = zeros(Int64, n)
    #start dfs at 1
    for sv=1:n
        v=sv
        if dt[v]>0
            continue
        end
        rss=rss+1
        rs[2*rss-1]=v
        rs[2*rss]=rp[v] #add v,rp[v] to stack
        low[v]=t
        dt[v]=t;
        t=t+1;

        while rss>0
            v=rs[2*rss-1]
            ri=rs[2*rss]
            rss=rss-1
            while ri<rp[v+1]
                w=ci[ri]
                ri=ri+1
                if dt[w] == 0
                    low[w]=w; dt[w]=t; t=t+1;
                    css=css+1
                    cs[2*css-1]=v
                    cs[2*css]=w
                    children[v]=children[v]+1
                    pred[w]=v
                    rss=rss+1
                    rs[2*rss-1]=v
                    rs[2*rss]=ri
                    v=w
                    ri=rp[w]
                    continue
                elseif(w != pred[v] && dt[w] < low[v])
 		    if low[v] > dt[w]
			low[v] = dt[w]
		    end
                    css=css+1
                    cs[2*css-1]=v
                    cs[2*css]=w
                end
            end
             for ri=rp[v]:rp[v+1]-1
                w=ci[ri]
                if components[w] == -1
                    if low[v] > low[w]
                        low[v]=low[w]
                    end
                end
                if (( (dt[v] == 1 ) && ( children[v] > 1 ) ) || ( dt[v] > 1 && low[w] >= dt[v]))
                    while css > 0
                        bcc_edges_count = bcc_edges_count + 1
                        bcc_edges[2*bcc_edges_count-1] = cs[2*css-1]
                        bcc_edges[2*bcc_edges_count] = cs[2*css]
                        components[bcc_edges_count] = cn
			if (cs[2*css-1] == v && cs[2*css] ==w)
				css=css-1
				break
			end
                        css = css-1
                    end

                    cn=cn+1
                end
            end

          end
	end
	print(components)
          return bcc_edges
end


text:
#=
Functions to compute chemical potential-related characteristics
=#

"""
    log_c_activity(mixture, nmol, volume, RT[; buf])

Return vector of ln(c_a) - logarithm of activity coefficients
for components of `mixture` at given `nmol`, `volume`, `RT`.
If buffers are provided as keyword arguments, their contents
are modified during the intermediate calculations.

# Arguments

- `mix`: mixture
- `nmol::AbstractVector`: amount of each component (mol)
- `volume`: volume of the mixture (m³)
- `RT`: thermal energy (J/mol)

# Keywords

- `buf::Union{NamedTuple, AbstractDict, BrusilovskyThermoBuffer}`: buffers for intermediate
    calculations (see [`pressure`](@ref))

# Returns

- `AbstractVector`: the logarithms of activity coefficients of the components at given
    number of moles, volume and temperature

See (Jirí Mikyska, Abbas Firoozabadi // 10.1002/aic.12387)

See also: [`log_c_activity!`](@ref), [`log_c_activity_wj`](@ref), [`log_c_activity_wj!`](@ref)
"""
function log_c_activity(
    mix::BrusilovskyEoSMixture{T},
    nmol::AbstractVector,
    volume::Real,
    RT::Real;
    buf = NamedTuple(),
) where {T}

    nc = ncomponents(mix)
    log_ca = Vector{T}(undef, nc)

    return log_c_activity!(log_ca, mix, nmol, volume, RT; buf = buf)
end

"""
    log_c_activity!(log_ca, mixture, nmol, volume, RT[; buf])

Return vector of ln(c_a) - logarithm of activity coefficients
for components of `mixture` at given `nmol`, `volume`, `RT`.
If buffers are provided as keyword arguments, their contents
are modified during the intermediate calculations. The answer is stored
in `log_ca`.

# Arguments

- `log_ca::AbstractVector`: buffer to store the result
- `mix::BrusilovskyEoSMixture`: mixture
- `nmol::AbstractVector`: amount of each component (mol)
- `volume`: volume of the mixture (m³)
- `RT`: thermal energy (J/mol)

# Keywords

- `buf::Union{BrusilovskyThermoBuffer, NamedTuple, AbstractDict}`: buffers for intermediate
    calculations (see [`pressure`](@ref))

# Returns

- `AbstractVector`: the logarithms of activity coefficients of the components at given
    number of moles, volume and temperature

See (Jirí Mikyska, Abbas Firoozabadi // 10.1002/aic.12387)

See also: [`log_c_activity`](@ref), [`log_c_activity_wj`](@ref), [`log_c_activity_wj!`](@ref)
"""
function log_c_activity!(
    log_ca::AbstractVector,
    mix::BrusilovskyEoSMixture,
    nmol::AbstractVector,
    volume,
    RT;
    buf = NamedTuple(),
)
    nc = ncomponents(mix)
    ∑mol = sum(nmol)

    # коэффициенты вычислены для состава в [моль]
    A, B, C, D, aij = eos_parameters(mix, nmol, RT; buf = buf)

    VmB = volume - B
    VpC = volume + C
    VpD = volume + D
    CmD = C - D

    log_VmBbyV = log(VmB/volume)
    ∑molbyVmB = ∑mol / VmB
    AbyRTCmD = A / (RT * CmD)
    log_VpCbyVpD = log(VpC / VpD)

    @inbounds for i in 1:nc
        subst = mix.components[i]
        b = subst.b
        c = subst.c
        d = subst.d

        # ∂A_i = 2 ∑_j (N_j * a_ij)
        # прежний код использовал sum( генератор с for in eachindex ) - давал 2 аллокации
        ∂A = 2 * dot(nmol, @view aij[:,i])

        log_ca[i] = -log_VmBbyV + b * ∑molbyVmB
        brackets = (∂A / A - (c - d) / CmD) * log_VpCbyVpD + (c / VpC - d / VpD)
        log_ca[i] -= AbyRTCmD * brackets
    end
    return log_ca
end

"""
    log_c_activity_wj(mixture, nmol, volume, RT[; buf])

Return vector of ln(c_a) - logarithm of activity coefficient
for components of `mixture` at given `nmol`, `volume`, `RT` -
and the jacobian ∂ln(c_a[i]) / ∂n[j].
If buffers are provided as keyword arguments, their contents
are modified during the intermediate calculations.

# Arguments

- `mix`: mixture
- `nmol::AbstractVector`: amount of each component (mol)
- `volume`: volume of the mixture (m³)
- `RT`: thermal energy (J/mol)

# Keywords

- `buf::BrusilovskyThermoBuffer`: buffers for intermediate calculations
    (see [`thermo_buffer`](@ref))

# Returns

- `Tuple{AbstractVector,AbstractMatrix}`: the logarithms of activity coefficients
of the components at given number of moles, volume and temperature, and the jacobian
matrix ∂ln(c_a[i]) / ∂n[j].

See also: [`log_c_activity`](@ref), [`log_c_activity!`](@ref), [`log_c_activity_wj!`](@ref)
"""
function log_c_activity_wj(
    mix::BrusilovskyEoSMixture{T},
    nmol::AbstractVector,
    volume::Real,
    RT::Real;
    buf::BrusilovskyThermoBuffer = thermo_buffer(mix),
) where {T}
    nc = ncomponents(mix)
    log_ca = Vector{T}(undef, nc)
    jacobian = similar(log_ca, (nc, nc))
    return log_c_activity_wj!(log_ca, jacobian, mix, nmol, volume, RT; buf = buf)
end

"""
    log_activity_wj!(log_ca, jacobian, mixture, nmol, volume, RT[; buf])

Return vector of ln(c_a) - logarithm of activity coefficient
for components of `mixture` at given `nmol`, `volume`, `RT` -
and the jacobian ∂ln(c_a[i]) / ∂n[j]. The first two arguments get overwritten
by the result.
If buffers are provided as keyword arguments, their contents
are modified during the intermediate calculations.

# Arguments

- `log_ca::AbstractVector`: a vector to store the activity coefficients
- `jacobian::AbstractMatrix`: a matrix to store ∂ln(c_a[i]) / ∂n[j]
- `mix`: mixture
- `nmol::AbstractVector`: amount of each component (mol)
- `volume`: volume of the mixture (m³)
- `RT`: thermal energy (J/mol)

# Keywords

- `buf::BrusilovskyThermoBuffer`: buffers for intermediate calculations
    (see [`thermo_buffer`](@ref))

# Returns

- `Tuple{AbstractVector,AbstractMatrix}`: the logarithms of activity coefficients
of the components at given number of moles, volume and temperature, and the jacobian
matrix ∂ln(c_a[i]) / ∂n[j]. The values are aliases for the first two arguments.

See also: [`log_c_activity`](@ref), [`log_c_activity!`](@ref), [`log_c_activity_wj`](@ref)
"""
function log_c_activity_wj!(
    log_ca::AbstractVector,
    jacobian::AbstractMatrix,
    mix::BrusilovskyEoSMixture,
    nmol::AbstractVector,
    volume::Real,
    RT::Real;
    buf::BrusilovskyThermoBuffer = thermo_buffer(mix),
)
    return __log_c_activity_wj_impl__(log_ca, jacobian, mix, nmol, volume, RT, buf)
end

function __log_c_activity_wj_impl__(
    log_ca::AbstractVector,
    jacobian::AbstractMatrix,
    mix::BrusilovskyEoSMixture,
    nmol::AbstractVector,
    volume,
    RT,
    buf::BrusilovskyThermoBuffer,
)
    A, B, C, D, aij = eos_parameters(mix, nmol, RT; buf = buf)

    aux1 = buf.vec1
    aux2 = buf.vec2
    nc = ncomponents(mix)
    V = volume
    ntotal = sum(nmol)

    VmB = volume - B
    VpC = volume + C
    VpD = volume + D
    CmD = C - D

    log_VmBbyV = log(VmB/volume)
    AbyRTCmD = A / (RT * CmD)
    log_VpCbyVpD = log(VpC / VpD)
    log2 = log1p(CmD / (V + D)) / (RT * CmD)

    @inbounds map!(aux1, 1:nc) do i
        mix[i].b / VmB
    end
    log_ca .= aux1 .* ntotal .- log_VmBbyV
    jacobian .= aux1 .+ aux1' .+ ntotal .* aux1 .* aux1'
    mul!(aux1, aij, nmol)
    @inbounds map!(aux2, 1:nc) do i
        mix[i].c - mix[i].d
    end
    log_ca .-= AbyRTCmD .* log_VpCbyVpD .* (2 .* aux1 ./ A .- aux2 ./ CmD)
    jacobian .-= (2 * log2) .* ((A / CmD^2) .* aux2 .* aux2' .+ aij)
    jacobian .+= (2 * log2 / CmD) .* (aux2 .* aux1' .+ aux1 .* aux2')

    aux1 .= 2 .* aux1 ./ (RT * CmD) .- aux2 .* (A / (RT * CmD^2))
    @inbounds map!(aux2, 1:nc) do i
        mix[i].c / VpC - mix[i].d / VpD
    end
    log_ca .-= AbyRTCmD .* aux2
    jacobian .-= aux1 .* aux2'.+ aux2 .* aux1'
    @inbounds map!(aux1, 1:nc) do i
        mix[i].d / VpD
    end
    @inbounds map!(aux2, 1:nc) do i
        mix[i].c / VpC
    end

    jacobian .-= AbyRTCmD .* (aux1 .* aux1' .- aux2 .* aux2')
    return log_ca, jacobian
end

max_stars_repo_name: vvpisarev/CubicEoS

max_issues_repo_name: vvpisarev/CubicEoS

max_forks_repo_name: vvpisarev/CubicEoS


text:
using DataFrames
using CSV

data = CSV.read("data/multilane_Saturday_28_Apr_16_17.csv")

means = by(data, :solver) do df
    n = size(df, 1)
    return DataFrame(reward=mean(df[:reward]),
                     reward_sem=std(df[:reward])/sqrt(n)
                    )
end

println(means)


text:
@enum LambdaFunctionState Pending Active Inactive Failed

struct LambdaException <: Exception
  msg::String
end

"""
    @with_kw mutable struct LambdaFunction
        FunctionName::Union{Missing, String} = missing
        FunctionArn::Union{Missing, String} = missing
        Runtime::Union{Missing, String} = missing
        Role::Union{Missing, String} = missing
        Handler::Union{Missing, String} = missing
        CodeSize::Union{Missing, Int64} = missing
        Description::Union{Missing, String} = missing
        Timeout::Union{Missing, Int64} = missing
        MemorySize::Union{Missing, Int64} = missing
        LastModified::Union{Missing, String} = missing
        CodeSha256::Union{Missing, String} = missing
        Version::Union{Missing, String} = missing
        TracingConfig::Union{Missing, Dict{String, Any}} = missing
        RevisionId::Union{Missing, String} = missing
        PackageType::Union{Missing, String} = missing
        exists::Bool = true
    end

Represents a Lambda function, hosted on AWS. Should not be instantiated directly. If `exists` is
`true`, then the image is assumed to exit and so should be visible from utilities such as `aws
lambda list-functions`
"""
@with_kw mutable struct LambdaFunction <: LambdaComponent
  FunctionName::Union{Missing, String} = missing
  FunctionArn::Union{Missing, String} = missing
  Runtime::Union{Missing, String} = missing
  Role::Union{Missing, String} = missing
  Handler::Union{Missing, String} = missing
  CodeSize::Union{Missing, Int64} = missing
  Description::Union{Missing, String} = missing
  Timeout::Union{Missing, Int64} = missing
  MemorySize::Union{Missing, Int64} = missing
  LastModified::Union{Missing, String} = missing
  CodeSha256::Union{Missing, String} = missing
  Version::Union{Missing, String} = missing
  TracingConfig::Union{Missing, Dict{String, Any}} = missing
  RevisionId::Union{Missing, String} = missing
  PackageType::Union{Missing, String} = missing
  exists::Bool = true
end
StructTypes.StructType(::Type{LambdaFunction}) = StructTypes.Mutable()
Base.:(==)(a::LambdaFunction, b::LambdaFunction) = (a.CodeSha256 == b.CodeSha256)

"""
    get_all_lambda_functions(jot_generated_only::Bool = true)::Vector{LambdaFunction}

Returns a vector of `LambdaFunction`s, representing all AWS-hosted Lambda functions.

`jot_generated_only` specifies whether to filter for jot-generated lambda functions only.
"""
function get_all_lambda_functions(jot_generated_only::Bool = true)::Vector{LambdaFunction}
  all_json = readchomp(`aws lambda list-functions`)
  all_lfs = JSON3.read(all_json, Dict{String, Vector{LambdaFunction}})["Functions"]
  jot_generated_only ? filter(is_jot_generated, all_lfs) : all_lfs
end

"""
    get_lambda_function(function_name::String)::Union{Nothing, LambdaFunction}

Queries AWS and returns a `LambdaFunction` object, representing a Lambda Function hosted on AWS.
"""
function get_lambda_function(function_name::String)::Union{Nothing, LambdaFunction}
  all = get_all_lambda_functions()
  index = findfirst(x -> x.FunctionName == function_name, all)
  isnothing(index) ? nothing : all[index]
end

"""
    get_lambda_function(repo::ECRRepo)::Union{Nothing, LambdaFunction}

Queries AWS and returns a `LambdaFunction` object, representing a Lambda Function hosted on AWS.
The Lambda function returned is based off the given `ECRRepo` instance.
"""
function get_lambda_function(repo::ECRRepo)::Union{Nothing, LambdaFunction}
  all = get_all_lambda_functions()
  index = findfirst(x -> x.FunctionName == function_name, all)
  isnothing(index) ? nothing : all[index]
end

"""
    get_remote_image(lambda_function::LambdaFunction)::RemoteImage

Queries AWS and returns a `RemoteImage` object, representing a docker image hosted on AWS ECR.
The RemoteImage returned provides the code for the provided `lambda_function`.
"""
function get_remote_image(lambda_function::LambdaFunction)::RemoteImage
  remote_images = get_all_remote_images()
  index = findfirst(x -> matches(x, lambda_function), remote_images)
  isnothing(index) && error("Unable to find RemoteImage associated with LambdaFunction")
  remote_images[index]
end

"""
    delete!(func::LambdaFunction)

Deletes a Lambda function hosted on AWS. The LambdaFunction instance continues to exist, but has its
`exists` attribute set to `false`.
"""
function delete!(func::LambdaFunction)
  func.exists || error("Function does not exist")
  delete_script = get_delete_lambda_function_script(func.FunctionArn)
  output = readchomp(`bash -c $delete_script`)
  func.exists = false
  associated_role = get_aws_role(create_role_name(func.FunctionName))
  !isnothing(associated_role) && delete!(associated_role)
  nothing
end

function get_function_state(func_name::String)::LambdaFunctionState
  get_state_script = get_lambda_function_status(func_name)
  state_json = readchomp(`bash -c $get_state_script`)
  state_data = JSON3.read(state_json)
  if state_data["State"] == "Pending" Pending
  elseif state_data["State"] == "Active" Active
  end
end

function get_function_state(func::LambdaFunction)::LambdaFunctionState
  get_function_state(func.FunctionArn)
end

"""
    invoke_function(
        request::Any,
        lambda_function::LambdaFunction;
        check_state::Bool=false,
      )::Any

Invokes a Lambda function, hosted on AWS. `request` is the argument that it will be called with.
This will be automatically converted to JSON before sending, so should match the
`response_function_param_type` of the responder used to create the function.

Returns the invoked Lambda function response, or throws an error if the invoked Lambda function has
returned an error status.

If `check_state` is `true`, the function will wait for the AWS Lambda function to become available
before sending the request. This can be useful if the Lambda function has been created within the
last few seconds, since there is a short set-up time before it can be called.
"""
function invoke_function(
    request::Any,
    lambda_function::LambdaFunction;
    check_state::Bool=false,
  )::Any
  if check_state
    while true
      Jot.get_function_state(lambda_function) == Active && break
    end
  end
  request_json = JSON3.write(request)
  outfile_path = tempname()
  invoke_script = get_invoke_lambda_function_script(lambda_function.FunctionArn,
                                                    request_json,
                                                    outfile_path)
  status = readchomp(`bash -c $invoke_script`) |> JSON3.read
  response = open(outfile_path, "r") do f
    read(f, String) |> JSON3.read
  end
  if haskey(status, "FunctionError")
    throw(LambdaException("$response"))
  else
    response
  end
end


text:
# test backup_ancestors and backup_descendants
let
    NODE_TYPE=Int
    ID_TYPE=Symbol
    my_tree = GraphUtils.NTree{NODE_TYPE,ID_TYPE}()
    add_node!(my_tree,1,:ONE)
    add_node!(my_tree,2,:TWO)
    add_node!(my_tree,3,:THREE)
    @test add_edge!(my_tree,:ONE,:TWO)
    @test add_edge!(my_tree,:TWO,:THREE)

    d_map = backup_descendants(my_tree,n->node_id(n)==:TWO)
    @test d_map[:THREE] === nothing
    @test d_map[:TWO] == :TWO
    @test d_map[:ONE] == :TWO
    d_map = GraphUtils.backup_multiple_descendants(my_tree,n->node_id(n)==:TWO)
    @test isempty(d_map[:THREE])
    @test :TWO in d_map[:TWO]
    @test :TWO in d_map[:ONE]
    d_map = backup_ancestors(my_tree,n->node_id(n)==:TWO)
    @test d_map[:THREE] == :TWO
    @test d_map[:TWO] == :TWO
    @test d_map[:ONE] === nothing
    d_map = GraphUtils.backup_multiple_ancestors(my_tree,n->node_id(n)==:TWO)
    @test :TWO in d_map[:THREE]
    @test :TWO in d_map[:TWO]
    @test isempty(d_map[:ONE])

end
let
    NODE_TYPE=Int
    ID_TYPE=Symbol
    my_tree = GraphUtils.NTree{NODE_TYPE,ID_TYPE}()
    add_node!(my_tree,1,:ONE)
    add_node!(my_tree,2,:TWO)
    add_node!(my_tree,3,:THREE)
    @test add_edge!(my_tree,:ONE,:TWO)
    # Vertex :TWO already has a parent--can't have another.
    @test !add_edge!(my_tree,:THREE,:TWO)
    add_child!(my_tree,:THREE,4,:FOUR)
    @test has_edge(my_tree,:THREE,:FOUR)

    d_map = backup_descendants(my_tree,n->node_id(n)==:FOUR)
    @test d_map[:FOUR] == :FOUR
    @test d_map[:THREE] == :FOUR
    @test d_map[:TWO] === nothing
    @test d_map[:ONE] === nothing

    d_map = backup_descendants(my_tree,n->node_id(n)==:THREE)
    @test d_map[:FOUR] === nothing
    @test d_map[:THREE] == :THREE
    @test d_map[:TWO] === nothing
    @test d_map[:ONE] === nothing

    d_map = backup_ancestors(my_tree,n->node_id(n)==:THREE)
    @test d_map[:FOUR] === :THREE
    @test d_map[:THREE] == :THREE
    @test d_map[:TWO] === nothing
    @test d_map[:ONE] === nothing
end
let
    NODE_TYPE = Int
    ID_TYPE = Symbol
    EDGE_TYPE = GraphUtils.CustomEdge{Int,ID_TYPE}
    g = GraphUtils.CustomNEGraph{DiGraph,NODE_TYPE,EDGE_TYPE,ID_TYPE}()
    e = GraphUtils.CustomEdge(:ONE,:TWO,1+2)
    e2 = GraphUtils.CustomEdge(:TWO,:THREE,2+3)
    # EDGE_TYPE = Int
    # g = GraphUtils.NEGraph{DiGraph,NODE_TYPE,EDGE_TYPE,ID_TYPE}()
    # e = 3
    # e2 = 5
    isa(g,GraphUtils.AbstractCustomNGraph)
    add_node!(g,1,:ONE)
    add_node!(g,2,:TWO)
    add_node!(g,3,:THREE)
    @test add_edge!(g,:ONE,:TWO,e)
    @test !has_edge(g,:TWO,:ONE)
    @test add_edge!(g,:TWO,:THREE,e2)
    rem_node!(g,:THREE)
    @test !has_vertex(g,:THREE)
    @test !has_edge(g,:TWO,:THREE)
    @test !has_edge(g,e2)

    make_edge(g,:TWO,:ONE,-1)
    GraphUtils.make_edge(g,u::Symbol,v::Symbol) = GraphUtils.make_edge(g,u,v,-1)
    @test add_edge!(g,:TWO,:ONE)
    @test edge_val(get_edge(g,:TWO,:ONE)) == -1
    rem_edge!(g,:TWO,:ONE)
    @test !has_edge(g,:TWO,:ONE)
    @test add_edge!(g,:TWO,:ONE,-2)
    @test edge_val(get_edge(g,:TWO,:ONE)) == -2
end
let
    G = DiGraph(3)
    add_edge!(G,1,2)
    add_edge!(G,2,3)
    @test depth_first_search(G,1,
        v->v==3,
        v->true,
        outneighbors
    ) == 3
    @test depth_first_search(G,1,
        v->v==3,
        v->false,
        outneighbors
    ) == -1
    @test depth_first_search(G,3,
        v->v==1,
        v->true,
        inneighbors
    ) == 1
    @test depth_first_search(G,1,
        v->v>=1,
        v->true,
        outneighbors
    ) == 1
    @test depth_first_search(G,1,
        v->v>=1,
        v->true,
        outneighbors;
        skip_first=true
    ) == 2

    @test collect_ancestors(G,2) == Set(1)
    @test collect_ancestors(G,3) == Set([1,2])
    @test collect_descendants(G,2) == Set(3)
    @test collect_descendants(G,1) == Set([2,3])
    @test collect_descendants(G,[1]) == Set([2,3])

    keep = true
    @test collect_descendants(G,[1],keep) == Set([1,2,3])

end
let
    graph = DiGraph(7)
    add_edge!(graph,1,3)
    add_edge!(graph,2,3)
    add_edge!(graph,4,6)
    add_edge!(graph,5,6)
    add_edge!(graph,3,7)
    add_edge!(graph,6,7)
    for (g,target) in [
        (graph,[1,1,2,1,1,2,3]),
        (reverse(graph),[3,3,2,3,3,2,1])
    ]
        for iter in [
                GraphUtils.BFSIterator(g),
                GraphUtils.SortedBFSIterator(g),
                ]
            depths = ones(Int,nv(g))
            while !isempty(iter)
                v = pop!(iter)
                if indegree(g,v) > 0
                    depths[v] = max(depths[v],1+maximum(depths[vp] for vp in inneighbors(g,v)))
                end
                GraphUtils.update_iterator!(iter,v)
            end
            @test all(depths .== target)
        end
    end
    @test length(collect(BFSIterator(graph))) == nv(graph)
    @test collect(SortedBFSIterator(graph)) == [1, 2, 4, 5, 3, 6, 7]

end

let
    G = NGraph{DiGraph,String,String}()
    for s in ["A","B","C","D","E"]
        add_node!(G,s,s)
    end
    edge_set = Set([("A","B"),("B","C"),("C","D"),("B","E")])
    for e in edge_set
        add_edge!(G,e[1],e[2])
    end
    for v in [3,2,1,1,1]
        rem_node!(G,v)
        for e in edges(G)
            @test (
                node_val(get_node(G,e.src)),
                node_val(get_node(G,e.dst))) in edge_set
        end
    end

end


text:
#!/bin/julia
#
# An implementation of the lattice boltzmann method
#
# @author arbennett
##

# Constant parameters
const T   = 9000          # Number of timesteps to run
const Re  = 20.0           # Reynolds number
const nx  = 150              # Size of the domain in x direction
const ny  = 50             # Size of the domain in y direction
const q   = 9               # number of flow directions
const v   = 0.05            # velocity
const nu  = (v*ny)/(10*Re)  # viscosity
const tau = (3.0*nu+0.5)    # 1/relaxation
const plot_frequency = T/5

## Set up the flow grids, weights, and indices
#
#  Explain here
#
##
dirs =     [  0   1   0  -1   0   1   -1   -1    1 ;
              0   0   1   0  -1   1    1   -1   -1 ]
weights =  [4/9 1/9 1/9 1/9 1/9 1/36 1/36 1/36 1/36]
indices =  [  4   5   6;  1   2   3;   7    8    9 ]
no_slip =  [  1   4   5   2   3   8    9    6    7 ]

## Create the obstacle array
#
#  The obst array is a boolean array that uses true to denote where
#  obstacles to flow are.  When the "fluid" encounters one of these
#  cells it will reverse the flow direction using the `no_slip`
#  condition detailed above
#
##
obst_coords = [nx/2, ny/2, 3]  # [x, y, radius]
obst_coords2 = [nx/3,ny/3,4]
y = collect(1:ny)
x = collect(1:nx)'
obst = (x.-obst_coords[1]).^2 .+ (y.-obst_coords[2]).^2 .<= obst_coords[3].^2
obst = (obst + ((x.-obst_coords2[1]).^2 .+ (y.-obst_coords2[2]).^2 .<= obst_coords2[3].^2) .> 0)
#obst[1,:] = true
#obst[end,:] = true
#obst[:,:] = false

##  Set up initial flow profile
#
#  Allocate the velocity array, u, which has 2 points for each grid
#  point, corresponding to the x and y velocity
#
#  Allocate the distribution array, f, which has 9 points for each grid
#  point, corresponding to the directions in dirs
##
u = zeros(ny,nx,2)
vel = zeros(ny,nx,2)
fEq = zeros(ny,nx,9)
rho = ones(ny,nx)

for j=1:ny, i=1:nx
    vel[j,i,2] = v*(1+0.1*sin((j/ny)/(2*pi)))
end

## Calculate equilibrium for some distribution
#
##
function equilib(u, rho)
    dist = zeros(ny,nx,9)
    for i = 1:9
        for x = 1:nx
            for y = 1:ny
                tmp = u[y,x,1] * dirs[1,i] + u[y,x,2] * dirs[2,i]
                dist[y,x,i] = rho[y,x] * weights[i] * (1 + 3*tmp + 4.5*tmp*tmp - 1.5*(u[y,x,1]*u[y,x,1] + u[y,x,2]*u[y,x,2]))
            end
        end
    end
    return dist
end
fIn = equilib(u, rho)
fLeft = fIn[:,1,:]

## Propogate time
#
##
for t=1:T
    if mod(t, plot_frequency) == 0
        println(t)
        writedlm("out/solution.dat."*lpad(string(t),length(string(int(T))),"0"),sqrt(u[:,:,1].*u[:,:,1] + u[:,:,2].*u[:,:,2]))
    end

    # Right wall conditions
    fIn[:,end,:] = fIn[:,end-1,:]
    rho = fIn[:,:,1] + fIn[:,:,2] + fIn[:,:,3] + fIn[:,:,4] + fIn[:,:,5] + fIn[:,:,6] + fIn[:,:,7] + fIn[:,:,8] + fIn[:,:,9]

    u[:,:,1] = (fIn[:,:,2] - fIn[:,:,4] + fIn[:,:,6] - fIn[:,:,7] - fIn[:,:,8] + fIn[:,:,9]) ./ rho[:,:]
    u[:,:,2] = (fIn[:,:,3] - fIn[:,:,5] + fIn[:,:,6] + fIn[:,:,7] - fIn[:,:,8] - fIn[:,:,9]) ./ rho[:,:]

    # Left wall conditions
    u[:,1,:] = vel[:,1,:]
    #rho[:,1] = 1./(1-u[:,1,2]) .* (fIn[:,1,3] + fIn[:,1,6] + fIn[:,1,7] + 2*(fIn[:,1,5] + fIn[:,1,8] + fIn[:,1,9]))
    #println(rho[:,1])

    fEq = fIn - (1/tau) * (fIn - equilib(u, rho))

    for i=1:q
        for x=1:nx
            for y=1:ny
                if obst[y,x]
                    fEq[y,x,i] = fIn[y, x, no_slip[i]]
                    u[y,x,:] = 0
                end
            end
        end
    end

    for i=1:q
        fIn[:,:,i] = circshift(fEq[:,:,i],dirs[:,i])
    end
end
println("done!")

max_stars_repo_name: arbennett/Simple-Simulations-In-Julia

max_issues_repo_name: arbennett/Simple-Simulations-In-Julia

max_forks_repo_name: arbennett/Simple-Simulations-In-Julia


text:
"""
Temperature correction parameters
"""
abstract type AbstractTemperatureCorrection{T} end

# Temperature mixins
@mix @columns struct Tbase{T}
    # Field       | Default | Unit | Bounds          | Log | Description
    reftemp::T    | 300.0   | K    | (273.0, 325.0)  | _   | "Reference temperature for all rate parameters"
    arrtemp::T    | 2000.0  | K    | (200.0, 4000.0) | _   | "Arrhenius temperature"
end

@mix @columns struct Tlow{T}
    tbelow::T     | -30.0   | K | (0.0, -40.0)       | _ | "Lower boundary of tolerance range"
    arrlower::T   | 20000.0 | K | (2000.0, 40000.0)  | _ | "Arrhenius temperature for lower boundary"
end

@mix @columns struct Tup{T}
    tabove::T     | 5.0     | K | (0.0, 20)          | _ | "Upper boundary of tolerance range"
    arrupper::T   | 70000.0 | K | (7000.0, 140000.0) | _ | "Arrhenius temperature for upper boundary"
end

"""
    tempcorr(f, t)

Temperature related correction of growth rate.

-`f`: a formulation struct or `nothing` for no temperature correction.
-`t`: the current temperature, in degrees Celcius or Kelvin (Unitful.jl u"°C" or u"K")
"""
function tempcorr end
tempcorr(f::Nothing, t) = 1

"""
    TempCorr(reftemp, arrtemp)

Simple temperature correction parameters

$(FIELDDOCTABLE)
"""
@Tbase struct TempCorr{T} <: AbstractTemperatureCorrection{T} end

tempcorr(f::TempCorr, t) = tempcorr(t |> K, f.reftemp, f.arrtemp)

"""
    TempCorrLower(reftemp, arrtemp, tbelow, arrlower)

Temperature correction with lower bounds parameters.

$(FIELDDOCTABLE)
"""
@Tbase @Tlow struct TempCorrLower{T} <: AbstractTemperatureCorrection{T} end

tempcorr(f::TempCorrLower, t) =
    tempcorr(t |> K, f.reftemp, f.arrtemp, f.tbelow + f.reftemp, f.arrlower)

"""
    TempCorrLowerUpper(reftemp, arrtemp, tbelow, arrlower, tabove, arrupper)

Temperature correction with lower and upper bound parameters.

$(FIELDDOCTABLE)
"""
@Tbase @Tlow @Tup struct TempCorrLowerUpper{T} <: AbstractTemperatureCorrection{T} end

tempcorr(f::TempCorrLowerUpper, t) =
    tempcorr(t |> K, f.reftemp, f.arrtemp, f.tbelow + f.reftemp, f.arrlower, f.tabove + f.reftemp, f.arrupper)

tempcorr(t, tref, a) = exp(a/tref - a/t)
tempcorr(t, tref, a, l, al) =
    exp(a/tref - a/t) * (1.0 + exp(al/tref - al/l)) / (1.0 + exp(al/t - al/l))
tempcorr(t, tref, a, l, al, h, ah) =
    exp(a/tref - a/t) * (1.0 + exp(al/tref - al/l) + exp(ah/h - ah/tref)) /
    (1.0 + exp(al/t - al/l) + exp(ah/h - ah/t))

"""
    ParentTardieu(ΔH_A, α, t0)

Simple 3 parameter temperature correction method.
Growth response to temperature has smoother transients in plants than in animals,
and a simpler formulaiton is more applicable.

$(FIELDDOCTABLE)
"""
@flattenable @columns struct ParentTardieu{Δ,Al,T,A} <: AbstractTemperatureCorrection{T}
    ΔH_A::Δ | true  | 63.5  | kJ*mol^-1 | (55.0, 65.0)   | _ | "The enthalpy of activation of the reaction. Determines the curvature at low temperature"
    α::Al   | true  | 3.5   | _         | (1.0, 10.0)    | _ | "The ratio ΔH_D / ΔH_A"
    t0::T   | true  | 300.0 | K         | (273.0, 325.0) | _ | "Reference temperature"
    A::A    | false | 0.0   | K^-1      | (0.0, 1.0)     | _ | "Trait scaling coefficient, calculated from other params"
end
# Allways recalculate A
ParentTardieu(ΔH_A, α, t0, A) = ParentTardieu(ΔH_A, α, t0)
# Precalculate scalar A using the reference temperature
# to normalise the temperature response with the reference at 1.0
ParentTardieu(ΔH_A, α, t0) = begin
    A = 1 / parent_tardieua_unscaled(ΔH_A, α, t0, t0)
    ParentTardieu{typeof.((ΔH_A, α, t0, A))...}(ΔH_A, α, t0, A)
end

parent_tardieua_unscaled(ΔH_A, α, t0, t) = begin
    ex = exp(-ΔH_A / (R * t))
    t * ex / (1 + ex^(α * (1 - (t / t0))))
end

tempcorr(f::ParentTardieu, t) = f.A * parent_tardieua_unscaled(f.ΔH_A, f.α, f.t0, t)


text: stripChar = (s, r) -> replace(s, Regex("[$r]") => "")

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
type Axis
	axis::RectangleShape
	marks::Array{RectangleShape}
	numbers::Array{RenderText}
	xaxis::Bool
	number::RenderText
	tic::RectangleShape
	tic_period::Real
end


text:
import Base: rand, length
import Distributions: logpdf, AbstractMvNormal

export elbo, dtc
export SparseFiniteGP

"""
    SparseFiniteGP{T1<:FiniteGP, T2<:FiniteGP}

A finite-dimensional projection of an `AbstractGP` `f` at locations `x`, which uses a second
`FiniteGP` defined at a sparse set of inducing points [1] to do approximate inference.

This object has similar methods to an ordinary `FiniteGP`, but note that the methods for
`logpdf` and `←` are just convenience wrappers around `elbo` and `PseudoObs`.

[1] - M. K. Titsias. "Variational learning of inducing variables in sparse Gaussian
processes". In: Proceedings of the Twelfth International Conference on Artificial
Intelligence and Statistics. 2009.

```jldoctest
julia> f = GP(Matern32(), GPC());

julia> fobs = f(rand(100));

julia> finducing = f(0:0.1:1);

julia> fxu = SparseFiniteGP(fobs, finducing);

julia> y = rand(fxu);

julia> logpdf(fxu, y) < logpdf(fobs, y)
true

julia> fpost = f | (fxu ← y);
```
"""
struct SparseFiniteGP{T1<:FiniteGP, T2<:FiniteGP} <: AbstractMvNormal
    fobs::T1
    finducing::T2
end

length(f::SparseFiniteGP) = length(f.fobs)
mean(f::SparseFiniteGP) = mean(f.fobs)
covariance_error = "The covariance matrix of a sparse GP can often be dense and " *
    "can cause the computer to run out of memory. If you are sure you have enough " *
    "memory, you can use `cov(f.fobs)`."
cov(f::SparseFiniteGP) = error(covariance_error) #cov(f.fobs)
# Not sure how to implement the following...
# cov(fx::FiniteGP, gx::FiniteGP) = cov(fx.f, gx.f, fx.x, gx.x)
marginals(f::SparseFiniteGP) = marginals(f.fobs)

rand(rng::AbstractRNG, f::SparseFiniteGP, N::Int) = rand(rng, f.fobs, N)
rand(f::SparseFiniteGP, N::Int) = rand(Random.GLOBAL_RNG, f, N)
rand(rng::AbstractRNG, f::SparseFiniteGP) = vec(rand(rng, f, 1))
rand(f::SparseFiniteGP) = vec(rand(f, 1))

elbo(f::SparseFiniteGP, y::AV{<:Real}) = elbo(f.fobs, y, f.finducing)
logpdf(f::SparseFiniteGP, y::AV{<:Real}) = elbo(f, y)
logpdf(f::SparseFiniteGP, Y::AbstractMatrix{<:Real}) = map(y -> logpdf(f, y), eachcol(Y))

PseudoObs(fxu::SparseFiniteGP, y) = PseudoObs(Obs(fxu.fobs, y), fxu.finducing)
←(fxu::SparseFiniteGP, y) = PseudoObs(fxu, y)


text:
## Layer types
"""Union type for convolutional layers."""
const ConvLayer = Union{Conv} # TODO: DepthwiseConv, ConvTranspose, CrossCor

"""Union type for dropout layers."""
const DropoutLayer = Union{Dropout,typeof(Flux.dropout),AlphaDropout}

"""Union type for reshaping layers such as `flatten`."""
const ReshapingLayer = Union{typeof(Flux.flatten)}

"""Union type for max pooling layers."""
const MaxPoolLayer = Union{MaxPool,AdaptiveMaxPool,GlobalMaxPool}

"""Union type for mean pooling layers."""
const MeanPoolLayer = Union{MeanPool,AdaptiveMeanPool,GlobalMeanPool}

"""Union type for pooling layers."""
const PoolingLayer = Union{MaxPoolLayer,MeanPoolLayer}

# Activation functions
"""Union type for ReLU-like activation functions."""
const ReluLikeActivation = Union{typeof(relu),typeof(gelu),typeof(swish),typeof(mish)}

"""Union type for softmax activation functions."""
const SoftmaxActivation = Union{typeof(softmax),typeof(softmax!)}

# Layers & activation functions supported by LRP
"""Union type for layers that are allowed by default in "deep rectifier networks"."""
const LRPSupportedLayer = Union{Dense,ConvLayer,DropoutLayer,ReshapingLayer,PoolingLayer}

"""Union type for activation functions that are allowed by default in "deep rectifier networks"."""
const LRPSupportedActivation = Union{typeof(identity),ReluLikeActivation}


text:
# Box type and related methods

struct Box <: Shape
    id::Int64
    obj_to_world::Transformation
    world_to_obj::Transformation
end

Box() = Box(1, Transformation(), Transformation())
Box(T::Transformation) = Box(1, T, inv(T))
Box(id::Int64, T::Transformation) = Box(id, T, inv(T))

can_intersect(b::Box) = true

function (T::Transformation)(B::Box)
    T2 = T * B.obj_to_world
    Box(B.id, T2, inv(T2))
end

obj_bounds(b::Box) = BoundingBox(Point3(-1,-1,-1), Point3(1,1,1))

function world_bounds(B::Box)
    B.obj_to_world(obj_bounds(B))
end

function shape_intersect(r::Ray, B::Box)
    ray = B.world_to_obj(r)

    t0 = ray.tmin
    t1 = ray.tmax
    ibest = -1
    for i = 1:3
        tNear = (-1.0 - ray.origin[i]) / ray.direction[i]
        tFar = (1.0 - ray.origin[i]) / ray.direction[i]
        tNear, tFar = min(tNear, tFar), max(tNear, tFar)
        if tNear > t0
            t0 = tNear
            ibest = i
        end
        t1 = tFar < t1 ? tFar : t1
        if t0 > t1
            # Cannot hit the box
            return nothing, NaN, NaN
        end
    end
    # We hit the box and (t0, t1) are the distances where this happens.
    # These are also already guaranteed to be within (ray.tmin, ray.tmax)

    t = t0
    P = ray(t)
    if ibest == 1
        n = Normal3(1,0,0)
        s = Vector3(0,1,0)
    elseif ibest == 2
        n = Normal3(0,1,0)
        s = Vector3(0,0,1)
    elseif ibest == 3
        n = Normal3(0,0,1)
        s = Vector3(1,0,0)
    end

    n = flip_normal(n, ray.direction)
    DG = DifferentialGeometry(B.obj_to_world(P), B.obj_to_world(n), B.obj_to_world(s))
    return DG, t, 5e-4 * t
end

# Quick intersect function
function intersectP(r::Ray, B::Box)
    ray = B.world_to_obj(r)

    t0 = ray.tmin
    t1 = ray.tmax
    for i = 1:3
        tNear = (-1.0 - ray.origin[i]) / ray.direction[i]
        tFar = (1.0 - ray.origin[i]) / ray.direction[i]
        tNear, tFar = min(tNear, tFar), max(tNear, tFar)
        t0 = tNear > t0 ? tNear : t0
        t1 = tFar < t1 ? tFar : t1
        if t0 > t1
            # Cannot hit the box
            return false
        end
    end
    return true
end

max_stars_repo_name: dronir/GRAYTR

max_issues_repo_name: dronir/GRAYTR

max_forks_repo_name: dronir/GRAYTR


text:
function update_global_vars!(sys::PSY.System, x::AbstractArray)
    index = PSY.get_ext(sys)[GLOBAL_VARS][:ω_sys_index]
    index == 0 && return
    PSY.get_ext(sys)[GLOBAL_VARS][:ω_sys] = x[index]
    return
end

function system!(out::Vector{<:Real}, dx, x, sys::PSY.System, t::Float64)

    I_injections_r = PSY.get_ext(sys)[AUX_ARRAYS][1]
    I_injections_i = PSY.get_ext(sys)[AUX_ARRAYS][2]
    injection_ode = PSY.get_ext(sys)[AUX_ARRAYS][3]
    branches_ode = PSY.get_ext(sys)[AUX_ARRAYS][4]

    #Index Setup
    bus_size = get_bus_count(sys)
    bus_vars_count = 2 * bus_size
    bus_range = 1:bus_vars_count
    injection_start = get_injection_pointer(sys)
    injection_count = 1
    branches_start = get_branches_pointer(sys)
    branches_count = 1
    update_global_vars!(sys, x)

    #Network quantities
    V_r = @view x[1:bus_size]
    V_i = @view x[(bus_size + 1):bus_vars_count]
    Sbase = PSY.get_base_power(sys)
    fill!(I_injections_r, 0.0)
    fill!(I_injections_i, 0.0)

    for d in PSY.get_components(PSY.DynamicInjection, sys)
        bus_n = PSY.get_number(PSY.get_bus(d))
        bus_ix = PSY.get_ext(sys)[LOOKUP][bus_n]
        n_states = PSY.get_n_states(d)
        ix_range = range(injection_start, length = n_states)
        ode_range = range(injection_count, length = n_states)
        injection_count = injection_count + n_states
        injection_start = injection_start + n_states
        device!(
            x,
            injection_ode,
            view(V_r, bus_ix),
            view(V_i, bus_ix),
            view(I_injections_r, bus_ix),
            view(I_injections_i, bus_ix),
            ix_range,
            ode_range,
            d,
            sys,
        )
        out[ix_range] = injection_ode[ode_range] - dx[ix_range]
    end

    for d in PSY.get_components(PSY.ElectricLoad, sys)
        bus_n = PSY.get_number(PSY.get_bus(d))
        bus_ix = PSY.get_ext(sys)[LOOKUP][bus_n]
        device!(
            view(V_r, bus_ix),
            view(V_i, bus_ix),
            view(I_injections_r, bus_ix),
            view(I_injections_i, bus_ix),
            d,
            sys,
        )
    end

    for d in PSY.get_components(PSY.Source, sys)
        bus_n = PSY.get_number(PSY.get_bus(d))
        bus_ix = PSY.get_ext(sys)[LOOKUP][bus_n]
        device!(
            view(V_r, bus_ix),
            view(V_i, bus_ix),
            view(I_injections_r, bus_ix),
            view(I_injections_i, bus_ix),
            d,
            sys,
        )
    end

    if PSY.get_ext(sys)[DYN_LINES]
        dyn_branches = PSY.get_components(PSY.DynamicBranch, sys)
        for br in dyn_branches
            arc = PSY.get_arc(br)
            n_states = PSY.get_n_states(br)
            from_bus_number = PSY.get_number(arc.from)
            to_bus_number = PSY.get_number(arc.to)
            bus_ix_from = PSY.get_ext(sys)[LOOKUP][from_bus_number]
            bus_ix_to = PSY.get_ext(sys)[LOOKUP][to_bus_number]
            ix_range = range(branches_start, length = n_states)
            ode_range = range(branches_count, length = n_states)
            branches_count = branches_count + n_states
            branch!(
                x,
                dx,
                branches_ode,
                #Get Voltage data
                view(V_r, bus_ix_from),
                view(V_i, bus_ix_from),
                view(V_r, bus_ix_to),
                view(V_i, bus_ix_to),
                #Get Current data
                view(I_injections_r, bus_ix_from),
                view(I_injections_i, bus_ix_from),
                view(I_injections_r, bus_ix_to),
                view(I_injections_i, bus_ix_to),
                ix_range,
                ode_range,
                br,
                sys,
            )
            out[ix_range] = branches_ode[ode_range] - dx[ix_range]
        end
    end

    kirchoff_laws!(sys, V_r, V_i, I_injections_r, I_injections_i, dx)
    out[bus_range] = PSY.get_ext(sys)[AUX_ARRAYS][6]

end


text:
# Find and add MDP's path to Julia search path
cdir = pwd(); ii = findlast("MDP", cdir);
if ii===nothing
    display("MDP's path is not found. Please uncomment the above line and set the path to MDP.");
else
    MDPpath = cdir[1:ii[end]] * "/";
    include(MDPpath * "Installation/setpath.jl");
    push!(LOAD_PATH, pwd());
end

using Revise, DelimitedFiles, Preprocessing, Potential, SnapCpp

app = Preprocessing.initializeapp(cdir,MDPpath);

app.unitstyle = "metal"; # unit system
app.pbc = [1 1 1];       # periodic boundary conditions
app.appname = "lj1";

# data input file
app.datapath = "/Users/ngoccuongnguyen/Dropbox (MIT)/Research/Software/MDP/examples/data";
app.datafile = "tiltbox.config"; # configuration file
app.dataformat = "mdp";  # mdp, lammps, json, extxyz
app.datamode = 1;      # 0 => binary, 1 => text
config = Preprocessing.readconfigfile(app);

# # parametric potential descriptors
# epsilon = 0.01029849;
# sigma = 3.4;
# A = 4*epsilon*sigma^12;
# B = 4*epsilon*sigma^6;
# rcut = 8.5;

# potential = Array{Any}(undef, 1)
# potential[1] = Potential.addpotential("pair", "nonbonded", [0], "ljpotential", [A; B], rcut);
# for i = 1:length(potential)
#     include(potential[i].potentialfunction * ".jl");  # include the potential file
# end

# rcutmax = rcut; eta = [0]; kappa = [0];
# etot, e, f = Potential.empiricalpotential(config.x, config.q, config.t, config.a,
#                         config.b, config.c, config.pbc, rcutmax, eta, kappa, potential);

libpath = "cpuSnap.dylib"
SnapCpp.loadsnap(libpath)

twojmax = Int32(6)
ntypes = Int32(1)
ncoeffall = Int32(31)
rcutfac = 4.67637;
rfac0 = 0.99363;
rmin0 = 0;
bzeroflag = 0;
quadraticflag = 0;
switchflag = 1;
chemflag = 0;
bnormflag = 0;
wselfallflag = 0;

snaparam = [ntypes, twojmax, rcutfac, rfac0, rmin0, bzeroflag, switchflag, quadraticflag, chemflag, bnormflag, wselfallflag]
elemradius = 0.5*ones(Float64, ntypes)
elemweight = ones(Float64, ntypes)
snapcoeff = zeros(Float64, ntypes*ncoeffall)

sna = SnapCpp.initsna(snaparam, elemradius, elemweight, snapcoeff)
e, f, v = SnapCpp.snappotential(config.x, config.t, config.a, config.b, config.c, config.pbc, sna)

sna = SnapCpp.initsna(snaparam, elemradius, elemweight)
bi, bd, bv = SnapCpp.snapdescriptors(config.x, config.t, config.a, config.b, config.c, config.pbc, sna)

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
function write(io::IO, x::Any, verbose=false)
    e = Encoder(io, verbose)
    encode_top_level(e, x)
end

function to_transit(x::Any, verbose=false)
    let buf = IOBuffer()
        write(buf, x, verbose)
        s = takebuf_string(buf)
	s
    end
end


text:
#Contains all the evaluation stuff
using JLD2
using PyPlot
using FileIO
using Statistics

a = load("Results/Videos/11_03_ballbot/1_all.jld2")["result"]
b = load("Results/Videos/11_03_ballbot/2_all.jld2")["result"]
c = load("Results/Videos/11_03_ballbot/3_all.jld2")["result"]
d = load("Results/Videos/11_03_ballbot/4_all.jld2")["result"]
e = load("Results/Videos/11_03_ballbot/5_all.jld2")["result"]
#d = load("Results/Videos/10_31_swingup/1_all.jld2")["result"]
#e = load("Results/Videos/10_31_swingup/1_all.jld2")["result"]

a_all = 4
iters = 500

RE  = zeros(a_all,iters)
A2C = zeros(a_all,iters)
TNPG= zeros(a_all,iters)
TRPO= zeros(a_all,iters)
PPO = zeros(a_all,iters)
TD3 = zeros(a_all,iters)

RE[1,:]   = a[1,:]
A2C[1,:]  = a[2,:]
TNPG[1,:] = a[3,:]
TRPO[1,:] = a[4,:]
PPO[1,:]  = a[5,:]
TD3[1,:]  = a[6,:]

RE[2,:]   = b[1,:]
A2C[2,:]  = b[2,:]
TNPG[2,:] = b[3,:]
TRPO[2,:] = b[4,:]
PPO[2,:]  = b[5,:]
TD3[2,:]  = b[6,:]

RE[3,:]   = c[1,:]
A2C[3,:]  = c[2,:]
TNPG[3,:] = c[3,:]
TRPO[3,:] = c[4,:]
PPO[3,:]  = c[5,:]
TD3[3,:]  = c[6,:]

RE[4,:]   = d[1,:]
A2C[4,:]  = d[2,:]
TNPG[4,:] = d[3,:]
TRPO[4,:] = d[4,:]
PPO[4,:]  = d[5,:]
TD3[4,:]  = d[6,:]

#=
RE[5,:]   = e[1,:]
A2C[5,:]  = e[2,:]
TNPG[5,:] = e[3,:]
TRPO[5,:] = e[4,:]
PPO[5,:]  = e[5,:]
TD3[5,:]  = e[6,:]
=#

scale  = 1
vis    = 0.5
offset = 0

re_mean = vec(mean(RE,dims = 1)) .+ offset
re_std  = scale*vec(std(RE,dims=1))

a2c_mean = vec(mean(A2C,dims = 1)) .+ offset
a2c_std  = scale*vec(std(A2C,dims=1))

tnpg_mean = vec(mean(TNPG,dims = 1)) .+ offset
tnpg_std  = scale*vec(std(TNPG,dims=1))

trpo_mean = vec(mean(TRPO,dims = 1)) .+ offset
trpo_std  = scale*vec(std(TRPO,dims=1))

ppo_mean = vec(mean(PPO,dims = 1)) .+ offset
ppo_std  = scale*vec(std(PPO,dims=1))

td3_mean = vec(mean(TD3,dims = 1)) .+ offset
td3_std  = scale*vec(std(TD3,dims=1))

figure("ball")
clf()
fill_between(1:iters,re_mean .+ re_std, re_mean .- re_std, alpha=vis)
plot(re_mean, label="REINFORCE")
fill_between(1:iters,a2c_mean .+ a2c_std, a2c_mean .- a2c_std, alpha=vis)
plot(a2c_mean, label="A2C")
fill_between(1:iters,tnpg_mean .+ tnpg_std, tnpg_mean .- tnpg_std, alpha=vis)
plot(tnpg_mean, label="TNPG")
fill_between(1:iters,trpo_mean .+ trpo_std, trpo_mean .- trpo_std, alpha=vis)
plot(trpo_mean, label="TRPO")
fill_between(1:iters,ppo_mean .+ ppo_std, ppo_mean .- ppo_std, alpha=vis)
plot(ppo_mean, label="PPO")
fill_between(1:iters,td3_mean .+ td3_std, td3_mean .- td3_std, alpha=vis)
plot(td3_mean, label="PPO-soft")
xlabel("Iterations")
ylabel("Costs")
legend()
#ylim([0,15])

max_stars_repo_name: NikonPic/REINFORCEMENT-JULIA

max_issues_repo_name: NikonPic/REINFORCEMENT-JULIA

max_forks_repo_name: NikonPic/REINFORCEMENT-JULIA


text:
using RungeKutta
using Base.Test

## Solves the test problem dy/dt = -t^3 y^3
##
##   which has the exact solution y = 1/sqrt(C + t^4/2)
##   where C = 1/y_0^2 - t_0^4/2

summer = +
scaler = *

con = -0.4
deriv = (t,y) ->
begin
#    println (string("deriv: t = ", t, " y = ", y, " y' = ", con*y))
    con*y
end
t0 = 0.0
y0 = 1.75
exact = ((t) -> y0*exp(con*(t - t0)))

showst = ((t,y) -> string (y, "\t", (y - exact (t))))

function gen_soln1 (integrator,h,t,st)
    stn = st
    tn = t
    while tn < 5.0
        stn = integrator (tn,stn)
        tn  = tn+h
    end
    stn = integrator (tn,stn)
    tn  = tn+h
    println (showst (tn,stn))
end

function do_case1 (integrator, n)
    h = (n < 0) ? 2^(-n) : 0.5^(n)
    print (string (h, "\t\t"))
    gen_soln1 (integrator (h),h,t0,y0)
end

function solver1 (integrator,stats)
    println (stats)
    println ("# step yf err")
    for h in { x - 2 for x = 0:14 }
        do_case1 (integrator (scaler,summer,deriv), h)
    end
    println ("# All done!")
    end

solver1 (make_rkfe(), show_rkfe)
solver1 (make_rk3(), show_rk3)
solver1 (make_rk4a(), show_rk4a)
solver1 (make_rk4b(), show_rk4b)


text:
using HDF5, DataFrames, CSV, JLD
using StatsBase, Polynomials, CurveFit
using PyPlot

kep_filename = joinpath(Pkg.dir("ExoplanetsSysSim"), "data", "q1_q17_dr25_stellar.csv")
gaia_filename = joinpath(Pkg.dir("ExoplanetsSysSim"), "data", "gaiadr2_keplerdr25_crossref.csv")
mast_filename = joinpath(Pkg.dir("ExoplanetsSysSim"), "data", "KeplerMAST_TargetProperties.csv")
stellar_catalog_file_out = joinpath(Pkg.dir("ExoplanetsSysSim"), "data", "q1q17_dr25_gaia_fgk.jld")

kep_df = CSV.read(kep_filename)
gaia_df = CSV.read(gaia_filename, rows_for_type_detect=30000)
mast_df = CSV.read(mast_filename)
tmass_df = CSV.read(tmass_filename, rows_for_type_detect=30000)

dup_gaiaid = find(nonunique(DataFrame(x = gaia_df[:source_id])))
deleterows!(gaia_df, dup_gaiaid)

println("Total crossref target stars = ", length(gaia_df[:kepid]))

mag_diff = gaia_df[:phot_g_mean_mag].-gaia_df[:kepmag]
quant_arr = quantile(mag_diff, [0.067,0.933])   # 1.5-sigma cut
mag_match = find(x->quant_arr[1]<=x<=quant_arr[2], mag_diff)
gaia_df = gaia_df[mag_match,:]

gaia_col = [:kepid,:source_id,:parallax,:parallax_error,:astrometric_gof_al,:astrometric_excess_noise_sig,:phot_g_mean_mag,:bp_rp,:a_g_val,:priam_flags,:teff_val,:teff_percentile_lower,:teff_percentile_upper,:radius_val,:radius_percentile_lower,:radius_percentile_upper,:lum_val,:lum_percentile_lower,:lum_percentile_upper]
df = join(kep_df, gaia_df[:,gaia_col], on=:kepid)
df = join(df, mast_df, on=:kepid)
kep_df = nothing
gaia_df = nothing

println("Total target stars (KOIs) matching magnitude = ", length(df[:kepid]), " (", sum(df[:nkoi]),")")

df[:teff] = df[:teff_val]
df[:teff_err1] = df[:teff_percentile_upper].-df[:teff_val]
df[:teff_err2] = df[:teff_percentile_lower].-df[:teff_val]
delete!(df, :teff_val)
delete!(df, :teff_percentile_upper)
delete!(df, :teff_percentile_lower)
for x in 1:length(df[:kepid])
    if !isnan(df[x,:radius_val])
        df[x,:radius_err1] = df[x,:radius_percentile_upper]-df[x,:radius_val]
        df[x,:radius_err2] = df[x,:radius_percentile_lower]-df[x,:radius_val]
        df[x,:radius] = df[x,:radius_val]
    end
end
delete!(df, :radius_val)
delete!(df, :radius_percentile_upper)
delete!(df, :radius_percentile_lower)

not_binary = (df[:astrometric_gof_al] .<= 20) .& (df[:astrometric_excess_noise_sig] .<= 5)
astrometry_good = []
for x in 1:length(df[:kepid])
    if !(ismissing(df[x,:priam_flags]))
        pflag = string(df[x,:priam_flags])
         if (pflag[2] == '0') & (pflag[3] == '0') # WARNING: Assumes flag had first '0' removed by crossref script
             push!(astrometry_good, true)
         else
             push!(astrometry_good, false)
         end
     else
         push!(astrometry_good, false)
     end
end
astro_good = astrometry_good .& (df[:parallax_error] .< 0.1*df[:parallax])
planet_search = df[:kepmag] .<= 16.

has_mass = .! (ismissing.(df[:mass]) .| ismissing.(df[:mass_err1]) .| ismissing.(df[:mass_err2]))
has_radius = .! (ismissing.(df[:radius]) .| ismissing.(df[:radius_err1]) .| ismissing.(df[:radius_err2]))# .| isnan.(df[:radius]) .| isnan.(df[:radius_err1]) .| isnan.(df[:radius_err2]))
has_dens = .! (ismissing.(df[:dens]) .| ismissing.(df[:dens_err1]) .| ismissing.(df[:dens_err2]))
has_cdpp = .! (ismissing.(df[:rrmscdpp01p5]) .| ismissing.(df[:rrmscdpp02p0]) .| ismissing.(df[:rrmscdpp02p5]) .| ismissing.(df[:rrmscdpp03p0]) .| ismissing.(df[:rrmscdpp03p5]) .| ismissing.(df[:rrmscdpp04p5]) .| ismissing.(df[:rrmscdpp05p0]) .| ismissing.(df[:rrmscdpp06p0]) .| ismissing.(df[:rrmscdpp07p5]) .| ismissing.(df[:rrmscdpp09p0]) .| ismissing.(df[:rrmscdpp10p5]) .| ismissing.(df[:rrmscdpp12p0]) .| ismissing.(df[:rrmscdpp12p5]) .| ismissing.(df[:rrmscdpp15p0]))
has_rest = .! (ismissing.(df[:dataspan]) .| ismissing.(df[:dutycycle]))
has_limbd = .! (ismissing.(df[:limbdark_coeff1]) .| ismissing.(df[:limbdark_coeff2]) .| ismissing.(df[:limbdark_coeff3]) .| ismissing.(df[:limbdark_coeff4]))

mast_cut =.&(df[:numLCEXqtrs].>0,df[:numLCqtrs].>4)

is_usable = .&(has_radius, has_mass, has_rest, has_cdpp, has_dens, astro_good, not_binary, planet_search, mast_cut)

df = df[find(is_usable),:]
println("Total stars (KOIs) with valid parameters = ", length(df[:kepid]), " (", sum(df[:nkoi]),")")

fgk_color = (0.5 .<= df[:bp_rp] .<= 1.7)
df = df[find(fgk_color),:]
println("Total FGK colored stars (KOIs) with valid parameters = ", length(find(fgk_color)), " (", sum(df[:, :nkoi]),")")

near_FGK_MS = []
coeff = [2.5,-3.6,0.9]
for i in 1:6
    near_FGK_MS = (log10.(df[:lum_val]).< map(x->polyval(Poly(coeff),x),df[:bp_rp]) + log10(1.75))
    coeff = poly_fit(df[near_FGK_MS,:bp_rp],log10.(df[near_FGK_MS,:lum_val]),2)
end
FGK = find(near_FGK_MS)

println("Total FGK stars (KOIs) with valid parameters = ", length(FGK), " (", sum(df[FGK, :nkoi]),")")

plot_samp = rand(1:length(df[:kepid]), 5000)
plt[:scatter](df[FGK, :bp_rp], df[FGK, :lum_val], s=3, label = "FGK", color="red")
plt[:scatter](df[plot_samp, :bp_rp], df[plot_samp, :lum_val], s=3, label = "All", color="black")
plt[:ylabel](L"$L (L_\odot)$")
plt[:xlabel](L"$B_p - R_p$")
plt[:yscale]("log")
plt[:ylim](reverse(plt[:ylim]()))
plt[:legend]()
plt[:savefig]("fgk-dwarf_samp.png")

# See options at: http://exoplanetarchive.ipac.caltech.edu/docs/API_keplerstellar_columns.html
# TODO SCI DETAIL or IMPORTANT?: Read in all CDPP's, so can interpolate?
symbols_to_keep = [ :kepid, :source_id, :mass, :mass_err1, :mass_err2, :radius, :radius_err1, :radius_err2, :dens, :dens_err1, :dens_err2, :teff, :phot_g_mean_mag, :bp_rp, :lum_val, :rrmscdpp01p5, :rrmscdpp02p0, :rrmscdpp02p5, :rrmscdpp03p0, :rrmscdpp03p5, :rrmscdpp04p5, :rrmscdpp05p0, :rrmscdpp06p0, :rrmscdpp07p5, :rrmscdpp09p0, :rrmscdpp10p5, :rrmscdpp12p0, :rrmscdpp12p5, :rrmscdpp15p0, :dataspan, :dutycycle, :limbdark_coeff1, :limbdark_coeff2, :limbdark_coeff3, :limbdark_coeff4, :contam]
# delete!(df, [~(x in symbols_to_keep) for x in names(df)])    # delete columns that we won't be using anyway
df = df[FGK, symbols_to_keep]
tmp_df = DataFrame()
for col in names(df)
    tmp_df[col] = collect(skipmissing(df[col]))
end
df = tmp_df

save(stellar_catalog_file_out,"stellar_catalog", df)


text:
import PyCall: pyimport

function installmaclinux()
   # See https://stackoverflow.com/questions/12332975/installing-python-module-within-code.
   PIP_PACKAGES = ["lale","imbalanced-learn","imblearn","openml"]
   try
      pyimport("lale")
      @info "lale succesfully installed"
   catch
      try
         sys = pyimport("sys")
         subprocess = pyimport("subprocess")
         subprocess.check_call([sys.executable, "-m", "pip", "install", "--user", "--upgrade", "--force-reinstall", PIP_PACKAGES...])
         @info "lale succesfully installed"
      catch
         println("lale failed to install")
      end
   end
end

function installwindows()
   nothing
end

if Sys.iswindows()
   installwindows()
else
   installmaclinux()
end