text:
import Base.getproperty
import Base.findmin
import Base.findmax
import Base.push!

using Zygote: @adjoint, @nograd

import Zygote.literal_getproperty

# ---- #
# Vec3 #
# ---- #

# Hack to avoid nothing in the gradient due to typemax
@adjoint bigmul(x::T) where {T} = bigmul(x), Δ -> (zero(T),)

@adjoint Vec3(a, b, c) = Vec3(a, b, c), Δ -> (Δ.x, Δ.y, Δ.z)

@adjoint literal_getproperty(v::Vec3, ::Val{:x}) =
    getproperty(v, :x), Δ -> (Vec3(Δ, zero(v.y), zero(v.z)), nothing)

@adjoint literal_getproperty(v::Vec3, ::Val{:y}) =
    getproperty(v, :y), Δ -> (Vec3(zero(v.x), Δ, zero(v.z)), nothing)

@adjoint literal_getproperty(v::Vec3, ::Val{:z}) =
    getproperty(v, :z), Δ -> (Vec3(zero(v.x), zero(v.y), Δ), nothing)

@adjoint function dot(a::Vec3, b::Vec3)
    dot(a, b), Δ -> begin
        Δ = map(x -> isnothing(x) ? zero(eltype(a.x)) : x, Δ)
        t1 = Δ * b
        t2 = Δ * a
        if length(a.x) != length(t1.x)
            t1 = Vec3(sum(t1.x), sum(t1.y), sum(t1.z))
        end
        if length(b.x) != length(t2.x)
            t2 = Vec3(sum(t2.x), sum(t2.y), sum(t2.z))
        end
        return (t1, t2)
    end
end

# The purpose of this adjoint is to ensure type inference works
# in the backward pass
@adjoint function cross(a::Vec3{T}, b::Vec3{T}) where {T}
    cross(a, b), Δ -> begin
        ∇a = zero(a)
        ∇b = zero(b)
        x = (Δ.z .* b.y) .- (Δ.y .* b.z)
        y = (Δ.x .* b.z) .- (Δ.z .* b.x)
        z = (Δ.y .* b.x) .- (Δ.x .* b.y)
        if length(a.x) == 1
            ∇a.x .= sum(x)
            ∇a.y .= sum(y)
            ∇a.z .= sum(z)
        else
            ∇a.x .= x
            ∇a.y .= y
            ∇a.z .= z
        end
        x = (Δ.y .* a.z) .- (Δ.z .* a.y)
        y = (Δ.z .* a.x) .- (Δ.x .* a.z)
        z = (Δ.x .* a.y) .- (Δ.y .* a.x)
        if length(b.x) == 1
            ∇b.x .= sum(x)
            ∇b.y .= sum(y)
            ∇b.z .= sum(z)
        else
            ∇b.x .= x
            ∇b.y .= y
            ∇b.z .= z
        end
        return (∇a, ∇b)
    end
end

@adjoint place(a::Vec3, cond) = place(a, cond), Δ -> (Vec3(Δ.x[cond], Δ.y[cond], Δ.z[cond]), nothing)

@adjoint place(a::Array, cond) = place(a, cond), Δ -> (Δ[cond], nothing)

@adjoint place_idx!(a::Vec3, b::Vec3, idx) = place_idx!(a, b, idx), Δ -> (zero(Δ), Vec3(Δ[idx]...), nothing)

# ----- #
# Light #
# ----- #

# -------------- #
# - PointLight - #
# -------------- #

@adjoint PointLight(color::Vec3, intensity, pos::Vec3) =
    PointLight(color, intensity, pos), Δ -> (Δ.color, Δ.intensity, Δ.position)

@adjoint literal_getproperty(p::PointLight, ::Val{:color}) =
    getproperty(p, :color), Δ -> (PointLight(Δ, zero(p.intensity), zero(p.position)), nothing)

@adjoint literal_getproperty(p::PointLight, ::Val{:intensity}) =
    getproperty(p, :intensity), Δ -> (PointLight(zero(p.color), Δ, zero(p.position)), nothing)

@adjoint literal_getproperty(p::PointLight, ::Val{:position}) =
    getproperty(p, :position), Δ -> (PointLight(zero(p.color), zero(p.intensity), Δ), nothing)

# ---------------- #
# - DistantLight - #
# ---------------- #

@adjoint DistantLight(color::Vec3, intensity, direction::Vec3) =
    DistantLight(color, intensity, direction), Δ -> (Δ.color, Δ.intensity, Δ.direction)

@adjoint literal_getproperty(d::DistantLight, ::Val{:color}) =
    getproperty(d, :color), Δ -> (DistantLight(Δ, zero(d.intensity), zero(d.direction)), nothing)

@adjoint literal_getproperty(d::DistantLight, ::Val{:intensity}) =
    getproperty(d, :intensity), Δ -> (DistantLight(zero(d.color), Δ, zero(d.direction)), nothing)

@adjoint literal_getproperty(d::DistantLight, ::Val{:direction}) =
    getproperty(d, :direction), Δ -> (DistantLight(zero(d.color), zero(d.intensity), Δ), nothing)

# -------- #
# Material #
# -------- #

@adjoint Material(color_ambient, color_diffuse, color_specular, specular_exponent,
                  reflection, texture_ambient, texture_diffuse, texture_specular,
                  uv_coordinates) =
    Material(color_ambient, color_diffuse, color_specular, specular_exponent, reflection,
             texture_ambient, texture_diffuse, texture_specular, uv_coordinates),
    Δ -> (Δ.color_ambient, Δ.color_diffuse, Δ.color_specular, Δ.specular_exponent,
          Δ.reflection, Δ.texture_ambient, Δ.texture_diffuse, Δ.texture_specular,
          Δ.uv_coordinates)

@adjoint literal_getproperty(m::Material, ::Val{:color_ambient}) =
    getproperty(m, :color_ambient), Δ -> (Material(Δ, zero(m.color_diffuse), zero(m.color_specular),
                                                   zero(m.specular_exponent), zero(m.reflection),
                                                   isnothing(m.texture_ambient) ? nothing : zero(m.texture_ambient),
                                                   isnothing(m.texture_diffuse) ? nothing : zero(m.texture_diffuse),
                                                   isnothing(m.texture_specular) ? nothing : zero(m.texture_specular),
                                                   isnothing(m.uv_coordinates) ? nothing : zero.(m.uv_coordinates)),
                                          nothing)

@adjoint literal_getproperty(m::Material, ::Val{:color_diffuse}) =
    getproperty(m, :color_diffuse), Δ -> (Material(zero(m.color_ambient), Δ, zero(m.color_specular),
                                                   zero(m.specular_exponent), zero(m.reflection),
                                                   isnothing(m.texture_ambient) ? nothing : zero(m.texture_ambient),
                                                   isnothing(m.texture_diffuse) ? nothing : zero(m.texture_diffuse),
                                                   isnothing(m.texture_specular) ? nothing : zero(m.texture_specular),
                                                   isnothing(m.uv_coordinates) ? nothing : zero.(m.uv_coordinates)),
                                          nothing)

@adjoint literal_getproperty(m::Material, ::Val{:color_specular}) =
    getproperty(m, :color_specular), Δ -> (Material(zero(m.color_ambient), zero(m.color_diffuse), Δ,
                                                   zero(m.specular_exponent), zero(m.reflection),
                                                   isnothing(m.texture_ambient) ? nothing : zero(m.texture_ambient),
                                                   isnothing(m.texture_diffuse) ? nothing : zero(m.texture_diffuse),
                                                   isnothing(m.texture_specular) ? nothing : zero(m.texture_specular),
                                                   isnothing(m.uv_coordinates) ? nothing : zero.(m.uv_coordinates)),
                                          nothing)

@adjoint literal_getproperty(m::Material, ::Val{:specular_exponent}) =
    getproperty(m, :specular_exponent), Δ -> (Material(zero(m.color_ambient), zero(m.color_diffuse),
                                                       zero(m.color_specular), Δ, zero(m.reflection),
                                                       isnothing(m.texture_ambient) ? nothing : zero(m.texture_ambient),
                                                       isnothing(m.texture_diffuse) ? nothing : zero(m.texture_diffuse),
                                                       isnothing(m.texture_specular) ? nothing : zero(m.texture_specular),
                                                       isnothing(m.uv_coordinates) ? nothing : zero.(m.uv_coordinates)),
                                              nothing)

@adjoint literal_getproperty(m::Material, ::Val{:reflection}) =
    getproperty(m, :reflection), Δ -> (Material(zero(m.color_ambient), zero(m.color_diffuse),
                                                zero(m.color_specular), zero(m.specular_exponent), Δ,
                                                isnothing(m.texture_ambient) ? nothing : zero(m.texture_ambient),
                                                isnothing(m.texture_diffuse) ? nothing : zero(m.texture_diffuse),
                                                isnothing(m.texture_specular) ? nothing : zero(m.texture_specular),
                                                isnothing(m.uv_coordinates) ? nothing : zero.(m.uv_coordinates)),
                                       nothing)

@adjoint literal_getproperty(m::Material, ::Val{:texture_ambient}) =
    getproperty(m, :texture_ambient), Δ -> (Material(zero(m.color_ambient), zero(m.color_diffuse),
                                                     zero(m.color_specular), zero(m.specular_exponent),
                                                     zero(m.reflection), Δ,
                                                     isnothing(m.texture_diffuse) ? nothing : zero(m.texture_diffuse),
                                                     isnothing(m.texture_specular) ? nothing : zero(m.texture_specular),
                                                     isnothing(m.uv_coordinates) ? nothing : zero.(m.uv_coordinates)),
                                            nothing)

@adjoint literal_getproperty(m::Material, ::Val{:texture_diffuse}) =
    getproperty(m, :texture_diffuse), Δ -> (Material(zero(m.color_ambient), zero(m.color_diffuse),
                                                     zero(m.color_specular), zero(m.specular_exponent),
                                                     zero(m.reflection),
                                                     isnothing(m.texture_ambient) ? nothing : zero(m.texture_ambient), Δ,
                                                     isnothing(m.texture_specular) ? nothing : zero(m.texture_specular),
                                                     isnothing(m.uv_coordinates) ? nothing : zero.(m.uv_coordinates)),
                                            nothing)

@adjoint literal_getproperty(m::Material, ::Val{:texture_specular}) =
    getproperty(m, :texture_specular), Δ -> (Material(zero(m.color_ambient), zero(m.color_diffuse),
                                                      zero(m.color_specular), zero(m.specular_exponent),
                                                      zero(m.reflection),
                                                      isnothing(m.texture_ambient) ? nothing : zero(m.texture_ambient),
                                                      isnothing(m.texture_diffuse) ? nothing : zero(m.texture_diffuse), Δ,
                                                      isnothing(m.uv_coordinates) ? nothing : zero.(m.uv_coordinates)),
                                             nothing)

@adjoint literal_getproperty(m::Material, ::Val{:uv_coordinates}) =
    getproperty(m, :uv_coordinates), Δ -> (Material(zero(m.color_ambient), zero(m.color_diffuse),
                                                    zero(m.color_specular), zero(m.specular_exponent),
                                                    zero(m.reflection),
                                                    isnothing(m.texture_ambient) ? nothing : zero(m.texture_ambient),
                                                    isnothing(m.texture_diffuse) ? nothing : zero(m.texture_diffuse),
                                                    isnothing(m.texture_specular) ? nothing : zero(m.texture_specular),
                                                    Δ),
                                             nothing)
# ------- #
# Objects #
# ------- #

# ---------- #
# - Sphere - #
# ---------- #

@adjoint Sphere(center, radius, material::Material) =
    Sphere(center, radius, material), Δ -> (Δ.sphere, Δ.radius, Δ.material)

@adjoint literal_getproperty(s::Sphere, ::Val{:center}) =
    getproperty(s, :center), Δ -> (Sphere(Δ, zero(s.radius), zero(s.material)), nothing)

@adjoint literal_getproperty(s::Sphere, ::Val{:radius}) =
    getproperty(s, :radius), Δ -> (Sphere(zero(s.center), Δ, zero(s.material)), nothing)

@adjoint literal_getproperty(s::Sphere, ::Val{:material}) =
    getproperty(s, :material), Δ -> (Sphere(zero(s.center), zero(s.radius), Δ), nothing)

# ------------ #
# - Triangle - #
# ------------ #

@adjoint Triangle(v1, v2, v3, material::Material) =
    Triangle(v1, v2, v3, material), Δ -> (Δ.v1, Δ.v2, Δ.v3, Δ.material)

@adjoint literal_getproperty(t::Triangle, ::Val{:v1}) =
    getproperty(t, :v1), Δ -> (Triangle(Δ, zero(t.v2), zero(t.v3), zero(t.material)), nothing)

@adjoint literal_getproperty(t::Triangle, ::Val{:v2}) =
    getproperty(t, :v2), Δ -> (Triangle(zero(t.v1), Δ, zero(t.v3), zero(t.material)), nothing)

@adjoint literal_getproperty(t::Triangle, ::Val{:v3}) =
    getproperty(t, :v3), Δ -> (Triangle(zero(t.v1), zero(t.v2), Δ, zero(t.material)), nothing)

@adjoint literal_getproperty(t::Triangle, ::Val{:material}) =
    getproperty(t, :material), Δ -> (Triangle(zero(t.v1), zero(t.v2), zero(t.v3), Δ), nothing)

# ---------------- #
# - TriangleMesh - #
# ---------------- #

#=
@adjoint TriangleMesh(tm, mat, ftmp) =
    TriangleMesh(tm, mat, ftmp), Δ -> (Δ.triangulated_mesh, Δ.material, Δ.ftmp)

@adjoint function literal_getproperty(t::TriangleMesh, ::Val{:triangulated_mesh})
    tm = getproperty(t, :triangulated_mesh)
    z = eltype(tm[1].v1.x)
    mat = Material(PlainColor(rgb(z)), z)
    return tm, Δ -> (TriangleMesh(Δ, mat, FixedTriangleMeshParams(IdDict(), [Vec3(z)])), nothing)
end

@adjoint function literal_getproperty(t::TriangleMesh, ::Val{:material})
    mat = getproperty(t, :material)
    z = eltype(t.triangulated_mesh[1].v1.x)
    tm = [Triangle([Vec3(z)]...) for _ in 1:length(t.triangulated_mesh)]
    return mat, Δ -> (TriangleMesh(tm, Δ, FixedTriangleMeshParams(IdDict(), [Vec3(z)])), nothing)
end

@adjoint function literal_getproperty(t::TriangleMesh, ::Val{:ftmp})
    z = eltype(t.triangulated_mesh[1].v1.x)
    mat = Material(PlainColor(rgb(z)), z)
    tm = [Triangle([Vec3(z)]...) for _ in 1:length(t.triangulated_mesh)]
    return getproperty(t, :ftmp), Δ -> (TriangleMesh(tm, mat, Δ), nothing)
end

@adjoint FixedTriangleMeshParams(isect, n) =
    FixedTriangleMeshParams(isect, n), Δ -> (Δ.isect, Δ.n)

# The gradients for this params are never used so fill them with anything as long
# as they are consistent with the types
@adjoint literal_getproperty(ftmp::FixedTriangleMeshParams, ::Val{f}) where {f} =
    getproperty(ftmp, f), Δ -> (FixedTriangleMeshParams(IdDict(), ftmp.normals[1:1]))
=#

# ------ #
# Camera #
# ------ #

@adjoint Camera(lf, la, vfov, focus, fp) =
    Camera(lf, la, vfov, focus, fp), Δ -> (Δ.lookfrom, Δ.lookat, Δ.vfov,
                                           Δ.focus, Δ.fixedparams)

@adjoint literal_getproperty(c::Camera, ::Val{:lookfrom}) =
    getproperty(c, :lookfrom), Δ -> (Camera(Δ, zero(c.lookat), zero(c.vfov), zero(c.focus),
                                            zero(c.fixedparams)), nothing)

@adjoint literal_getproperty(c::Camera, ::Val{:lookat}) =
    getproperty(c, :lookat), Δ -> (Camera(zero(c.lookfrom), Δ, zero(c.vfov), zero(c.focus),
                                          zero(c.fixedparams)), nothing)

@adjoint literal_getproperty(c::Camera, ::Val{:vfov}) =
    getproperty(c, :vfov), Δ -> (Camera(zero(c.lookfrom), zero(c.lookat), Δ, zero(c.focus),
                                        zero(c.fixedparams)), nothing)

@adjoint literal_getproperty(c::Camera, ::Val{:focus}) =
    getproperty(c, :focus), Δ -> (Camera(zero(c.lookfrom), zero(c.lookat), zero(c.vfov), Δ,
                                         zero(c.fixedparams)), nothing)

@adjoint literal_getproperty(c::Camera, ::Val{:fixedparams}) =
    getproperty(c, :fixedparams), Δ -> (Camera(zero(c.lookfrom), zero(c.lookat), zero(c.vfov),
                                               zero(c.focus), Δ), nothing)

@adjoint FixedCameraParams(vup, w, h) =
    FixedCameraParams(vup, w, h), Δ -> (Δ.vup, Δ.width, Δ.height)

@adjoint literal_getproperty(fcp::FixedCameraParams, ::Val{f}) where {f} =
    getproperty(fcp, f), Δ -> (zero(fcp), nothing)

# ------- #
# ImUtils #
# ------- #

@adjoint function zeroonenorm(x)
    mini, indmin = findmin(x)
    maxi, indmax = findmax(x)
    res = (x .- mini) ./ (maxi - mini)
    function ∇zeroonenorm(Δ)
        ∇x = similar(x)
        fill!(∇x, 1 / (maxi - mini))
        res1 = (x .- maxi) ./ (maxi - mini)^2
        ∇x[indmin] = sum(res1) - minimum(res1)
        res2 = - res ./ (maxi - mini)
        ∇x[indmax] = sum(res2) - minimum(res2)
        return (∇x .* Δ, )
    end
    return res, ∇zeroonenorm
end

# ----------------- #
# General Functions #
# ----------------- #

for func in (:findmin, :findmax)
    @eval begin
        @adjoint function $(func)(xs::AbstractArray; dims = :)
            y = $(func)(xs, dims = dims)
            function dfunc(Δ)
                res = zero(xs)
                res[y[2]] .= Δ[1]
                return (res, nothing)
            end
            return y, dfunc
        end
    end
end

@adjoint reducehcat(x) = reduce(hcat, x), Δ -> ([Δ[:, i] for i in 1:length(x)], )

@adjoint push!(arr, val) = push!(arr, val), Δ -> (Δ[1:end-1], Δ[end])

@nograd fill

@nograd function update_index!(arr, i, j, val)
    arr[i, j] = val
end


text:
@testset "pipe" begin

# test pipe

JULIA = joinpath(get(ENV, "JULIA_HOME", ""), "julia")

READLINE_CMD = "readline(stdin; keep=true)"

@static if Sys.iswindows()
    PIPETEST = joinpath(dirname(@__FILE__), "pipetest.cmd")
    TESTJL = "\"while !eof(stdin); print(\\\"Julia\\\", $READLINE_CMD); end\""
    TESTJLENV = "\"while !eof(stdin); print(ENV[\\\"HEADER\\\"], $READLINE_CMD); end\""
else
    PIPETEST = joinpath(dirname(@__FILE__), "pipetest.sh")
    TESTJL = """while !eof(stdin); print("Julia", $READLINE_CMD); end"""
    TESTJLENV = """while !eof(stdin); print(ENV["HEADER"], $READLINE_CMD); end"""
end

sc = SparkContext(master="local")

nums1 = parallelize(sc, 1:30)

nums11 = collect(pipe(nums1, PIPETEST))
@test length(nums11) > 30

nums11 = collect(pipe(nums1, [JULIA, "-e", TESTJL]))
@test length(nums11) == 30
for l in nums11
    @test startswith(l, "Julia")
end

nums11 = collect(pipe(nums1, [JULIA, "-e", TESTJLENV], Dict{String,String}("HEADER"=>"Julia")))
@test length(nums11) == 30
for l in nums11
    @test startswith(l, "Julia")
end

pnums1 = cartesian(nums1, nums1)

pnums11 = collect(pipe(pnums1, PIPETEST))
@test length(pnums11) > 900

pnums11 = collect(pipe(pnums1, [JULIA, "-e", TESTJL]))
@test length(pnums11) == 900
for l in pnums11
    @test startswith(l, "Julia")
end

pnums11 = collect(pipe(pnums1, [JULIA, "-e", TESTJLENV], Dict{String,String}("HEADER"=>"Julia")))
@test length(pnums11) == 900
for l in pnums11
    @test startswith(l, "Julia")
end

close(sc)

end

ext: jl


text:
# Some modified wrappers for Lapack
import LinearAlgebra: BlasFloat, BlasInt, LAPACKException,
    DimensionMismatch, SingularException, PosDefException, chkstride1, checksquare
import LinearAlgebra.BLAS: @blasfunc, libblas, BlasReal, BlasComplex
import LinearAlgebra.LAPACK: liblapack, chklapackerror

@static if isdefined(Base, :require_one_based_indexing)
    import Base: require_one_based_indexing
else
    require_one_based_indexing(A...) = !Base.has_offset_axes(A...) || throw(ArgumentError("offset arrays are not supported but got an array with index other than 1"))
end

struct RowIterator{A<:AbstractMatrix,R<:IndexRange}
    a::A
    r::R
end
rows(a::AbstractMatrix, r::IndexRange = axes(a,1)) = RowIterator(a, r)

function Base.iterate(iter::RowIterator)
    next = iterate(iter.r)
    if next === nothing
        return nothing
    else
        i, s = next
        return view(iter.a, i, :), s
    end
end
function Base.iterate(iter::RowIterator, s)
    next = iterate(iter.r, s)
    if next === nothing
        return nothing
    else
        i, s = next
        return view(iter.a, i, :), s
    end
end

Base.IteratorSize(::Type{<:RowIterator}) = Base.HasLength()
Base.IteratorEltype(::Type{<:RowIterator}) = Base.HasEltype()
Base.length(iter::RowIterator) = length(iter.r)
Base.eltype(iter::RowIterator{A}) where {T,A<:DenseArray{T}} =
    SubArray{T,1,A,Tuple{Int,Base.Slice{Base.OneTo{Int}}},true}

struct ColumnIterator{A<:AbstractMatrix,R<:IndexRange}
    a::A
    r::R
end
cols(a::AbstractMatrix, r::IndexRange = axes(a,2)) = ColumnIterator(a, r)

function Base.iterate(iter::ColumnIterator)
    next = iterate(iter.r)
    if next === nothing
        return nothing
    else
        i, s = next
        return view(iter.a, :, i), s
    end
end
function Base.iterate(iter::ColumnIterator, s)
    next = iterate(iter.r, s)
    if next === nothing
        return nothing
    else
        i, s = next
        return view(iter.a, :, i), s
    end
end

Base.IteratorSize(::Type{<:ColumnIterator}) = Base.HasLength()
Base.IteratorEltype(::Type{<:ColumnIterator}) = Base.HasEltype()
Base.length(iter::ColumnIterator) = length(iter.r)
Base.eltype(iter::ColumnIterator{A}) where {T,A<:DenseArray{T}} =
    SubArray{T,1,A,Tuple{Base.Slice{Base.OneTo{Int}},Int},true}

# # QR decomposition
# function qr!(A::StridedMatrix{<:BlasFloat})
#     m, n = size(A)
#     A, T = LAPACK.geqrt!(A, min(minimum(size(A)), 36))
#     Q = LAPACK.gemqrt!('L', 'N', A, T, eye(eltype(A), m, min(m,n)))
#     R = triu!(A[1:min(m,n), :])
#     return Q, R
# end

# Triangular division: for some reason this is faster than LAPACK's trsv
function ldiv!(A::UpperTriangular, y::AbstractVector, r::UnitRange{Int} = 1:length(y))
    R = A.data
    @inbounds for j in reverse(r)
        R[j,j] == zero(R[j,j]) && throw(SingularException(j))
        yj = (y[j] = R[j,j] \ y[j])
        @simd for i in first(r):j-1
            y[i] -= R[i,j] * yj
        end
    end
    return y
end

# Eigenvalue decomposition of SymTridiagonal matrix
tridiageigh!(A::SymTridiagonal{T}, Z::StridedMatrix{T} = one(A)) where {T<:BlasFloat} =
    stegr!(A.dv, A.ev, Z) # redefined

# Generalized eigenvalue decomposition of symmetric / Hermitian problem
geneigh!(A::StridedMatrix{T}, B::StridedMatrix{T}) where {T<:BlasFloat} =
    LAPACK.sygvd!(1, 'V', 'U', A, B)

# Singular value decomposition of a Bidiagonal matrix
function bidiagsvd!(B::Bidiagonal{T}, U::StridedMatrix{T} = one(B),
                    VT::StridedMatrix{T} = one(B)) where {T<:BlasReal}
    s, Vt, U, = LAPACK.bdsqr!(B.uplo, B.dv, B.ev, VT, U, similar(U, (size(B,1), 0)))
    return U, s, Vt
end

function reversecols!(U::AbstractMatrix)
    n = size(U, 2)
    @inbounds for j = 1:div(n, 2)
        @simd for i = 1:size(U, 1)
            U[i,j], U[i, n+1-j] = U[i,n+1-j], U[i,j]
        end
    end
    return U
end
function reverserows!(V::AbstractVecOrMat)
    m = size(V, 1)
    @inbounds for j = 1:size(V,2)
        @simd for i = 1:div(m, 2)
            V[i, j], V[m+1-i, j] = V[m+1-i,j], V[i,j]
        end
    end
    return V
end

# Schur factorization of a Hessenberg matrix
hschur!(H::StridedMatrix{T}, Z::StridedMatrix{T} = one(H)) where {T<:BlasFloat} =
    hseqr!(H, Z)

schur2eigvals(T::StridedMatrix{<:BlasFloat}) = schur2eigvals(T, 1:size(T,1))

function schur2eigvals(T::StridedMatrix{<:BlasComplex}, which::AbstractVector{Int})
    n = checksquare(T)
    which2 = unique(which)
    length(which2) == length(which) ||
        throw(ArgumentError("which should contain unique values"))
    return [T[i,i] for i in which2]
end

function schur2eigvals(T::StridedMatrix{<:BlasReal}, which::AbstractVector{Int})
    n = checksquare(T)
    which2 = unique(which)
    length(which2) == length(which) ||
        throw(ArgumentError("which should contain unique values"))
    D = zeros(Complex{eltype(T)}, length(which2))
    for k = 1:length(which)
        i = which[k]
        if i < n && !iszero(T[i+1,i])
            halftr = (T[i,i]+T[i+1,i+1])/2
            diff = (T[i,i]-T[i+1,i+1])/2
            d = diff*diff + T[i,i+1]*T[i+1,i]  # = hafltr*halftr - det
            D[i] = halftr + im*sqrt(-d)
        elseif i > 1 && !iszero(T[i,i-1])
            halftr = (T[i,i]+T[i-1,i-1])/2
            diff = -(T[i,i]-T[i-1,i-1])/2
            d = diff*diff + T[i,i-1]*T[i-1,i]  # = hafltr*halftr - det
            D[i] = halftr - im*sqrt(-d)
        else
            D[i] = T[i,i]
        end
    end
    return D
end

function _normalizevecs!(V)
    @inbounds for k = 1:size(V,2)
        normalize!(view(V, :, k))
    end
    return V
end
function schur2eigvecs(T::StridedMatrix{<:BlasComplex})
    n = checksquare(T)
    VR = similar(T, n, n)
    VL = similar(T, n, 0)
    select = Vector{BlasInt}(undef, 0)
    trevc!('R','A', select, T, VL, VR)
    return _normalizevecs!(VR)
end
function schur2eigvecs(T::StridedMatrix{<:BlasComplex}, which::AbstractVector{Int})
    n = checksquare(T)
    which2 = unique(which)
    length(which2) == length(which) ||
        throw(ArgumentError("which should contain unique values"))
    m = BlasInt(length(which2))
    VR = similar(T, n, m)
    VL = similar(T, n, 0)

    select = zeros(BlasInt, n)
    for k = 1:length(which2)
        i = which2[k]
        select[i] = one(BlasInt)
        trevc!('R','S', select, T, VL, view(VR,:,k:k))
        select[i] = zero(BlasInt)
    end
    return _normalizevecs!(VR)
end
function schur2eigvecs(T::StridedMatrix{<:BlasReal})
    n = checksquare(T)
    VR = similar(T, Complex{eltype(T)}, n, n)
    VR′ = similar(T, n, n)
    VL′ = similar(T, n, 0)
    select = Vector{BlasInt}(undef, 0)
    trevc!('R','A', select, T, VL′, VR′)
    i = 1
    while i <= n
        if i == n || iszero(T[i+1,i])
            @inbounds @simd for k = 1:n
                VR[k,i]= VR′[k,i]
            end
            i += 1
        else
            @inbounds @simd for k = 1:n
                VR[k,i]= VR′[k,i] + im*VR′[k,i+1]
                VR[k,i+1]= VR′[k,i] - im*VR′[k,i+1]
            end
            i += 2
        end
    end
    return _normalizevecs!(VR)
end
function schur2eigvecs(T::StridedMatrix{<:BlasReal}, which::AbstractVector{Int})
    n = checksquare(T)
    which2 = unique(which)
    length(which2) == length(which) ||
        throw(ArgumentError("which should contain unique values"))
    m = length(which2)
    VR = similar(T, Complex{eltype(T)}, n, m)
    VR′ = similar(T, n, 2)
    VL′ = similar(T, n, 0)

    select = zeros(BlasInt, n)
    i = 1
    while i <= n
        if i == n || iszero(T[i+1,i])
            j = findfirst(isequal(i), which2)
            if j !== nothing
                select[i] = one(BlasInt)
                trevc!('R','S', select, T, VL′, VR′)
                @inbounds @simd for k = 1:n
                    VR[k,j]= VR′[k,1]
                end
                select[i] = zero(BlasInt)
            end
            i += 1
        else
            j1 = findfirst(isequal(i), which2)
            j2 = findfirst(isequal(i+1), which2)
            if j1 !== nothing || j2 !== nothing
                select[i] = one(BlasInt)
                select[i+1] = one(BlasInt)
                trevc!('R','S', select, T, VL′, VR′)
                @inbounds @simd for k = 1:n
                    if j1 !== nothing
                        VR[k,j1]= VR′[k,1] + im*VR′[k,2]
                    end
                    if j2 !== nothing
                        VR[k,j2]= VR′[k,1] - im*VR′[k,2]
                    end
                end
                select[i] = zero(BlasInt)
                select[i+1] = zero(BlasInt)
            end
            i += 2
        end
    end
    return _normalizevecs!(VR)
end

function permuteeig!(D::StridedVector{S}, V::StridedMatrix{S},
                        perm::AbstractVector{Int}) where {S}
    n = checksquare(V)
    p = collect(perm) # makes copy cause will be overwritten
    isperm(p) && length(p) == n ||
        throw(ArgumentError("not a valid permutation of length $n"))
    i = 1
    @inbounds while true
        if p[i] == i
            i = 1
            while i <= n && p[i] == i
                i += 1
            end
            i > n && break
        else
            iprev = findfirst(isequal(i), p)
            inext = p[i]
            p[iprev] = inext
            p[i] = i

            D[i], D[inext] = D[inext], D[i]
            for j = 1:n
                V[j,i], V[j,inext] = V[j,inext], V[j,i]
            end
            i = inext
        end
    end
    return D, V
end

permuteschur!(T::StridedMatrix{<:BlasFloat}, p::AbstractVector{Int}) =
    permuteschur!(T, one(T), p)
function permuteschur!(T::StridedMatrix{S}, Q::StridedMatrix{S},
        perm::AbstractVector{Int}) where {S<:BlasComplex}
    n = checksquare(T)
    p = collect(perm) # makes copy cause will be overwritten
    isperm(p) && length(p) == n ||
        throw(ArgumentError("not a valid permutation of length $n"))
    @inbounds for i = 1:n
        ifirst::BlasInt = p[i]
        ilast::BlasInt = i
        T, Q = LAPACK.trexc!(ifirst, ilast, T, Q)
        for k = (i+1):n
            if p[k] < p[i]
                p[k] += 1
            end
        end
    end
    return T, Q
end

function permuteschur!(T::StridedMatrix{S}, Q::StridedMatrix{S},
        perm::AbstractVector{Int}) where {S<:BlasReal}
    n = checksquare(T)
    p = collect(perm) # makes copy cause will be overwritten
    isperm(p) && length(p) == n ||
        throw(ArgumentError("not a valid permutation of length $n"))
    i = 1
    @inbounds while i <= n
        ifirst::BlasInt = p[i]
        ilast::BlasInt = i
        if ifirst == n || iszero(T[ifirst+1,ifirst])
            T, Q = LAPACK.trexc!(ifirst, ilast, T, Q)
            @inbounds for k = (i+1):n
                if p[k] < p[i]
                    p[k] += 1
                end
            end
            i += 1
        else
            p[i+1] == ifirst+1 ||
                error("cannot split 2x2 blocks when permuting schur decomposition")
            T, Q = LAPACK.trexc!(ifirst, ilast, T, Q)
            @inbounds for k = (i+2):n
                if p[k] < p[i]
                    p[k] += 2
                end
            end
            i += 2
        end
    end
    return T,Q
end

# redefine LAPACK interface to tridiagonal eigenvalue problem
for (stegr, elty) in
    ((:dstegr_,:Float64),
     (:sstegr_,:Float32))
    @eval begin
        function stegr!(dv::AbstractVector{$elty}, ev::AbstractVector{$elty}, Z::AbstractMatrix{$elty})
            require_one_based_indexing(dv, ev, Z)
            chkstride1(dv, ev, Z)
            n = length(dv)
            if length(ev) == n - 1
                eev = [ev; zero($elty)]
            elseif length(ev) == n
                eev = ev
            else
                throw(DimensionMismatch("ev has length $(length(ev)) but needs one less than dv's length, $n)"))
            end
            checksquare(Z) == n || throw(DimensionMismatch())
            ldz = max(1, stride(Z,2))
            jobz = 'V'
            range = 'A'
            abstol = Vector{$elty}(undef, 1)
            il = 1
            iu = n
            vl = zero($elty)
            vu = zero($elty)
            m = Ref{BlasInt}()
            w = similar(dv, $elty, n)
            isuppz = similar(dv, BlasInt, 2*size(Z, 2))
            work = Vector{$elty}(undef, 1)
            lwork = BlasInt(-1)
            iwork = Vector{BlasInt}(undef, 1)
            liwork = BlasInt(-1)
            info = Ref{BlasInt}()
            for i = 1:2  # first call returns lwork as work[1] and liwork as iwork[1]
                ccall((@blasfunc($stegr), liblapack), Cvoid,
                    (Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ptr{$elty},
                    Ptr{$elty}, Ref{$elty}, Ref{$elty}, Ref{BlasInt},
                    Ref{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty},
                    Ptr{$elty}, Ref{BlasInt}, Ptr{BlasInt}, Ptr{$elty},
                    Ref{BlasInt}, Ptr{BlasInt}, Ref{BlasInt}, Ptr{BlasInt},
                    Clong, Clong),
                    jobz, range, n, dv,
                    eev, vl, vu, il,
                    iu, abstol, m, w,
                    Z, ldz, isuppz, work,
                    lwork, iwork, liwork, info,
                    1, 1)
                chklapackerror(info[])
                if i == 1
                    lwork = BlasInt(work[1])
                    resize!(work, lwork)
                    liwork = iwork[1]
                    resize!(iwork, liwork)
                end
            end
            w, Z
        end
    end
end

# redefine LAPACK interface to schur
for (hseqr, trevc, trsen, elty) in
    ((:dhseqr_, :dtrevc_, :dtrsen_, :Float64),
     (:shseqr_, :strevc_, :stgsen_, :Float32))
    @eval begin
        function hseqr!(H::StridedMatrix{$elty}, Z::StridedMatrix{$elty} = one(H))
            require_one_based_indexing(H, Z)
            chkstride1(H, Z)
            n = checksquare(H)
            checksquare(Z) == n || throw(DimensionMismatch())
            job = 'S'
            compz = 'V'
            ilo = 1
            ihi = n
            ldh = stride(H, 2)
            ldz = stride(Z, 2)
            wr = similar(H, $elty, n)
            wi = similar(H, $elty, n)
            work = Vector{$elty}(undef, 1)
            lwork = BlasInt(-1)
            info = Ref{BlasInt}()
            for i = 1:2  # first call returns lwork as work[1]
                ccall((@blasfunc($hseqr), liblapack), Cvoid,
                    (Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ref{BlasInt},
                        Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ptr{$elty},
                        Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt},
                        Ptr{BlasInt}, Clong, Clong),
                        job, compz, n, ilo, ihi,
                        H, ldh, wr, wi,
                        Z, ldz, work, lwork,
                        info, 1, 1)
                chklapackerror(info[])
                if i == 1
                    lwork = BlasInt(real(work[1]))
                    resize!(work, lwork)
                end
            end
            H, Z, complex.(wr, wi)
        end
        function trevc!(side::Char, howmny::Char, select::StridedVector{BlasInt},
                        T::AbstractMatrix{$elty}, VL::AbstractMatrix{$elty},
                        VR::AbstractMatrix{$elty})
            # Extract
            if side ∉ ['L','R','B']
                throw(ArgumentError("side argument must be 'L' (left eigenvectors), 'R' (right eigenvectors), or 'B' (both), got $side"))
            end
            n = checksquare(T)
            mm = side == 'L' ? size(VL,2) : (side == 'R' ? size(VR,2) : min(size(VL,2), size(VR,2)))
            ldt, ldvl, ldvr = stride(T, 2), stride(VL, 2), stride(VR, 2)

            # Check
            chkstride1(T, select, VL, VR)

            # Allocate
            m = Ref{BlasInt}()
            work = Vector{$elty}(undef, 3n)
            info = Ref{BlasInt}()

            ccall((@blasfunc($trevc), liblapack), Cvoid,
                (Ref{UInt8}, Ref{UInt8}, Ptr{BlasInt}, Ref{BlasInt},
                 Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt},
                 Ptr{$elty}, Ref{BlasInt}, Ref{BlasInt}, Ptr{BlasInt},
                 Ptr{$elty}, Ptr{BlasInt}, Clong, Clong),
                side, howmny, select, n,
                T, ldt, VL, ldvl,
                VR, ldvr, mm, m,
                work, info, 1, 1)
            chklapackerror(info[])

            return VL, VR, m
        end
        function trsen!(job::Char, compq::Char, select::AbstractMatrix{BlasInt},
                        T::AbstractMatrix{$elty}, Q::AbstractMatrix{$elty})
            chkstride1(T, Q, select)
            n = checksquare(T)
            ldt = max(1, stride(T, 2))
            ldq = max(1, stride(Q, 2))
            wr = similar(T, $elty, n)
            wi = similar(T, $elty, n)
            m = sum(select)
            work = Vector{$elty}(undef, 1)
            lwork = BlasInt(-1)
            iwork = Vector{BlasInt}(undef, 1)
            liwork = BlasInt(-1)
            info = Ref{BlasInt}()
            select = convert(Array{BlasInt}, select)
            s = Ref{$elty}(zero($elty))
            sep = Ref{$elty}(zero($elty))
            for i = 1:2  # first call returns lwork as work[1] and liwork as iwork[1]
                ccall((@blasfunc($trsen), liblapack), Cvoid,
                    (Ref{UInt8}, Ref{UInt8}, Ptr{BlasInt}, Ref{BlasInt},
                    Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt},
                    Ptr{$elty}, Ptr{$elty}, Ref{BlasInt}, Ref{$elty}, Ref{$elty},
                    Ptr{$elty}, Ref{BlasInt}, Ptr{BlasInt}, Ref{BlasInt},
                    Ptr{BlasInt}, Clong, Clong),
                    job, compq, select, n,
                    T, ldt, Q, ldq,
                    wr, wi, m, s, sep,
                    work, lwork, iwork, liwork,
                    info, 1, 1)
                chklapackerror(info[])
                if i == 1 # only estimated optimal lwork, liwork
                    lwork  = BlasInt(real(work[1]))
                    resize!(work, lwork)
                    liwork = BlasInt(real(iwork[1]))
                    resize!(iwork, liwork)
                end
            end
            T, Q, complex.(wr, wi), s[], sep[]
        end
    end
end

 for (hseqr, trevc, trsen, elty, relty) in
    ((:zhseqr_, :ztrevc_, :ztrsen_, :ComplexF64, :Float64),
     (:chseqr_, :ctrevc_, :ctrsen_, :ComplexF32, :Float32))
    @eval begin
        function hseqr!(H::AbstractMatrix{$elty}, Z::AbstractMatrix{$elty} = one(H))
            require_one_based_indexing(H, Z)
            chkstride1(H, Z)
            n     = checksquare(H)
            checksquare(Z) == n || throw(DimensionMismatch())
            job = 'S'
            compz = 'V'
            ilo = 1
            ihi = n
            ldh = stride(H, 2)
            ldz = stride(Z, 2)
            w    = similar(H, $elty, n)
            work  = Vector{$elty}(undef, 1)
            lwork = BlasInt(-1)
            info  = Ref{BlasInt}()
            for i = 1:2  # first call returns lwork as work[1]
                ccall((@blasfunc($hseqr), liblapack), Cvoid,
                    (Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ref{BlasInt},
                        Ptr{$elty}, Ref{BlasInt}, Ptr{$elty},
                        Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt},
                        Ptr{BlasInt}, Clong, Clong),
                        job, compz, n, ilo, ihi,
                        H, ldh, w,
                        Z, ldz, work, lwork,
                        info, 1, 1)
                chklapackerror(info[])
                if i == 1
                    lwork = BlasInt(real(work[1]))
                    resize!(work, lwork)
                end
            end
            H, Z, w
        end
        function trevc!(side::Char, howmny::Char, select::AbstractVector{BlasInt},
                        T::AbstractMatrix{$elty}, VL::AbstractMatrix{$elty} = similar(T),
                        VR::AbstractMatrix{$elty} = similar(T))

            # Check
            require_one_based_indexing(select, T, VL, VR)
            chkstride1(T, select, VL, VR)

            # Extract
            if side ∉ ['L','R','B']
                throw(ArgumentError("side argument must be 'L' (left eigenvectors), 'R' (right eigenvectors), or 'B' (both), got $side"))
            end
            n = checksquare(T)
            mm = side == 'L' ? size(VL,2) : (side == 'R' ? size(VR,2) : min(size(VL,2), size(VR,2)))
            ldt, ldvl, ldvr = stride(T, 2), stride(VL, 2), stride(VR, 2)

            # Allocate
            m = Ref{BlasInt}()
            work = Vector{$elty}(undef, 2n)
            rwork = Vector{$relty}(undef, n)
            info = Ref{BlasInt}()
            ccall((@blasfunc($trevc), liblapack), Cvoid,
                (Ref{UInt8}, Ref{UInt8}, Ptr{BlasInt}, Ref{BlasInt},
                 Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt},
                 Ptr{$elty}, Ref{BlasInt}, Ref{BlasInt}, Ptr{BlasInt},
                 Ptr{$elty}, Ptr{$relty}, Ptr{BlasInt}, Clong, Clong),
                side, howmny, select, n,
                T, ldt, VL, ldvl,
                VR, ldvr, mm, m,
                work, rwork, info, 1, 1)
            chklapackerror(info[])

            return VL, VR, m
        end
        function trsen!(job::Char, compq::Char, select::StridedVector{BlasInt},
                        T::StridedMatrix{$elty}, Q::StridedMatrix{$elty})
            chkstride1(select, T, Q)
            n = checksquare(T)
            ldt = max(1, stride(T, 2))
            ldq = max(1, stride(Q, 2))
            w = similar(T, $elty, n)
            m = sum(select)
            work = Vector{$elty}(undef, 1)
            lwork = BlasInt(-1)
            info = Ref{BlasInt}()
            select = convert(Array{BlasInt}, select)
            s = Ref{$relty}(zero($relty))
            sep = Ref{$relty}(zero($relty))
            for i = 1:2  # first call returns lwork as work[1]
                ccall((@blasfunc($trsen), liblapack), Nothing,
                    (Ref{UInt8}, Ref{UInt8}, Ptr{BlasInt}, Ref{BlasInt},
                    Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt},
                    Ptr{$elty}, Ref{BlasInt}, Ref{$relty}, Ref{$relty},
                    Ptr{$elty}, Ref{BlasInt},
                    Ptr{BlasInt}),
                    job, compq, select, n,
                    T, ldt, Q, ldq,
                    w, m, s, sep,
                    work, lwork,
                    info)
                chklapackerror(info[])
                if i == 1 # only estimated optimal lwork, liwork
                    lwork  = BlasInt(real(work[1]))
                    resize!(work, lwork)
                end
            end
            T, Q, w, s[], sep[]
        end
    end
end


text:
@doc raw"""
θ-pinch equilibrium in (x,y,z) coordinates with covariant components of
the vector potential given by
```math
A (x,y,z) = \frac{B_0}{2} \big( - y , \, x , \, 0 \big)^T
```
resulting in the magnetic field with covariant components
```math
B (x,y,z) = \big( 0 , \, 0 , \, B_0 \big)^T .
```

Parameters:
    B₀: B-field at magnetic axis
"""
module ThetaPinch

    using RecipesBase
    using LaTeXStrings

    import ..ElectromagneticFields
    import ..ElectromagneticFields: CartesianEquilibrium, code
    import ..AnalyticCartesianField: X, Y, Z

    export ThetaPinchEquilibrium

    const DEFAULT_B₀ = 1.0

    struct ThetaPinchEquilibrium{T <: Number} <: CartesianEquilibrium
        name::String
        B₀::T

        function ThetaPinchEquilibrium{T}(B₀::T) where T <: Number
            new("ThetaPinchEquilibrium", B₀)
        end
    end

    ThetaPinchEquilibrium(B₀::T=DEFAULT_B₀) where T <: Number = ThetaPinchEquilibrium{T}(B₀)

    function init(B₀=DEFAULT_B₀)
        ThetaPinchEquilibrium(B₀)
    end

    macro code(B₀=DEFAULT_B₀)
        code(init(B₀); escape=true)
    end

    function Base.show(io::IO, equ::ThetaPinchEquilibrium)
        print(io, "θ-Pinch Equilibrium in (x,y,z) Coordinates with\n")
        print(io, "  B₀ = ", equ.B₀)
    end

    r²(x::AbstractVector, equ::ThetaPinchEquilibrium) = X(x,equ)^2 + Y(x,equ)^2
    r(x::AbstractVector, equ::ThetaPinchEquilibrium) = sqrt(r²(x,equ))
    R(x::AbstractVector, equ::ThetaPinchEquilibrium) = r(x,equ)
    θ(x::AbstractVector, equ::ThetaPinchEquilibrium) = atan(Y(x,equ), X(x,equ))
    ϕ(x::AbstractVector, equ::ThetaPinchEquilibrium) = θ(x,equ)

    ElectromagneticFields.A₁(x::AbstractVector, equ::ThetaPinchEquilibrium) = - equ.B₀ * Y(x,equ) / 2
    ElectromagneticFields.A₂(x::AbstractVector, equ::ThetaPinchEquilibrium) = + equ.B₀ * X(x,equ) / 2
    ElectromagneticFields.A₃(x::AbstractVector, equ::ThetaPinchEquilibrium) = zero(eltype(x))

    ElectromagneticFields.get_functions(::ThetaPinchEquilibrium) = (X=X, Y=Y, Z=Z, R=R, r=r, θ=θ, ϕ=ϕ, r²=r²)

    @recipe function f(equ::ThetaPinchEquilibrium;
                       nx = 100, ny = 100, levels = 20, size = (800,400),
                       xlims = (-1., +1.),
                       ylims = (-1., +1.))

        xgrid = LinRange(xlims[1], xlims[2], nx)
        ygrid = LinRange(ylims[1], ylims[2], ny)
        pot1  = [ElectromagneticFields.A₁([xgrid[i], ygrid[j], 0.0], equ) for i in eachindex(xgrid), j in eachindex(ygrid)]
        pot2  = [ElectromagneticFields.A₂([xgrid[i], ygrid[j], 0.0], equ) for i in eachindex(xgrid), j in eachindex(ygrid)]

        seriestype   := :contour
        aspect_ratio := :equal
        layout := (1,2)
        size   := size
        xlims  := xlims
        ylims  := ylims
        levels := levels
        legend := :none

        @series begin
            subplot := 1
            title  := L"A_x (x,y)"
            xguide := L"x"
            yguide := L"y"
            (xgrid, ygrid, pot1)
        end

        @series begin
            subplot := 2
            title  := L"A_y (x,y)"
            xguide := L"x"
            yguide := L"y"
            (xgrid, ygrid, pot2)
        end
    end

end


text:
include("prime.jl")

const N=10000000
primes=prime_table(N);

function is_truncatable_prime(n::Int)
	nn=n

	while n>=1
		if ! primes[n]
			return false
		end
		n=floor(Int,n/=10)
	end

	i=10
	while i<nn
		if !primes[nn%i]
			return false
		end

		i*=10
	end

	return true
end

function compute()
    ans=0
	index=11
	cnt=0
	while true
		if is_truncatable_prime(index)
			ans+=index
			cnt+=1
		end
		if cnt==11
			return ans
		end
		index+=1
	end

end

println(compute())


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

using Documenter, ClipData

makedocs(
	sitename = "ClipData.jl",
	pages = Any[
		"Introduction" => "index.md",
		"API" => "api/api.md"],
	format = Documenter.HTML(
		canonical = "https://pdeffebach.github.io/ClipData.jl/stable/"
	))

deploydocs(
    repo = "github.com/pdeffebach/ClipData.jl.git",
    target = "build",
    deps = nothing,
    make = nothing)

max_forks_repo_name: pdeffebach/DataClip


text:
using Printf
using LinearAlgebra
BLAS.set_num_threads(1)  # 1 thread only

import Tulip
using DoubleFloats

function tulip_d64(finst, ϵ)
    tlp = Tulip.Model{Double64}()
    Tulip.load_problem!(tlp, finst)

    # Set parameters
    tlp.params.OutputLevel = 1
    tlp.params.Presolve = 1
    tlp.params.BarrierIterationsLimit = 500

    # Set tolerances
    tlp.params.BarrierTolerancePFeas = ϵ
    tlp.params.BarrierToleranceDFeas = ϵ
    tlp.params.BarrierToleranceRGap = ϵ
    tlp.params.BarrierToleranceIFeas = ϵ / 100

    # Solve model
    t = @elapsed Tulip.optimize!(tlp)

    println("Solve time: $t")
    println("Barrier iterations: ", Tulip.get_attribute(tlp, Tulip.BarrierIterations()))
    return nothing
end

if abspath(PROGRAM_FILE) == @__FILE__
    ϵ = parse(Double64, ARGS[2])
    tulip_d64(ARGS[1], ϵ)
end

max_stars_repo_name: mtanneau/LPBenchmarks

max_issues_repo_name: mtanneau/LPBenchmarks

max_forks_repo_name: mtanneau/LPBenchmarks


text:
module DataWriter
  export writeSequences

  import DataReader.FastaRecord

  function writeSequences(output_file_name :: AbstractString, sequences :: Vector{FastaRecord})

    output_file = open(output_file_name, "w")
    for s in sequences
      println(output_file, ">", s.description)
      println(output_file, s.sequence)
    end
    close(output_file)
  end

end


text:
"""
    LIF

A leaky-integrate-fire neuron described by the following differential equation

``\\frac{\\mathrm{d}v}{\\mathrm{d}t} = \\frac{R}{\\tau} I - \\lambda``

# Fields
- `voltage::VT`: membrane potential
- `current::VT`: injected (unprocessed) current
- `lastt::IT`: the last time this neuron processed a spike
- `τm::VT`: membrane time constant
- `vreset::VT`: reset voltage potential
- `R::VT`: resistive constant (typically = 1)
"""
mutable struct LIF{VT<:Real, IT<:Integer} <: AbstractCell
    # required fields
    voltage::VT
    current::VT

    # model specific fields
    lastt::IT
    τm::VT
    vreset::VT
    R::VT
end

Base.show(io::IO, ::MIME"text/plain", neuron::LIF) =
    print(io, """LIF:
                     voltage: $(neuron.voltage)
                     current: $(neuron.current)
                     τm:      $(neuron.τm)
                     vreset:  $(neuron.vreset)
                     R:       $(neuron.R)""")
Base.show(io::IO, neuron::LIF) =
    print(io, "LIF(voltage = $(neuron.voltage), current = $(neuron.current))")

"""
    LIF(τm, vreset, vth, R = 1.0)

Create a LIF neuron with zero initial voltage and empty current queue.
"""
LIF(τm::Real, vreset::Real, R::Real = 1.0) = LIF{Float32, Int}(vreset, 0, 0, τm, vreset, R)

isactive(neuron::LIF, t::Integer; dt::Real = 1.0) = (neuron.current > 0)
getvoltage(neuron::LIF) = neuron.voltage

"""
    excite!(neuron::LIF, current)
    excite!(neurons::AbstractArray{<:LIF}, current)

Excite a `neuron` with external `current`.
"""
excite!(neuron::LIF, current) = (neuron.current += current)
excite!(neurons::T, current) where T<:AbstractArray{<:LIF} = (neurons.current .+= current)

"""
    spike!(neuron::LIF, t::Integer; dt::Real = 1.0)
    spike!(neurons::AbstractArray{<:LIF}, spikes; dt::Real = 1.0)

Record a output spike from the threshold function with the `neuron` body.
Sets `neuron.lastt`.
"""
spike!(neuron::LIF, t::Integer; dt::Real = 1.0) = (t > 0) && (neuron.voltage = neuron.vreset)
function spike!(neurons::T, spikes; dt::Real = 1.0) where T<:AbstractArray{<:LIF}
    map!((x, y, s) -> (s > 0) ? y : x, neurons.voltage, neurons.voltage, neurons.vreset, spikes)
end

"""
    lif(t::Real, I, V; R, tau)
    lif!(t::AbstractArray{<:Real}, I::AbstractArray{<:Real}, V::AbstractArray{<:Real}; R::AbstractArray{<:Real}, tau::AbstractArray{<:Real})
    lif!(t::CuVector{<:Real}, I::CuVector{<:Real}, V::CuVector{<:Real}; R::CuVector{<:Real}, tau::CuVector{<:Real})

Evaluate a leaky integrate-and-fire neuron.
Use `CuVector` instead of `Vector` to evaluate on GPU.

# Fields
- `t`: time since last evaluation in seconds
- `I`: external current
- `V`: current membrane potential
- `vrest`: resting membrane potential
- `R`: resistance constant
- `tau`: time constant
"""
function lif(t::Real, I, V; vrest, R, tau)
    V = V * exp(-t / tau) + vrest
    V += I * (R / tau)

    return V
end
function lif!(t::AbstractArray{<:Real}, I::AbstractArray{<:Real}, V::AbstractArray{<:Real}; vrest::AbstractArray{<:Real}, R::AbstractArray{<:Real}, tau::AbstractArray{<:Real})
    # account for leakage
    @. V = V * exp(-t / tau) + vrest

    # apply update step
    @. V += I * (R / tau)

    return V
end
function lif!(t::CuVector{<:Real}, I::CuVector{<:Real}, V::CuVector{<:Real}; vrest::CuVector{<:Real}, R::CuVector{<:Real}, tau::CuVector{<:Real})
    # account for leakage
    @. V = V * exp(-t / tau) + vrest

    # apply update step
    @. V += I * (R / tau)

    return V
end

"""
    evaluate!(neuron::LIF, t::Integer; dt::Real = 1.0)
    (neuron::LIF)(t::Integer; dt::Real = 1.0)
    evaluate!(neurons::AbstractArray{<:LIF}, t::Integer; dt::Real = 1.0)

Evaluate the neuron model at time `t`. Return the resulting membrane potential.
"""
function evaluate!(neuron::LIF, t::Integer; dt::Real = 1.0)
    neuron.voltage = lif((t - neuron.lastt) * dt, neuron.current, neuron.voltage; vrest = 0, R = neuron.R, tau = neuron.τm)
    neuron.lastt = t
    neuron.current = 0

    return neuron.voltage
end
(neuron::LIF)(t::Integer; dt::Real = 1.0) = evaluate!(neuron, t; dt = dt)
function evaluate!(neurons::T, t::Integer; dt::Real = 1.0) where T<:AbstractArray{<:LIF}
    lif!((t .- neurons.lastt) .* dt, neurons.current, neurons.voltage;
         vrest = zero(neurons.voltage), R = neurons.R, tau = neurons.τm)
    neurons.lastt .= t
    neurons.current .= 0

    return neurons.voltage
end

"""
    reset!(neuron::LIF)
    reset!(neurons::AbstractArray{<:LIF})

Reset `neuron` by setting the membrane potential to `neuron.vreset`.
"""
function reset!(neuron::LIF)
    neuron.lastt = 0
    neuron.voltage = neuron.vreset
    neuron.current = 0
end
function reset!(neurons::T) where T<:AbstractArray{<:LIF}
    neurons.lastt .= 0
    neurons.voltage .= neurons.vreset
    neurons.current .= 0
end


text:
using LinearAlgebra

@inline function exp_gen!(cache,A,::Val{7})
    T=promote_type(eltype(A),Float64) # Make it work for many 'bigger' types (matrices and scalars)
    # max_memslots=6
    n=size(A,1)
    # The first slots are precomputed nodes [:A]
    memslots2 = getmem(cache,2)
    memslots3 = getmem(cache,3)
    memslots4 = getmem(cache,4)
    memslots5 = getmem(cache,5)
    memslots6 = getmem(cache,6)
    # Assign precomputed nodes memslots
    memslots1=A # overwrite A
    # Uniform scaling is exploited.
    # No matrix I explicitly allocated.
    # Computation order: C A2 A4 A6 Ub3 Ub Uc U Vb3 Vb V Z X P S1 S2
    # Computing C = x*A+x*I
    coeff1=0.25
    coeff2=0.0
    # Smart lincomb recycle A
    memslots1 .= coeff1.*memslots1
    inplace_add!(memslots1,I*coeff2)
    # Computing A2 with operation: mult
    mul!(memslots2,memslots1,memslots1)
    # Computing A4 with operation: mult
    mul!(memslots3,memslots2,memslots2)
    # Computing A6 with operation: mult
    mul!(memslots4,memslots2,memslots3)
    # Computing Ub3 = x*A2+x*A4+x*A6
    coeff1=4.08408e7
    coeff2=16380.0
    coeff3=1.0
    memslots5 .= coeff1.*memslots2 .+ coeff2.*memslots3 .+ coeff3.*memslots4
    # Computing Ub with operation: mult
    mul!(memslots6,memslots5,memslots4)
    # Deallocating Ub3 in slot 5
    # Computing Uc = x*I+x*A2+x*A4+x*A6+x*Ub
    coeff1=3.238237626624e16
    coeff2=1.1873537964288e15
    coeff3=1.05594705216e13
    coeff4=3.352212864e10
    coeff5=1.0
    # Smart lincomb recycle Ub
    memslots6 .= coeff2.*memslots2 .+ coeff3.*memslots3 .+ coeff4.*memslots4 .+ coeff5.*memslots6
    inplace_add!(memslots6,I*coeff1)
    # Computing U with operation: mult
    mul!(memslots5,memslots1,memslots6)
    # Deallocating C in slot 1
    # Deallocating Uc in slot 6
    # Computing Vb3 = x*A2+x*A4+x*A6
    coeff1=1.32324192e9
    coeff2=960960.0
    coeff3=182.0
    memslots1 .= coeff1.*memslots2 .+ coeff2.*memslots3 .+ coeff3.*memslots4
    # Computing Vb with operation: mult
    mul!(memslots6,memslots1,memslots4)
    # Deallocating Vb3 in slot 1
    # Computing V = x*I+x*A2+x*A4+x*A6+x*Vb
    coeff1=6.476475253248e16
    coeff2=7.7717703038976e15
    coeff3=1.29060195264e14
    coeff4=6.704425728e11
    coeff5=1.0
    # Smart lincomb recycle A2
    memslots2 .= coeff2.*memslots2 .+ coeff3.*memslots3 .+ coeff4.*memslots4 .+ coeff5.*memslots6
    inplace_add!(memslots2,I*coeff1)
    # Deallocating A4 in slot 3
    # Deallocating A6 in slot 4
    # Deallocating Vb in slot 6
    # Computing Z = x*V+x*U
    coeff1=1.0
    coeff2=-1.0
    memslots1 .= coeff1.*memslots2 .+ coeff2.*memslots5
    # Computing X = x*V+x*U
    coeff1=1.0
    coeff2=1.0
    # Smart lincomb recycle V
    memslots2 .= coeff1.*memslots2 .+ coeff2.*memslots5
    # Deallocating U in slot 5
    # Computing P with operation: ldiv
    ldiv_for_generated!(memslots3, memslots1, memslots2)
    # Deallocating Z in slot 1
    # Deallocating X in slot 2
    # Computing S1 with operation: mult
    mul!(memslots1,memslots3,memslots3)
    # Deallocating P in slot 3
    # Computing S2 with operation: mult
    mul!(memslots2,memslots1,memslots1)
    # Deallocating S1 in slot 1
    copyto!(A,memslots2) # Returning S2
end

max_stars_repo_name: JuliaDiffEq/ExponentialUtilities

max_issues_repo_name: JuliaDiffEq/ExponentialUtilities

max_forks_repo_name: JuliaDiffEq/ExponentialUtilities


text:
function integrate_trapz_7(f, dr)
    N = size(dr, 1)
    g = f .* dr
    s = (  36799 * ( g[1] + g[N  ] ) +
          176648 * ( g[2] + g[N-1] ) +
           54851 * ( g[3] + g[N-2] ) +
          177984 * ( g[4] + g[N-3] ) +
           89437 * ( g[5] + g[N-4] ) +
          130936 * ( g[6] + g[N-5] ) +
          119585 * ( g[7] + g[N-6] )
    ) / 120960
    for i in 8:N-7
        s = s + g[i]
    end
    return s
end


text:
###############
# Status type #
###############

type Status <: GitHubType
    id::Nullable{Int}
    total_count::Nullable{Int}
    state::Nullable{GitHubString}
    description::Nullable{GitHubString}
    context::Nullable{GitHubString}
    sha::Nullable{GitHubString}
    url::Nullable{HttpCommon.URI}
    target_url::Nullable{HttpCommon.URI}
    created_at::Nullable{Dates.DateTime}
    updated_at::Nullable{Dates.DateTime}
    creator::Nullable{Owner}
    repository::Nullable{Repo}
    statuses::Nullable{Vector{Status}}
end

Status(data::Dict) = json2github(Status, data)
Status(id::Real) = Status(Dict("id" => id))

namefield(status::Status) = status.id

###############
# API Methods #
###############

function create_status(repo, sha; options...)
    result = gh_post_json("/repos/$(name(repo))/statuses/$(name(sha))"; options...)
    return Status(result)
end

function statuses(repo, ref; options...)
    results, page_data = gh_get_paged_json("/repos/$(name(repo))/commits/$(name(ref))/statuses"; options...)
    return map(Status, results), page_data
end

function status(repo, ref; options...)
    result = gh_get_json("/repos/$(name(repo))/commits/$(name(ref))/status"; options...)
    return Status(result)
end


text:
using Logging: debug, info, warn
include("log_test.jl")
include("macro_test.jl")
include("test_hierarchy.jl")
include("macro_scope_test.jl")


text:
import PyPlot

for p in [:plot, :semilogx, :semilogy, :loglog]
    q = Expr(:., :PyPlot, QuoteNode(p))
    fn = quote
        function $q(trs::Vector{TestResults}, args...; xfield::Symbol=:scd, yfield::Symbol=:walltime,
                    xunit="", yunit="(s)")
            x = getfield_vec(trs, xfield)
            y = getfield_vec(trs, yfield)
            $q(x,y, args...)
            PyPlot.xlabel(string(xfield)*" "*xunit)
            PyPlot.ylabel(string(yfield)*" "*yunit)
        end
    end

    eval(fn)
end


text:
module FluspectMod
#using GSL
using Polynomials
using Statistics
# Matlab reading
using MAT
# Numerical integration package (Simson rule)
using QuadGK

# Is this OK?
file_Opti = joinpath(dirname(pathof(FluspectMod)), "Optipar2017_ProspectD.mat")

const minwle  = 400.; # PAR range
const maxwle  = 700.;
const minwlf  = 650.; # SIF range
const maxwlf  = 850.;

# Doubling Adding layers
const ndub = 10

# Canopy Layers:
const nl = 60
const xl = collect(0:-1/nl:-1)
const lazitab = collect(5:10:355)
const dx = 1.0/nl

# Read in all optical data:
opti = matread(file_Opti)["optipar"]
nr_     =  convert(Array{Float32}, opti["nr"]);nr = nr_
Km_     =  convert(Array{Float32}, opti["Kdm"]);Km = Km_
Kab_    =  convert(Array{Float32}, opti["Kab"]);Kab = Kab_
Kant_   =  convert(Array{Float32}, opti["Kant"]);Kant =Kant_
Kcar_   =  convert(Array{Float32}, opti["Kca"]);Kcar= Kcar_
Kw_     =  convert(Array{Float32}, opti["Kw"]);Kw=Kw_
KBrown_ =  convert(Array{Float32}, opti["Ks"]);KBrown=KBrown_
phi_    =  convert(Array{Float32}, opti["phi"]);phi=phi_
KcaV_   =  convert(Array{Float32}, opti["KcaV"]);KcaV=KcaV_
KcaZ_   =  convert(Array{Float32}, opti["KcaZ"]);KcaZ =KcaZ_
lambda_ =  convert(Array{Float32}, opti["wl"]);lambda = lambda_

# Enable downsampling spectral resolution to arbitrary grid (specify array with boundaries here)
# Don't like global arrays yet, need to make sure this is not creating performance issues.
function init(swl)
    global WL = swl
    # Stupid...
    global nr = zeros(length(swl)-1)
    global Km = zeros(length(swl)-1)
    global Kab = zeros(length(swl)-1)
    global Kant = zeros(length(swl)-1)
    global Kcar = zeros(length(swl)-1)
    global Kw = zeros(length(swl)-1)
    global KBrown = zeros(length(swl)-1)
    global phi = zeros(length(swl)-1)
    global KcaV = zeros(length(swl)-1)
    global KcaZ = zeros(length(swl)-1)
    global lambda = zeros(length(swl)-1)
    for i in 1:length(swl)-1
        wo = findall((lambda_.>=swl[i]).&(lambda_.<swl[i+1]) )
        #println(mean(nr_[wo]))
        nr[i]   =  mean(nr_[wo])
        Km[i]  =  mean(Km_[wo])
        Kab[i]  =  mean(Kab_[wo])
        Kant[i] =  mean(Kant_[wo])
        Kcar[i]  =  mean(Kcar_[wo])
        Kw[i]   =  mean(Kw_[wo])
        KBrown[i] = mean(KBrown_[wo])
        phi[i]  =  mean(phi_[wo])
        KcaV[i] = mean(KcaV_[wo])
        KcaZ[i] = mean(KcaZ_[wo])
        lambda[i] = mean(lambda_[wo])
    end
end

function fluspect(x::Vector; fqe::Number=0.0, Kab__=Kab, Kant__=Kant, KBrown__=KBrown, Kw__=Kw, Km__=Km, KcaV__=KcaV, KcaZ__=KcaZ)
    #println(fqe)
# , Kab__=Kab, Kant__=Kant, KBrown__=KBrown, Kw__=Kw, Km__=Km, KcaV__=KcaV, KcaZ__=KcaZ)
# ***********************************************************************
# Jacquemoud S., Baret F. (1990), PROSPECT: a model of leaf optical
# properties spectra; Remote Sens. Environ.; 34:75-91.
# Reference:
# Féret, Gitelson, Noble & Jacquemoud [2017]. PROSPECT-D: Towards modeling
# leaf optical properties through a complete lifecycle
# Remote Sensing of Environment; 193:204215
# DOI: http://doi.org/10.1016/j.rse.2017.03.004
# The specific absorption coefficient corresponding to brown pigment is()
# provided by Frederic Baret [EMMAH, INRA Avignon, baret@avignon.inra.fr]
# & used with his autorization.
# ***********************************************************************
N = x[1]
Cab = x[2]
Car = x[3]
Ant = x[4]
Brown = x[5]
Cw = x[6]
Cm = x[7]
Cx = x[8]
Kcaro = (1.0-Cx).* KcaV__ + Cx .* KcaZ__;

Kall    = (Cab*Kab__.+Car*Kcaro.+Ant*Kant__.+Brown*KBrown__.+Cw*Kw__.+Cm*Km__)./N
# Relative absorption by Chlorophyll only (drives SIF and GPP eventually)
kChlrel  = Cab*Kab./(Kall.*N.+eps());

# Adding eps() here to keep it stable and NOT set to 1 manually when Kall=0 (ForwardDiff won't work otherwise)
tau = real.((1.0.-Kall).*exp.(-Kall) .+ Kall.^2.0.*expint.(Kall.+eps()))

# ***********************************************************************
# reflectance & transmittance of one layer
# ***********************************************************************
# Allen W.A., Gausman H.W., Richardson A.J., Thomas J.R. (1969)
# Interaction of isotropic ligth with a compact plant leaf; J. Opt.
# Soc. Am., 59[10]:1376-1379.
# ***********************************************************************
# reflectivity & transmissivity at the interface
#-------------------------------------------------
talf    = calctav.(40.,nr)
ralf    = 1.0.-talf
t12     = calctav.(90.,nr)
r12     = 1.0.-t12
t21     = t12./(nr.^2)
r21     = 1.0.-t21

# top surface side
denom   = 1.0.-r21.*r21.*tau.^2
Ta      = talf.*tau.*t21./denom
Ra      = ralf.+r21.*tau.*Ta

# bottom surface side
t       = t12.*tau.*t21./denom
r       = r12.+r21.*tau.*t

# ***********************************************************************
# reflectance & transmittance of N layers
# Stokes equations to compute properties of next N-1 layers [N real]
# Normal case()
# ***********************************************************************
# Stokes G.G. (1862), On the intensity of the light reflected from
# | transmitted through a pile of plates; Proc. Roy. Soc. Lond.
# 11:545-556.
# ***********************************************************************
D       = sqrt.((1.0.+r.+t).*(1.0.+r.-t).*(1.0.-r.+t).*(1.0.-r.-t))
rq      = r.^2
tq      = t.^2
a       = (1.0.+rq.-tq.+D)./(2r)
b       = (1.0.-rq.+tq.+D)./(2t)

bNm1    = b.^(N.-1);                  #
bN2     = bNm1.^2
a2      = a.^2
denom   = a2.*bN2.-1
Rsub    = a.*(bN2.-1)./denom
Tsub    = bNm1.*(a2.-1)./denom

# Case of zero absorption
j       = findall(r.+t .>= 1)
Tsub[j] = t[j]./(t[j]+(1.0.-t[j])*(N-1))
Rsub[j]	= 1.0.-Tsub[j]

# Reflectance & transmittance of the leaf: combine top layer with next N-1 layers
denom   = 1.0.-Rsub.*r
tran    = Ta.*Tsub./denom
refl    = Ra.+Ta.*Rsub.*t./denom

RT     = [refl tran]
if fqe ==0.0
    return RT
end
# FROM SCOPE notes:
# From here a new path is taken: The doubling method used to calculate
# fluoresence is now only applied to the part of the leaf where absorption
# takes place, that is, the part exclusive of the leaf-air interfaces. The
# reflectance (rho) and transmittance (tau) of this part of the leaf are
# now determined by "subtracting" the interfaces
# CF Note: All of the below takes about 10 times more time than the RT above. Need to rething speed and accuracy. (10nm is bringing it down a lot!)

Rb  = (refl.-ralf)./(talf.*t21+(refl.-ralf).*r21);  # Remove the top interface
Z   = tran.*(1.0.-Rb.*r21)./(talf.*t21);            # Derive Z from the transmittance

rho = (Rb.-r21.*Z.^2)./(1.0.-(r21.*Z).^2);      # Reflectance and transmittance
tau = (1.0.-Rb.*r21)./(1.0.-(r21.*Z).^2).*Z;    # of the leaf mesophyll layer
t   =   tau;
r   =   max.(rho,0.0);                       # Avoid negative r

# Derive Kubelka-Munk s and k
I_rt     =   findall((r.+t).<1);
D[I_rt]  =   sqrt.((1 .+ r[I_rt] .+ t[I_rt]) .* (1 .+ r[I_rt] .- t[I_rt]) .* (1 .- r[I_rt] .+ t[I_rt]) .*  (1 .- r[I_rt] .- t[I_rt]));
a[I_rt]  =   (1 .+ r[I_rt].^2 .- t[I_rt].^2 .+ D[I_rt]) ./ (2r[I_rt]);
b[I_rt]  =   (1 .- r[I_rt].^2 + t[I_rt].^2 .+ D[I_rt]) ./ (2t[I_rt]);
a[(r.+t).>=1] .=   1.0;
b[(r.+t).>=1] .=   1.0;

s        =   r./t;
I_a      =   findall((a.>1).&(a.!=Inf));
s[I_a]   =   2 .*a[I_a] ./ (a[I_a].^2 .- 1) .* log.(b[I_a]);

k        =   log.(b);
k[I_a]   =   (a[I_a].-1) ./ (a[I_a].+1) .* log.(b[I_a]);
kChl     =   kChlrel .* k;

# indices of wle and wlf within wlp
Iwle   = findall((lambda.>=minwle) .& (lambda.<=maxwle));
Iwlf   = findall((lambda.>=minwlf) .& (lambda.<=maxwlf));
wle    = lambda[Iwle];    # excitation wavelengths,  column
wlf    = lambda[Iwlf];    # fluorescence wavelengths,  column
epsi         = 2.0^(-ndub);

# initialisations
te     = 1 .-(k[Iwle].+s[Iwle]) .* epsi;
tf     = 1 .-(k[Iwlf].+s[Iwlf]) .* epsi;
re     = s[Iwle] .* epsi;
rf     = s[Iwlf] .* epsi;

sigmoid     = 1 ./(1 .+exp.(-wlf./10).*exp.(wle'./10));  # matrix computed as an outproduct
#println(size(sigmoid)," ", size(phi), " ", size(kChl)," ", size(Iwle), " ", size(Iwlf), " ", size(kChl[Iwle]))

Mf = Mb = fqe .* ((0.5*phi[Iwlf]).*epsi) .* kChl[Iwle]'.*sigmoid

Ih          = ones(1,length(te));     # row of ones
Iv          = ones(length(tf),1);     # column of ones

# Doubling routine
for i = 1:ndub
    xe = te./(1 .-re.*re);  ten = te.*xe;  ren = re.*(1 .+ten);
    xf = tf./(1 .-rf.*rf);  tfn = tf.*xf;  rfn = rf.*(1 .+tfn);

    A11  = xf*Ih + Iv*xe';           A12 = (xf*xe').*(rf*Ih .+ Iv*re');
    A21  = 1 .+(xf*xe').*(1 .+rf*re');   A22 = (xf.*rf)*Ih+Iv*(xe.*re)';
    #println(size(A11)," ", size(A12), " ", size(Mf)," ", size(Mb), " ")
    Mfn   = Mf  .* A11 .+ Mb  .* A12;
    Mbn   = Mb  .* A21 .+ Mf  .* A22;

    te   = ten;  re  = ren;   tf   = tfn;   rf   = rfn;
    Mf  = Mfn; Mb = Mbn;
end

# Here we add the leaf-air interfaces again for obtaining the final
# leaf level fluorescences.
g = Mb; f = Mf;

Rb = rho .+ tau.^2 .*r21./(1 .-rho.*r21);

Xe = Iv * (talf[Iwle]./(1 .-r21[Iwle].*Rb[Iwle]))';
Xf = t21[Iwlf]./(1 .-r21[Iwlf].*Rb[Iwlf]) * Ih;
Ye = Iv * (tau[Iwle].*r21[Iwle]./(1 .-rho[Iwle].*r21[Iwle]))';
Yf = tau[Iwlf].*r21[Iwlf]./(1 .-rho[Iwlf].*r21[Iwlf]) * Ih;

A = Xe .* (1 .+ Ye.*Yf) .* Xf;
B = Xe .* (Ye .+ Yf) .* Xf;

gn = A .* g + B .* f;
fn = A .* f + B .* g;

Mb  = gn;
Mf  = fn;
return  RT,Mf,Mb
end

function RTM_sail(x::Vector; LIDFa=-0.35, LIDFb=-0.15,q=0.05, tts=30, tto=0, psi=90, wl=lambda, TypeLidf=1)

    # State Vector X includes (in that order):
    #N,Cab,Car,Ant,Cbrown,Cw,Cm,Cx,lai,rsoil
    LAI = x[9]

    # Define soil as polynomial (depends on state vector size):
    pSoil =  Polynomials.Poly(x[10:end])
    rsoil = Polynomials.polyval(pSoil,wl.-mean(wl));

    iLAI    = LAI/nl;               # [1] LAI of elementary layer (guess we can change that)

    # Size of wavelength array:
    nwl = length(wl)

    # Call Fluspect for Leaf optical properties (can be done outside later if needed)
    LRT = fluspect(x[1:8], fqe=0.0)
    ρ=LRT[:,1]
    τ=LRT[:,2]

    # Geometric quantities (need to check allocation cost!)
    cts = cos(deg2rad(tts))
    cto = cos(deg2rad(tto))
    sin_tts  = sin(deg2rad(tts));             #           sin solar       angle
    ctscto = cts*cto;
    tants	= tan(deg2rad(tts));
    tanto	= tan(deg2rad(tto));
    cospsi	= cos(deg2rad(psi));

    dso		= sqrt(tants*tants+tanto*tanto-2.0*tants*tanto*cospsi);

    # Generate leaf angle distribution:
    if TypeLidf==1
        lidf,litab = dladgen(LIDFa,LIDFb);
    elseif TypeLidf==2
        lidf,litab = campbell(LIDFa);
    end
    #println(lidf)

    cos_ttlo    = cos.(deg2rad.(lazitab));  #   cos leaf azimuth angles
    cos_ttli    = cos.(deg2rad.(litab));    #   cosine of normal of upperside of leaf
    sin_ttli    = sin.(deg2rad.(litab));    #   sine   of normal of upperside of leaf

    # angular distance, compensation of shadow length
    # Calculate geometric factors associated with extinction and scattering
    #Initialise sums
    ks	= 0.0;
    ko	= 0.0;
    bf	= 0.0;
    sob	= 0.0;
    sof	= 0.0;

    #	Weighted sums over LIDF
	@simd for i=1:length(litab)
		# ttl = litab[i];	% leaf inclination discrete values
		ctl = cos(deg2rad(litab[i]));
		#	SAIL volume scattering phase function gives interception and portions to be
		#	multiplied by rho and tau
		chi_s,chi_o,frho,ftau=volscatt(tts,tto,psi,litab[i]);

		#********************************************************************************
		#*                   SUITS SYSTEM COEFFICIENTS
		#*
		#*	ks  : Extinction coefficient for direct solar flux
		#*	ko  : Extinction coefficient for direct observed flux
		#*	att : Attenuation coefficient for diffuse flux
		#*	sigb : Backscattering coefficient of the diffuse downward flux
		#*	sigf : Forwardscattering coefficient of the diffuse upward flux
		#*	sf  : Scattering coefficient of the direct solar flux for downward diffuse flux
		#*	sb  : Scattering coefficient of the direct solar flux for upward diffuse flux
		#*	vf   : Scattering coefficient of upward diffuse flux in the observed direction
		#*	vb   : Scattering coefficient of downward diffuse flux in the observed direction
		#*	w   : Bidirectional scattering coefficient
		#********************************************************************************

		#	Extinction coefficients
		ksli = chi_s./cts;
		koli = chi_o./cto;

		#	Area scattering coefficient fractions
		sobli	= frho*pi/ctscto;
		sofli	= ftau*pi/ctscto;
		bfli	= ctl*ctl;
		ks	= ks+ksli*lidf[i];
		ko	= ko+koli*lidf[i];
		bf	= bf+bfli*lidf[i];
		sob	= sob+sobli*lidf[i];
		sof	= sof+sofli*lidf[i];
    end
    #println(sob, " ", sof)
    #	Geometric factors to be used later with rho and tau
	sdb	= 0.5*(ks+bf);
	sdf	= 0.5*(ks-bf);
	dob	= 0.5*(ko+bf);
	dof	= 0.5*(ko-bf);
	ddb	= 0.5*(1 .+bf);
	ddf	= 0.5*(1 .-bf);
    # Skipped SCOPE lines 186-213 here (catch up later)
    # 1.4 solar irradiance factor for all leaf orientations
    # See eq 19 in vdT 2009
    Cs          = cos_ttli.*cts;             # [nli]     pag 305 modified by Joris
    Ss          = sin_ttli.*sin_tts;         # [nli]     pag 305 modified by Joris

    cos_deltas  = Cs*ones(1,length(lazitab)) .+ Ss*cos_ttlo';   # [nli,nlazi]
    fs          = abs.(cos_deltas./cts);         # [nli,nlazi] pag 305

    # 1.5 probabilities Ps, Po, Pso
    Ps          =   exp.(ks*xl*LAI);     # [nl+1]  p154{1} probability of viewing a leaf in solar dir
    Po          =   exp.(ko*xl*LAI);     # [nl+1]  p154{1} probability of viewing a leaf in observation dir

    Ps[1:nl]    =   Ps[1:nl] *(1 .-exp.(-ks*LAI*dx))/(ks*LAI*dx);   # Correct Ps/Po for finite dx
    Po[1:nl]    =   Po[1:nl] *(1 .-exp.(-ko*LAI*dx))/(ko*LAI*dx);   # Correct Ps/Po for finite dx

    #Pso: Probability of observing a sunlit leaf at depth x, see eq 31 in vdT 2009
    Pso         =   similar(Po);
    for j=1:length(xl)
        #println(size(a), " ", size(Pso), " ", size(Po))
        Pso[j] = quadgk(x -> Psofunction(ko,ks,LAI,q,dso,x), xl[j]-dx,xl[j], rtol=1e-2)[1]/dx
        #Pso[j,:]=   quad(@(y)Psofunction(K,k,LAI,q,dso,y),xl(j)-dx,xl(j))/dx; %#ok<FREMO>
    end

    Pso[Pso.>Po]= minimum([Po[Pso.>Po] Ps[Pso.>Po]],dims=2);    #takes care of rounding error
    Pso[Pso.>Ps]= minimum([Po[Pso.>Ps] Ps[Pso.>Ps]],dims=2);    #takes care of rounding error

    # All with length of wavelengths:
    sigb = ddb.*ρ.+ddf.*τ;
	sigf = ddf.*ρ.+ddb.*τ;
    sb   = sdb*ρ .+ sdf*τ;            # [nwl]     sb,     p305{1} diffuse     backscatter scattering coefficient for specular incidence
    sf   = sdf*ρ .+ sdb*τ;            # [nwl]     sf,     p305{1} diffuse     forward     scattering coefficient for specular incidence
    vb   = dob*ρ .+ dof*τ;            # [nwl]     vb,     p305{1} directional backscatter scattering coefficient for diffuse  incidence
    vf   = dof*ρ .+ dob*τ;            # [nwl]     vf,     p305{1} directional forward     scattering coefficient for diffuse  incidence
    w    = sob*ρ .+ sof*τ;            # [nwl]     w,      p309{1} bidirectional scattering coefficent (directional-directional)
    a    = 1 .-sigf;                  # [nwl]     attenuation
    m    = sqrt.(a.^2 .-sigb.^2);     # [nwl]
    rinf = (a.-m)./sigb;              # [nwl]
    rinf2= rinf.*rinf;                # [nwl]
    #println(minimum(m), " ", min(ks))
    # direct solar radiation
    J1k        = calcJ1.(-1, m,ks,LAI);          # [nwl]
    J2k        = calcJ2.( 0, m,ks,LAI);          # [nwl]
    J1K        = calcJ1.(-1, m,ko,LAI);          # [nwl]   % added for calculation of rdo
    J2K        = calcJ2.( 0, m,ko,LAI);          # [nwl]   % added for calculation of rdo
    e1          = exp.(-m.*LAI);                 # [nwl]
    e2          = e1.^2;                         # [nwl]
    re          = rinf.*e1;                      # [nwl]

    denom       = 1 .-rinf2.*e2;
    s1          = sf .+rinf.*sb;
    s2          = sf.*rinf+sb;
    v1          = vf.+rinf.*vb;
    v2          = vf.*rinf.+vb;

    Pss         = s1.*J1k;          # [nwl]
    Qss         = s2.*J2k;          # [nwl]

    Poo         = v1.*J1K;          # (nwl)   % added for calculation of rdo
    Qoo         = v2.*J2K;          # [nwl]   % added for calculation of rdo

    tau_ss      = exp(-ks*LAI);                  # [1]
    tau_oo      = exp(-ko*LAI);

    Z           = (1 - tau_ss * tau_oo)/(ks + ko);  # needed for analytic rso

    tau_dd      = (1 .-rinf2).*e1 ./denom;        # [nwl]
    rho_dd      = rinf.*(1 .-e2)  ./denom;        # [nwl]
    tau_sd      = (Pss.-re.*Qss) ./denom;         # [nwl]
    tau_do      = (Poo.-re.*Qoo) ./denom;         # [nwl]
    rho_sd      = (Qss.-re.*Pss) ./denom;         # [nwl]
    rho_do      = (Qoo.-re.*Poo) ./denom;         # (nwl)

    T1          = v2.*s1.*(Z.-J1k*tau_oo)./(ko.+m).+v1.*s2.*(Z.-J1K*tau_ss)./(ks.+m);
    T2          = -(Qoo.*rho_sd+Poo.*tau_sd).*rinf;
    rho_sod     = (T1+T2)./(1 .-rinf2);

    #	Bidirectional reflectance
    #	Single scattering contribution
    rho_sos = w.*iLAI.*sum(Pso[1:nl]);
    #	Total canopy contribution
    rho_so=rho_sos.+rho_sod;
    #println(rho_so[100:120])

    dn=1 .-rsoil.*rho_dd;
    # Total canopy contribution
    rso         = rho_so .+ rsoil .* Pso[nl+1] .+ ((tau_sd.+tau_ss*rsoil.*rho_dd).*tau_oo.+(tau_sd.+tau_ss).*tau_do).*rsoil./denom;
    # SAIL analytical reflectances
    # rsd: directional-hemispherical reflectance factor for solar incident flux
    rsd     = rho_sd .+ (tau_ss .+ tau_sd).*rsoil.*tau_dd./denom;
    # rdd: bi-hemispherical reflectance factor
    rdd     = rho_dd .+ tau_dd.*rsoil.*tau_dd./denom;
    # rdo: bi-directional reflectance factor
    rdo     = rho_do .+ (tau_oo .+ tau_do).*rsoil.*tau_dd./denom;
    #return [rso rsd rdd rdo

    # Dummy code here to track direct and diffuse light first, need to separate this into another function later
    Esun_ = zeros(nwl).+100
    Esky_ = zeros(nwl).+100
    Emin_ = zeros(nl+1,nwl)
    Eplu_ = zeros(nl+1,nwl)
    Eplu_1      = rsoil.*((tau_ss.+tau_sd).*Esun_.+tau_dd.*Esky_)./denom;
    Eplu0       = rho_sd.*Esun_ .+ rho_dd.*Esky_ .+ tau_dd.*Eplu_1;
    Emin_1      = tau_sd.*Esun_ .+ tau_dd.*Esky_ .+ rho_dd.*Eplu_1;
    delta1      = Esky_  .- rinf.*Eplu0;
    delta2      = Eplu_1 .- rinf.*Emin_1;

    # calculation of the fluxes in the canopy (this seems slow!)
    t1 = sf.+rinf.*sb
    t2 = sb+rinf.*sf
    # The order here mattered, now doing the loop over nl, not nwl! (faster)
    # This loop is time consuming, probably don't need high spectral resolution for the NIR part here, so it can be shortened a lot (say 100nm steps from 700-2500nm?)
    # We just need it for net energy balance and PAR here.
    return [rso rsd rdd rdo]
    for i = 1:nl+1
        J1kx  = calcJ1.(xl[i],m,ks,LAI);    #           [nl]
        J2kx  = calcJ2.(xl[i],m,ks,LAI);    #           [nl]
        F1    = Esun_.*J1kx.*t1 .+ delta1.*exp.( m.*LAI.*xl[i]);      #[nl]
        F2    = Esun_.*J2kx.*t2 .+ delta2.*exp.(-m.*LAI.*(xl[i].+1));  #[nl]
        Emin_[i,:]  = (F1.+rinf.*F2)./(1 .-rinf2);#        [nl,nwl]
        Eplu_[i,:]  = (F2.+rinf.*F1)./(1 .-rinf2);#        [nl,nwl]
    end

end

# Had to make this like the Prosail F90 function, forwardDiff didn't work otherwise.
function calcJ1(x,m,k,LAI)
    del=(k-m)*LAI
    if abs(del)>1E-3;
        J1 = (exp(m*LAI*x)-exp(k*LAI*x))/(k-m);
    else
        J1 = -0.5*LAI*x*(exp(m*LAI*x)+exp(k*LAI*x))*(1.0-del^2/12.0);
    end
    return J1
end

function calcJ2(x,m,k,LAI)
    return (exp(k*LAI*x)-exp(-k*LAI)*exp(-m*LAI*(1+x)))/(k+m);
end

# APPENDIX IV function Pso from SCOPE v1.73
function Psofunction(K,k,LAI,q,dso,xl)
    if dso!=0.0
        alf         =   (dso/q) *2/(k+K);
        pso         =   exp((K+k)*LAI*xl + sqrt(K*k)*LAI/(alf  )*(1-exp(xl*(alf  ))));# [nl+1]  factor for correlation of Ps and Po
    else
        pso         =   exp((K+k)*LAI*xl - sqrt(K*k)*LAI*xl);# [nl+1]  factor for correlation of Ps and Po
    end
    return pso
end

# FROM SCOPE v1.73 APPENDIX II function volscat
"""********************************************************************************
!*	tts		= solar zenith
!*	tto		= viewing zenith
!*	psi		= azimuth
!*	ttl		= leaf inclination angle
!*	chi_s	= interception functions
!*	chi_o	= interception functions
!*	frho	= function to be multiplied by leaf reflectance rho
!*	ftau	= functions to be multiplied by leaf transmittance tau
!********************************************************************************

!	Compute volume scattering functions and interception coefficients
!	for given solar zenith, viewing zenith, azimuth and leaf inclination angle.

!	chi_s and chi_o are the interception functions.
!	frho and ftau are the functions to be multiplied by leaf reflectance rho and
!	leaf transmittance tau, respectively, in order to obtain the volume scattering
!	function.
"""
function volscatt(tts,tto,psi,ttli)
    #Volscatt version 2.
    #created by W. Verhoef
    #edited by Joris Timmermans to matlab nomenclature.
    # date: 11 February 2008
    #tts    [1]         Sun            zenith angle in degrees
    #tto    [1]         Observation    zenith angle in degrees
    #psi    [1]         Difference of  azimuth angle between solar and viewing position
    #ttli   [ttli]      leaf inclination array

    nli     = length(ttli);

    psi_rad         = deg2rad(psi);
    cos_psi         = cos(deg2rad(psi));                #   cosine of relative azimuth angle
    cos_ttli        = cos(deg2rad(ttli));               #   cosine of normal of upperside of leaf
    sin_ttli        = sin(deg2rad(ttli));               #   sine   of normal of upperside of leaf
    cos_tts         = cos(deg2rad(tts));                #   cosine of sun zenith angle
    sin_tts         = sin(deg2rad(tts));                #   sine   of sun zenith angle
    cos_tto         = cos(deg2rad(tto));                #   cosine of observer zenith angle
    sin_tto         = sin(deg2rad(tto));                #   sine   of observer zenith angle
    Cs              = cos_ttli*cos_tts;                 #   p305{1}
    Ss              = sin_ttli*sin_tts;                 #   p305{1}
    Co              = cos_ttli*cos_tto;                 #   p305{1}
    So              = sin_ttli*sin_tto;                 #   p305{1}
    As              = maximum([Ss,Cs]);
    Ao              = maximum([So,Co]);
    #println(-Cs./As, " ", Ss, " ", Cs)
    bts             = acos.(-Cs./As);                    #   p305{1}
    bto             = acos.(-Co./Ao);                    #   p305{2}
    chi_o           = 2/pi*((bto-pi/2).*Co + sin(bto).*So);
    chi_s           = 2/pi*((bts-pi/2).*Cs + sin(bts).*Ss);
    delta1          = abs(bts-bto);                     #   p308{1}
    delta2          = pi-abs(bts + bto - pi);           #   p308{1}

    Tot             = psi_rad + delta1 + delta2;        #   pag 130{1}

    bt1             = minimum([psi_rad,delta1]);
    bt3             = maximum([psi_rad,delta2]);
    bt2             = Tot - bt1 - bt3;

    T1              = 2Cs.*Co + Ss.*So.*cos_psi;
    T2              = sin(bt2).*(2As.*Ao + Ss.*So.*cos(bt1).*cos(bt3));

    Jmin            = (   bt2).*T1 - T2;
    Jplus           = (pi-bt2).*T1 + T2;

    frho            =  Jplus/(2pi^2);
    ftau            = -Jmin /(2pi^2);

    # pag.309 wl-> pag 135{1}
    frho            = maximum([0.0,frho]);
    ftau            = maximum([0.0,ftau]);
    #println(tts, " ",tto, " ",psi, " ",ttli)
    #println(chi_s, " ", chi_o, " ",frho, " ",ftau)
    return chi_s,chi_o,frho,ftau
end

function dladgen(a::Number,b::Number)
    litab=[5.,15.,25.,35.,45.,55.,65.,75.,81.,83.,85.,87.,89.];
    freq = similar(litab)
    for i1=1:8
        t   = i1*10;
        freq[i1]=dcum(a,b,t);
    end
    for i2=9:12
        t   = 80.0+(i2-8)*2.;
        freq[i2]=dcum(a,b,t);
    end

    freq[13]=1;
    for i   = 13:-1:2
        freq[i]=freq[i]-freq[i-1];
    end
    return freq,litab
end

function dcum(a::Number,b::Number,t::Number)
    y = 0.0
    if a>=1
        f = 1-cos(deg2rad(t));
    else
        epsi=1e-8;
        delx=1;
        x=2*deg2rad(t);
        p=x;
    	while (delx >= epsi)
            #println(delx)
            y = a*sin(x)+0.5*b*sin(2.0*x);
            dx=0.5*(y-x+p);
            x=x+dx;
            delx=abs(dx);
        end
    	f = (2.0*y+p)/pi;
    end
    return f
end

"""
From SCOPE v1.73:
********************************************************************************
*                          Campbell.f
*
*    Computation of the leaf angle distribution function value (freq)
*    Ellipsoidal distribution function caracterised by the average leaf
*    inclination angle in degree (ala)
*    Campbell 1986
*
********************************************************************************
 edit 2017 12 28: change sampling of angles to match with dladgen.m
"""
function campbell(ala)
    tx1=[10.,20.,30.,40.,50.,60.,70.,80.,82.,84.,86.,88.,90.];
    tx2=[0.,10.,20.,30.,40.,50.,60.,70.,80.,82.,84.,86.,88.];

    litab   = (tx2.+tx1)./2.0;
    n=length(litab);
    tl1     = deg2rad(tx1)
    tl2     = deg2rad(tx2)
    excent  = exp(-1.6184e-5*ala^3+2.1145e-3*ala^2-1.2390e-1*ala+3.2491);
    sum0    = 0;

    freq=zeros(n);
    for i=1:n
        x1  = excent./(sqrt(1 .+excent^2 .*tan(tl1(i)).^2));
        x2  = excent./(sqrt(1 .+excent^2 .*tan(tl2(i)).^2));
        if (excent==1)
            freq[i] = abs(cos(tl1(i))-cos(tl2(i)));
        else
            alpha  = excent./sqrt(abs(1-excent.^2));
            alpha2 = alpha.^2;
            x12 = x1.^2;
            x22 = x2.^2;
            if (excent>1)
                alpx1 = sqrt(alpha2(excent>1)+x12(excent>1));
                alpx2[excent>1] = sqrt(alpha2(excent>1)+x22(excent>1));
                dum   = x1*alpx1+alpha2*log(x1+alpx1);
                freq[i] = abs(dum-(x2.*alpx2+alpha2.*log(x2+alpx2)));
            else
                almx1 = sqrt(alpha2-x12);
                almx2 = sqrt(alpha2-x22);
                dum   = x1.*almx1+alpha2.*asin(x1./alpha);
                freq[i] = abs(dum-(x2.*almx2+alpha2.*asin(x2./alpha)));
            end
        end
    end
    sum0 = sum(freq,dims=2);
    freq0=freq./sum0;
return freq0,litab
end

"""
    calctav(alfa, nr)

    ***********************************************************************
    From calctav.m in PROSPECT-D
    ***********************************************************************
    Stern F. (1964), Transmission of isotropic radiation across an
    interface between two dielectrics, Appl. Opt., 3(1):111-113.
    Allen W.A. (1973), Transmission of isotropic light across a
    dielectric surface in two and three dimensions, J. Opt. Soc. Am.,
    63(6):664-666.
    ***********************************************************************

"""
function calctav(α,nr)
    rd  = pi/180
    n2  = nr^2
    np  = n2+1
    nm  = n2-1
    a   = (nr+1)*(nr+1)/2
    k   = -(n2-1)*(n2-1)/4
    sa  = sin(α*rd)
    if α!=90.0
        b1  = sqrt((sa^2-np/2).*(sa^2-np/2)+k)
    else
        b1 = 0
    end
    b2  = sa^2-np/2
    b   = b1-b2
    b3  = b^3
    a3  = a^3
    ts  = (k^2.0/(6*b3)+k/b-b/2)-(k^2.0/(6*a3)+k/a-a/2)
    tp1 = -2*n2*(b-a)/(np^2)
    tp2 = -2*n2*np*log(b/a)/(nm^2)
    tp3 = n2*(1.0/b-1.0/a)/2
    tp4 = 16*n2^2.0*(n2^2+1)*log((2*np*b-nm^2)/(2*np*a-nm^2))/(np^3.0*nm^2)
    tp5 = 16*n2^3.0*(1.0/(2*np*b-nm^2)-1.0/(2*np*a-nm^2))/(np^3)
    tp  = tp1+tp2+tp3+tp4+tp5
    tav = (ts+tp)/(2*sa^2)
    return tav
end

#function expint(x)
#    p = Poly([8.267661952366478e+00, -7.773807325735529e-01, -3.012432892762715e-01, -7.811863559248197e-02, -1.019573529845792e-02,-6.973790859534190e-04,-2.569498322115933e-05, -4.819538452140960e-07,  -3.602693626336023e-09])
#    polyv = polyval(p,real(x));
#    k = findall( abs(imag(x)) <= polyv );
#    -GSL.sf_expint_Ei(-x)
#end
# From Matlab!
function expint(x)
    p = Polynomials.Poly([8.267661952366478e+00, -7.773807325735529e-01, -3.012432892762715e-01, -7.811863559248197e-02, -1.019573529845792e-02,-6.973790859534190e-04,-2.569498322115933e-05, -4.819538452140960e-07,  -3.602693626336023e-09])
    polyv = Polynomials.polyval(p,real(x));
    if abs(imag(x)) <= polyv
        #initialization
        egamma=0.57721566490153286061;
        xk = x;
        yk = -egamma - log(xk);
        j = 1;
        pterm = xk;
        term = xk;
        while abs(term) > (eps(yk))
            yk = yk + term;
            j = j + 1;
            pterm = -xk.*pterm/j;
            term = pterm/j;
        end # end of the while loop
        y = yk;
    else
        n = 1.0;
        xk = x;
        am2 = 0.0
        bm2 = 1.0
        am1 = 1.0
        bm1 = xk
        f = am1 / bm1;
        oldf = Inf;
        j = 2;
        while abs(f-oldf) > (100*eps()*abs(f))
            alpha = n-1+(j/2); # note: beta= 1
            #calculate A(j), B(j), and f(j)
            a = am1 + alpha * am2;
            b = bm1 + alpha * bm2;

            # save new normalized variables for next pass through the loop
            #  note: normalization to avoid overflow or underflow
            am2 = am1 / b;
            bm2 = bm1 / b;
            am1 = a / b;
            bm1 = 1.0;

            f = am1;
            j = j+1;

            # calculate the coefficients for j odd
            alpha = (j-1)/2;
            beta = xk;
            a = beta * am1 + alpha * am2;
            b = beta * bm1 + alpha * bm2;
            am2 = am1 / b;
            bm2 = bm1 / b;
            am1 = a / b;
            bm1 = 1;
            oldf = f;
            f = am1;
            j = j+1;
        end
        y= exp(-xk) * f - 1im*pi*((real(xk)<0)&(imag(xk)==0));
    end
    return y
end

end

max_stars_repo_name: ruillercoaster/LSM-SPAM

max_issues_repo_name: ruillercoaster/LSM-SPAM

max_forks_repo_name: ruillercoaster/LSM-SPAM


text:
using JuMP, Gurobi, Distributions, StatsFuns

const gurobi_env = Gurobi.Env()

# DIMENSIONS
const numAssets = 1000
const numRandomVariables = numAssets
const numVariables = numAssets

const returnLBD = 1.2 #lower bound on the return

# CASES
# number of samples, replicates for scenario approximations
const NumSamples = ceil.(Int64,logspace(1,5,50))
const NumReplicates = 20

#*============ OTHER MODEL PARAMETERS ================
const mean = Array{Float64}(numAssets)
const stdev = Array{Float64}(numAssets)
for i = 1:numAssets
	mean[i] = 1.05 + 0.3*(numAssets-i)/(numAssets-1)
	stdev[i] = (0.05 + 0.6*(numAssets-i)/(numAssets-1))/3
end
const sigma = zeros(Float64, numAssets, numAssets)
for i = 1:numAssets
	sigma[i,i] = stdev[i]
end
#*===============================================

# FILENAMES
const objFile = "objectiveValue.txt"
const solnFile = "solution.txt"
const riskFile = "riskLevel.txt"
const solnTimeFile = "solutionTime.txt"
const checkingTimeFile = "checkingTime.txt"
const sortingTimeFile = "sortingTime.txt"
const callbackFile = "numCallbacks.txt"
const numConstraintsFile = "numConstraintsEnforced.txt"

# compute negative probability of all constraints being satisfied
# for given scenarios of demands d and a given decision vector x
function computeRiskLevel(y_opt::Array{Float64})

	riskLevel::Float64 = 1.0 - normcdf((sum(mean[i]*y_opt[i] for i = 1:numAssets) - returnLBD)/norm(sigma*y_opt))

	return riskLevel
end

# generate random samples of the demands
function generateRandomSamples(numSamples::Int64)

	xi = Array{Float64}(numAssets,numSamples)
	for i = 1:numAssets
		xi[i,:] = rand(Normal(mean[i],stdev[i]), numSamples)
	end

	return xi
end

# determines if all constraints are satisfied for a given single scenario
# returns one if scenario constraints are all satisfied, zero otherwise
function getScenarioConstraintViolations(x::Array{Float64},xi::Array{Float64},pickFirstSetOfConstraints::Bool,maxNumConstraintsPerIteration::Int64)

	const numSamples = size(xi,2)
	numViolated::Int64 = 0
	constraintViolations = Float64[]
	constraintViolationIndices = Int64[]
	constraintViolationTolerance::Float64 = 1E-06

	for iter = 1:numSamples
		con::Float64 = returnLBD - xi[:,iter]'*x

		if(con >= constraintViolationTolerance)
			numViolated += 1
			push!(constraintViolations,con)
			push!(constraintViolationIndices,iter)
			if(pickFirstSetOfConstraints && numViolated == maxNumConstraintsPerIteration)
				break
			end
		end
	end

	return numViolated, constraintViolations, constraintViolationIndices
end

# write solution to file
function initializeFiles(details_file::String)

	writeMode::String = "w"

	open(details_file, writeMode) do f
		write(f,"Model: $modelName \n")
		write(f,"numAssets: $numAssets \n")
		write(f,"returnLBD: $returnLBD \n")
		write(f,"NumSamples: $NumSamples \n")
		write(f,"NumReplicates: $NumReplicates \n")
		write(f,"Method: Scenario approximations \n")
	end

end

# write solution to file
function writeSolutionToFile(dirName::String, status::Symbol, objvar_opt::Float64, x_opt::Array{Float64}, riskLevel::Float64, solutionTime::Float64, checkingTime::Float64, sortingTime::Float64, numCallbacks::Int64, numConstraintsAdded::Int64)

	if(status == Symbol("Optimal"))
		writeMode::String = "a"

		obj_file::String = dirName * objFile
		open(obj_file, writeMode) do f
			write(f,"$objvar_opt \n")
		end

		solution_file::String = dirName * solnFile
		open(solution_file, writeMode) do f
			for i = 1:numVariables
				write(f,"$(x_opt[i]) ")
			end
			write(f,"\n")
		end

		risk_file::String = dirName * riskFile
		open(risk_file, writeMode) do f
			write(f,"$riskLevel \n")
		end

		solntime_file::String = dirName * solnTimeFile
		open(solntime_file, writeMode) do f
			write(f,"$solutionTime \n")
		end

		checktime_file::String = dirName * checkingTimeFile
		open(checktime_file, writeMode) do f
			write(f,"$checkingTime \n")
		end

		sorttime_file::String = dirName * sortingTimeFile
		open(sorttime_file, writeMode) do f
			write(f,"$sortingTime \n")
		end

		callback_file::String = dirName * callbackFile
		open(callback_file, writeMode) do f
			write(f,"$numCallbacks \n")
		end

		constraints_file::String = dirName * numConstraintsFile
		open(constraints_file, writeMode) do f
			write(f,"$numConstraintsAdded \n")
		end
	end

end

# solve the scenario approximation model
function solveScenarioApproximationModel(numScenarios::Int64,dirName::String)

	tic()

	# sample the random variables
	xi = generateRandomSamples(numScenarios)

	maxNumConstraintsPerIteration::Int64 = 100000
	numConstraintsAdded::Int64 = 0

	pickFirstSetOfConstraints::Bool = false

	sortingTime::Float64 = 0.0
	checkingTime::Float64 = 0.0
	numCallbacks::Int64 = 0

	objvar_opt::Float64 = Inf
	y_opt = zeros(Float64,numVariables)
	status::Symbol = Symbol("Error")

	numViolations::Int64 = numScenarios
	constraintIndices = Int64[]

	while(numViolations > 0)

		numCallbacks += 1
		numViolations = 0

		numConstraintsConsidered::Int64 = size(constraintIndices)[1]

		# construct the scenario model
		mod = Model(solver=GurobiSolver(gurobi_env,Presolve=0,OutputFlag=0))

		@variable(mod, y[1:numAssets] >= 0)

		@objective(mod, Min, sum(abs2(50.0*sigma[i,i])*y[i]*y[i] for i = 1:numAssets))
		@constraint(mod, sum(y[i] for i = 1:numAssets) == 1)
		@constraint(mod, [scen=1:numConstraintsConsidered], sum(xi[i,scen]*y[i] for i = 1:numAssets) >= returnLBD)

		status = solve(mod)
		y_opt = getvalue(y)
		objvar_opt = sum(abs2(sigma[i,i]*y_opt[i]) for i = 1:numAssets)

		numViolatedConstraints::Int64 = 0
		constraintViolations = Float64[]
		constraintViolationIndices = Int64[]

		tic()

		numViolatedConstraints, constraintViolations, constraintViolationIndices = getScenarioConstraintViolations(y_opt,xi,pickFirstSetOfConstraints,maxNumConstraintsPerIteration)

		checkingTime += toq()

		numConstraintsEnforced::Int64 = min(size(constraintViolations,1),maxNumConstraintsPerIteration)

		tic()
		relevantIndices = 1:1:numConstraintsEnforced
		if(!pickFirstSetOfConstraints)
			relevantIndices = sortperm(constraintViolations, rev=true)
		end
		sortingTime += toq()

		for iter = 1:numConstraintsEnforced
			push!(constraintIndices,constraintViolationIndices[relevantIndices[iter]])
		end

		numConstraintsAdded += numConstraintsEnforced

		numViolations += numViolatedConstraints

	end

	riskLevel::Float64 = computeRiskLevel(y_opt)

	solutionTime::Float64 = toq()

	@printf "  Objective: %1.5f,  Risklevel: %.5f,  Time: %4.2f,  NumCallbacks: %d,  NumConstraints: %d \n" objvar_opt riskLevel solutionTime numCallbacks numConstraintsAdded

	writeSolutionToFile(dirName, status, objvar_opt, y_opt, riskLevel, solutionTime, checkingTime, sortingTime, numCallbacks, numConstraintsAdded)

	return status, objvar_opt, y_opt, riskLevel, solutionTime, checkingTime, sortingTime, numCallbacks, numConstraintsAdded
end

max_stars_repo_name: rohitkannan/SA-for-CCP

max_issues_repo_name: rohitkannan/SA-for-CCP

max_forks_repo_name: rohitkannan/SA-for-CCP