text:
function firstorder_optimality_tangnorm(pb::CompositeProblem, x, M, ∇f_x)
    x_man = x
    if !isa(x, SVDMPoint)
        x_man = project(M, x)
    end
    !is_manifold_point(M, x) && @debug "trouble here: !is manifold point" x M
    return firstorder_optimality_tangnorm(pb::CompositeProblem, pb.regularizer, x_man, M, ∇f_x)
end

function compute_minnorm_subgradient(pb::CompositeProblem, regularizer, x, M, ∇f_x)
    # model = Model(with_optimizer(
    #     Mosek.Optimizer,
    #     QUIET = true,
    #     INTPNT_CO_TOL_REL_GAP = 1e-15,
    #     INTPNT_CO_TOL_PFEAS = 1e-13,
    #     INTPNT_CO_TOL_DFEAS = 1e-13,
    # ))

    model = Model(with_optimizer(
        OSQP.Optimizer
    ))

    ḡ_normal = model_g_subgradient!(model, regularizer, M, x)
    ḡ = build_subgradient_from_normalcomp(regularizer, M, x, ḡ_normal)

    ḡres = zeros(pb.n1, pb.n2)
    if length(ḡ_normal) == 0
        ḡres .= convert.(Float64, ḡ)
        # return convert.(Float64, ḡ)
    else
        @objective(model, Min, 0.5*dot(ḡ+∇f_x, ḡ+∇f_x))

        optimize!(model)
        # @show termination_status(model)

        ḡres .= value.(ḡ)
        # return value.(ḡ)
    end

    @show manifold_dimension(M)
    display(ḡres)

    ḡexpl = zeros(pb.n1, pb.n2)
    # for i in 1:pb.n
    #     # @show i, i in M.nnz_coords, ∇f_x[i]
    #     if M.nnz_coords[i]
    #         # printstyled("$i\n", color=:red)
    #         ḡexpl[i] = pb.regularizer.λ * sign(x[i])
    #     else
    #         # printstyled("$i\n", color=:blue)
    #         ḡexpl[i] = min(pb.regularizer.λ, max(-pb.regularizer.λ, -∇f_x[i]))
    #     end
    # end

    # @show norm(ḡres - ḡexpl)
    # @show M
    # @show M.nnz_coords
    # display(ḡres)
    # display(ḡexpl)

    if norm(ḡres - ḡexpl) > 1e-7
        @error "--"
        @show norm(ḡres - ḡexpl)
        @show M
        display(ḡres)
        display(ḡexpl)
        display(-∇f_x)
    end

    if manifold_dimension(M) < 9
        @assert false
    end
    return ḡres
end

#
### l1 regularizer
#
function firstorder_optimality_tangnorm(pb::CompositeProblem, regularizer::regularizer_l1, x, M, ∇f_x)
    ḡ = compute_minnorm_subgradient(pb, regularizer, x, M, ∇f_x)

    fḡ_tangent = project(M, x, ∇f_x+ḡ)
    fḡ_normal = ∇f_x+ḡ - fḡ_tangent

    grad_fgₖ = egrad_to_rgrad(M, x, ∇f_x) + ∇M_g(pb, M, x)

    if !isapprox(norm(fḡ_tangent), norm(M, x, grad_fgₖ))
        @debug "riemannian gradient and tangent subdiff mismatch" norm(fḡ_tangent) norm(M, x, grad_fgₖ) norm(project(M, x, ḡ)) norm(M, x, ∇M_g(pb, M, x))
    end

    return norm(fḡ_tangent), dist_int_subdiff(pb, x, M, ḡ)
end

function compute_minnorm_subgradient(pb::CompositeProblem, regularizer::regularizer_l1, x, M, ∇f_x)
    ḡ = zeros(pb.n)
    for i in 1:pb.n
        if M.nnz_coords[i]
            ḡ[i] = pb.regularizer.λ * sign(x[i])
        else
            ḡ[i] = min(pb.regularizer.λ, max(-pb.regularizer.λ, -∇f_x[i]))
        end
    end
    return ḡ
end

function dist_int_subdiff(pb, x_man, M, ḡ)
    if manifold_dimension(M) == problem_dimension(pb)
        return 0
    end

    ḡ_normal = build_normalcomp_from_subgradient(pb.regularizer, M, x_man, ḡ)
    res = maximum(dist_int_subdiff_l1.(ḡ_normal, pb.regularizer.λ)) / pb.regularizer.λ

    return res
end

dist_int_subdiff_l1(x, λ) = max(x-λ, -x-λ)

#
### Nuclear regularizer
#
function firstorder_optimality_tangnorm(pb::CompositeProblem, regularizer::regularizer_lnuclear, x, M::FixedRankMatrices{m,n,k}, ∇f_x) where {m,n,k}
    #! normal component not implemented here.

    ḡ_tan = regularizer.λ .* x.U[:, 1:k] * x.Vt[1:k, :]
    # fḡ_tangent = project(M, x, ∇f_x+ḡ_tan)

    return norm(project(M, x, ∇f_x) + ḡ_tan), -1
end


text:
import Distributions: logpdf, pdf

struct SDT{T1,T2} <: ContinuousUnivariateDistribution
    d::T1
    c::T2
end

logpdf(d::SDT,data::Vector{Int64}) = logpdf(d,data...)

logpdf(d::SDT,data::Tuple{Vararg{Int64}}) = logpdf(d,data...)

function logpdf(d::SDT,hits,fas,Nd)
    @unpack d,c=d
    θhit=cdf(Normal(0,1),d/2-c)
    θfa=cdf(Normal(0,1),-d/2-c)
    loghits = logpdf(Binomial(Nd,θhit),hits)
    logfas = logpdf(Binomial(Nd,θfa),fas)
    return loghits+logfas
end

pdf(d::SDT,data::Vector{Int64}) = exp(logpdf(d,data...))

ext: jl


text:
## MCTS.m is a monte carlos tree search that runs in ColorGraph.m
# This version of MCTS keeps track of the best found so far and reports
# that at the end
# really good tutorial https://www.youtube.com/watch?v=onBYsen2_eA
# if policy is incomplete; it does a random rollout for the remaining nodes
## clean up
startTime=clock()
evalRuns=0
tic()
## Define the Problem
#Define the design tree that we are searching
# We start from a seed of no color, this is indexed as 1 in matlab (0
# in graphsynth)

# There are two alternating choices, 1) Select a node to add to, and 2)
# Select the color of the node ROYGBV [indexed 1-6 for simplicity]

# How big should out sample off of every node be?
sample=100

# Create Space for recipe rows=node index, columns=[origin,color]
#different from ColorGraph.m which has 3 columns
rainbow=zeros(Nodes,3)
BESTrainbow=[]
BESTscore=[]
currentTime=clock()
TIMEELAPSED=currentTime-startTime

    ## SELECT - Select nodes as informed by the current policy
while TIMEELAPSED[5]*60+TIMEELAPSED[6]<=TIMEOUT*60
    # We are starting from the beginning so reset n to 1
    n=1
    CHECK=1
    S=1
    rainbow=zeros(Nodes,3)
    rainbow[:,1]=[1:Nodes]'

    while n<Nodes+1
        #Make a Random Graph
        run PureRandom.m
        #Count
        n=n+1
    end

    # Evaluate the generated rainbow graph
    #run ColorScore_mcts.m
    run ColorScore.m

    # Distance in 3d from target coordinate
    Qp=sqrt(GraphScore[1]^2+GraphScore[2]^2+GraphScore[3]^2)
    # Is it the best score?
    [BESTscore,BESTrainbow]=isbest(Qp,rainbow,BESTscore,BESTrainbow)

    # TRACK TIME PER RUN
    currentTime=clock()
    TIMEELAPSED=currentTime-startTime
    if toc()>60
        disp(strcat("Time Elapsed: ",...
            num2str(TIMEELAPSED[3]),"-days, ','_",...
            num2str(TIMEELAPSED[4]),"-hours, ','_",...
            num2str(TIMEELAPSED[5]),"-mins, ','_",...
            num2str(TIMEELAPSED[6]),"-secs, ','_"))
        tic()
    end
end
#Take best
rainbow=BESTrainbow
Qp=BESTscore

## Best score function()
function() [BESTscore,BESTrainbow]=isbest(Qp,rainbow,BESTscore,BESTrainbow)
# Check if the best score
if isempty(BESTrainbow)||BESTscore>Qp
    BESTscore=Qp
    BESTrainbow=rainbow
#     disp(strcat("BS: ',num2str(BESTscore),' - ',' Qp: ",num2str(Qp)))
#     disp(rainbow)
end
end
end

max_stars_repo_name: DesignEngrLab/ColorGraph

max_issues_repo_name: DesignEngrLab/ColorGraph

max_forks_repo_name: DesignEngrLab/ColorGraph


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Auxiliary functions to print the parts of the table.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# Print the table header without the horizontal lines.
function _print_table_header!(buf::IO,
                              display::Display,
                              header::Any,
                              header_str::Matrix{String},
                              header_len::Matrix{Int},
                              id_cols::Vector{Int},
                              id_rows::Vector{Int},
                              num_printed_cols::Int,
                              Δc::Int,
                              cols_width::Vector{Int},
                              vlines::Vector{Int},
                              # Configurations.
                              alignment::Vector{Symbol},
                              header_alignment::Vector{Symbol},
                              header_cell_alignment::Tuple,
                              show_row_names::Bool,
                              show_row_number::Bool,
                              tf::TextFormat,
                              # Crayons.
                              border_crayon::Crayon,
                              header_crayon::Vector{Crayon},
                              subheader_crayon::Vector{Crayon},
                              rownum_header_crayon::Crayon,
                              row_name_header_crayon::Crayon)

    header_num_rows, ~ = size(header_str)

    @inbounds @views for i = 1:header_num_rows
        0 ∈ vlines && _p!(display, border_crayon, tf.column, false, 1)

        for j = 1:num_printed_cols
            # Get the information about the alignment and the crayon.
            if j ≤ Δc
                if show_row_number && (j == 1)
                    crayon_ij    = rownum_header_crayon
                    alignment_ij = alignment[1]
                elseif show_row_names
                    crayon_ij    = row_name_header_crayon
                    alignment_ij = alignment[Δc]
                end
            else
                jc = id_cols[j-Δc]
                crayon_ij    = (i == 1) ? header_crayon[jc] : subheader_crayon[jc]
                alignment_ij = header_alignment[jc]

                # Check for cell alignment override.
                for f in header_cell_alignment
                    aux = f(header, i, jc)

                    if aux ∈ (:l, :c, :r, :L, :C, :R, :s, :S)
                        alignment_ij = aux
                        break
                    end
                end

                # If alignment is `:s`, then we must use the column
                # alignment.
                alignment_ij ∈ (:s,:S) && (alignment_ij = alignment[jc+Δc])
            end

            # Prepare the text to be printed.
            header_ij_str, header_ij_len = _str_aligned(header_str[i,j],
                                                        alignment_ij,
                                                        cols_width[j],
                                                        header_len[i,j])
            header_ij_str  = " " * header_ij_str * " "
            header_ij_len += 2

            flp = j == num_printed_cols

            # Print the text.
            _p!(display, crayon_ij, header_ij_str, false, header_ij_len)

            # Check if we need to draw a vertical line here.
            _pc!(j ∈ vlines, display, border_crayon, tf.column, "", flp, 1, 0)

            _eol(display) && break
        end

        i != header_num_rows && _nl!(display, buf)
    end

    _nl!(display, buf)
end

# Print the entire table data.
function _print_table_data(buf::IO,
                           display::Display,
                           data::Any,
                           data_str::Matrix{Vector{String}},
                           data_len::Matrix{Vector{Int}},
                           id_cols::Vector{Int},
                           id_rows::Vector{Int},
                           Δc::Int,
                           cols_width::Vector{Int},
                           vlines::Vector{Int},
                           row_printing_recipe::Vector{Union{Tuple{Int, Int, Int, Int},Symbol}},
                           col_printing_recipe::Vector{Union{Int,Symbol}},
                           # Configurations.
                           alignment::Vector{Symbol},
                           body_hlines_format::Tuple,
                           cell_alignment::Tuple,
                           continuation_row_alignment::Symbol,
                           ellipsis_line_skip::Integer,
                           highlighters::Tuple,
                           noheader::Bool,
                           show_row_names::Bool,
                           show_row_number::Bool,
                           tf::TextFormat,
                           # Crayons.
                           border_crayon::Crayon,
                           row_name_crayon::Crayon,
                           text_crayon::Crayon)

    line_count = 0

    for r in row_printing_recipe
        if r isa Symbol
            if r == :row_line
                _draw_line!(display, buf, body_hlines_format..., border_crayon,
                            cols_width, vlines)
            elseif r == :continuation_line
                _draw_continuation_row(display, buf, tf, text_crayon,
                                       border_crayon, cols_width, vlines,
                                       continuation_row_alignment)
            else
                error("Internal error: wrong symbol in row printing recipe.")
            end

        else
            i, ir, num_lines, l₀ = r

            for l = l₀:(l₀+num_lines-1)

                # Check if we should print the ellipsis here.
                display.cont_char =
                    line_count % (ellipsis_line_skip + 1) == 0 ? '⋯' : ' '
                line_count += 1

                num_col_recipes = length(col_printing_recipe)

                for c_id = 1:num_col_recipes
                    c = col_printing_recipe[c_id]

                    flp = c_id == num_col_recipes

                    if c isa Symbol
                        if c == :left_line
                            _p!(display, border_crayon, tf.column, flp, 1)
                        elseif c == :column_line
                            _p!(display, border_crayon, tf.column, flp, 1)
                        elseif c == :right_line
                            _p!(display, border_crayon, tf.column, flp, 1)
                        else
                            error("Internal error: wrong symbol in column printing recipe.")
                        end
                    else
                        j = c

                        # Get the information about the alignment and the
                        # crayon.
                        #
                        # NOTE: The variable `data_cell` is used to avoid
                        # applying formatters and highlighters to the row number
                        # and row name columns.
                        if j ≤ Δc
                            if show_row_number && (j == 1)
                                crayon_ij    = text_crayon
                                alignment_ij = alignment[1]
                            elseif show_row_names
                                crayon_ij    = row_name_crayon
                                alignment_ij = alignment[Δc]
                            end
                            jc = 0
                            data_cell = false
                        else
                            jc = id_cols[j-Δc]
                            crayon_ij    = text_crayon
                            alignment_ij = alignment[jc+Δc]
                            data_cell    = true
                        end

                        # String to be processed.
                        if length(data_str[i,j]) >= l
                            data_ij_str = data_str[i,j][l]
                            data_ij_len = data_len[i,j][l]
                        else
                            data_ij_str = ""
                            data_ij_len = 0
                        end

                        # Process the string with the correct alignment and also
                        # apply the highlighters.
                        data_ij_str, data_ij_len, crayon_ij =
                            _process_cell_text(data,
                                               ir,
                                               jc,
                                               data_cell,
                                               data_ij_str,
                                               data_ij_len,
                                               cols_width[j],
                                               crayon_ij,
                                               alignment_ij,
                                               cell_alignment,
                                               highlighters)

                        _p!(display, crayon_ij, data_ij_str, false, data_ij_len)
                    end
                end
                _nl!(display, buf)
            end
        end
    end
end


text:
#batched CuArray gemm by BatchedRoutines.jl
function batched_gemm!(transA::AbstractChar,
                       transB::AbstractChar,
                       alpha::Float32,
                       A::CUDA.CuArray{Float32, 3},
                       B::CUDA.CuArray{Float32, 3},
                       beta::Float32,
                       C::CUDA.CuArray{Float32, 3})
    CUDA.CUBLAS.gemm_strided_batched!(transA, transB, alpha, A, B, beta, C)
end

function batched_gemm!(transA::AbstractChar,
                       transB::AbstractChar,
                       alpha::Float64,
                       A::CUDA.CuArray{Float64, 3},
                       B::CUDA.CuArray{Float64, 3},
                       beta::Float64,
                       C::CUDA.CuArray{Float64, 3})
    CUDA.CUBLAS.gemm_strided_batched!(transA, transB, alpha, A, B, beta, C)
end


text:
"""
Creation and manipulation of ∇² differential operator
"""
module Laplacians

export Laplacian

files = (
    "1D/Laplacians1D.jl",
    "2D/Symmetric/Laplacians2D.jl",
    # "3D/Laplacians3D.jl"
    )

using ..Lattices
using ..Points
using SparseArrays

import ..Symmetric, ..Unsymmetric
import LinearAlgebra: I

struct Laplacian{N,CLASS}
    l0::SparseMatrixCSC{ComplexF64,Int}
    complex_scaling::SparseMatrixCSC{ComplexF64,Int}
end

foreach(include,files)

function Base.conj!(lap::Laplacian{N}) where N
    conj!(lap.l0)
    return lap
end
Base.conj(lap::Laplacian{N,CLASS}) where {N,CLASS} = Laplacian{N,CLASS}(conj(lap.l0),lap.complex_scaling)

################################################################################
# PRETTY PRINTING

import ..PRINTED_COLOR_DARK
import ..PRINTED_COLOR_VARIABLE

function Base.show(io::IO,::Laplacian{N,CLASS}) where {N,CLASS}
    print(io,"$(N)D ")
    CLASS<:Symmetric ? print(io,"Symmetric ") : nothing
    CLASS<:Unsymmetric ? print(io,"Unsymmetric ") : nothing
    printstyled(io,"Laplacian", color=PRINTED_COLOR_DARK)
end

end


text:
"""
    git()

Return a `Cmd` for running Git.

## Example

```julia
julia> run(`\$(git()) clone https://github.com/JuliaRegistries/General`)
```

This can equivalently be written with explicitly split arguments as

```
julia> run(git(["clone", "https://github.com/JuliaRegistries/General"]))
```

to bypass the parsing of the command string.
"""
function git(; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    @static if Sys.iswindows()
        return Git_jll.git(; adjust_PATH, adjust_LIBPATH)::Cmd
    else
        root = Git_jll.artifact_dir

        libexec = joinpath(root, "libexec")
        libexec_git_core = joinpath(libexec, "git-core")

        share = joinpath(root, "share")
        share_git_core = joinpath(share, "git-core")
        share_git_core_templates = joinpath(share_git_core, "templates")

        ssl_cert = joinpath(dirname(Sys.BINDIR), "share", "julia", "cert.pem")

        env_mapping = Dict{String,String}()
        env_mapping["GIT_EXEC_PATH"]    = libexec_git_core
        env_mapping["GIT_SSL_CAINFO"]   = ssl_cert
        env_mapping["GIT_TEMPLATE_DIR"] = share_git_core_templates

        original_cmd = Git_jll.git(; adjust_PATH, adjust_LIBPATH)::Cmd
        return addenv(original_cmd, env_mapping...)::Cmd
    end
end

function git(args::AbstractVector{<:AbstractString}; kwargs...)
    cmd = git(; kwargs...)
    append!(cmd.exec, args)
    return cmd
end


text:
using MATLAB
using Test

is_ci() = lowercase(get(ENV, "CI", "false")) == "true"

if !is_ci() # only test if not CI
include("engine.jl")
include("matfile.jl")
include("matstr.jl")
include("mxarray.jl")
end


text:
"""
    V(N)

Vocabulary size. (Tweedie & Baayen, 1998: p. 326)
"""
V(N::Int) = unique() # TODO

"""
    V(i,N)

Frequency spectrum
"""
V(i::Int,N::Int) = unique() # TODO

# function frequency_spectrum(tokens)
#     fl = Dict{String,Int64}()
#     fs = Dict{Int64,Int64}()
#     for token in tokens
#         fl[token] = get(fl, token, 0) + 1
#     end
#     for f in values(fl)
#         fs[f] = get(fs, f, 0) + 1
#     end
#     fs
# end

"""
    P(N)

Vocabulary growth rate. (Tweedie & Baayen, 1998: p. 325)
"""
P(N::Int) = V(1, N)/N

"""
    MWF(N)

Mean word frequency (Eq. 2)

```math
MWF(N) = \\frac{N}{V(N)}
```
"""
MWF(N) = N/V(N)

"""
    TTR(N)

Type token ratio (Eq. 3)

```math
TTR(N) = \\frac{V(N)}{N}
```
"""
ttr(tokens::Array{T}) where T = length(unique(tokens))/length(tokens)

"""
    guiraud_r(N)

Guiraud (1954) (T&B Eq. 4)

```math
r = \\frac{V(N)}{\\sqrt{N}}
```
"""
guiraud_r(N) = V(N)/sqrt(N)

"""
    herdan_c(N,c)

Herdan (1960; 1964) (T&B Eq. 5)

```math
c = 0.0
```
"""
herdan_c(N,c) = log(N^c)/log(N)

"""
    dugast_k(N)

Dugast (1979) (T&B Eq. 6)

```math
k = \\frac{\\log V(N)}{\\log(\\log N)}
```
"""
dugast_k(N) = log(log(k,N))/log(log(N))

"""
    a²(N)

Maas (1972)

```math
V(N) = N \\log^{-a^2} N
```
"""
a²(N) = (log(N)-log(V(N)))/log(2,N)

"""
    dugast_U(N)

Dugast U (1978, 1979) (Eq. 8)

```math
V(N) = N \\log^{-1/U} N
```
"""
dugast_U(N) = 1/a²(N)

"""
    LN(N)

Lukjanenkov & Nesitoj (1975; in Tulvada, 1977) (Eq. 9)

```math
V(N) = \\frac{1}{\\sqrt{1 + LN \\log N}}
```
"""
LN(N) = (1-V(N)^2)/(V(N)^2*log(N))

"""
    brunet_W(N)

Brunet's parametric expression W (1978) (Eq. 10)

```math
W = N^{V(N)^{-a}},
```
where ```a = -0.172``` and
```math
V(N) = (\\frac{\\log W}{\\log N})^a = \\log^a_N W.
```
"""
brunet_W(N) = N^(V(N)^(-0.172))


text:
@testset "sort.jl" begin
    @testset "sort" begin
        a = rand(10)
        test_rrule(sort, a)
        test_rrule(sort, a; fkwargs=(;rev=true))
    end
    @testset "partialsort" begin
        a = rand(10)
        test_rrule(partialsort, a, 4)
        test_rrule(partialsort, a, 3:5 ⊢ DoesNotExist())
        test_rrule(partialsort, a, 1:2:6 ⊢ DoesNotExist())

        test_rrule(partialsort, a, 4, fkwargs=(;rev=true))
    end
end


text:
@testset "Safety" begin
  include("SimpleBoreCollisionModule.jl")
end


text:
"""
This function runs the actual energy and ancillary services market clearing module using endogeneous ED.
Returns the market clearing prices, capacity factors and reserve percentages.
"""
function energy_mkt_clearing(sys_UC::Nothing,
                             sys_ED::Nothing,
                             sys_local_ED::MarketClearingProblem,
                             load_growth::AxisArrays.AxisArray{Float64, 1},
                             zones::Vector{String},
                             num_days::Int64,
                             solver::JuMP.MOI.OptimizerWithAttributes)
    energy_price,
    reserve_up_price,
    reserve_down_price,
    capacity_factors,
    reserve_up_perc,
    reserve_down_perc = economicdispatch(sys_local_ED,
                                       solver,
                                       "C:/Users/manwar2/Documents/GitRepos/emt-tests/data/simulation_data/results/ed_results.txt")

    return energy_price, reserve_up_price, reserve_down_price, capacity_factors, reserve_up_perc, reserve_down_perc;
end

"""
This function runs the actual energy and ancillary services market clearing module using SIIP PSI Simualtion.
Returns the market clearing prices, capacity factors and reserve percentages.
"""
function energy_mkt_clearing(sys_UC::PSY.System,
                             sys_ED::PSY.System,
                             sys_local_ED::Union{Nothing, MarketClearingProblem},
                             load_growth::AxisArrays.AxisArray{Float64, 1},
                             zones::Vector{String},
                             num_days::Int64,
                             solver::JuMP.MOI.OptimizerWithAttributes)

    update_PSY_timeseries!(sys_UC, load_growth, num_days)

    energy_price,
    reserve_up_price,
    reserve_down_price,
    capacity_factors,
    reserve_up_perc,
    reserve_down_perc = create_simulation(sys_UC, sys_ED, zones, num_days, solver)

    return energy_price, reserve_up_price, reserve_down_price, capacity_factors, reserve_up_perc, reserve_down_perc;
end


text:
fname = :TRSpectral

c = Combi(hessian_dense,PDataSpectral,solve_modelTRDiag,preprocessSpectral,decreaseFact,Tparam())

include("Template.jl")


text:
module MD5

using SHA
using SHA: lrot, SHA_CTX
import SHA: update!, digest!, transform!, pad_remainder!, buffer_pointer
import SHA: blocklen, state_type, digestlen, short_blocklen
export md5

include("constants.jl")
include("types.jl")
include("core.jl")

# Our basic function is to process arrays of bytes
function md5(data::T) where T<:Union{AbstractVector{UInt8}, NTuple{N,UInt8} where N}
    ctx = MD5_CTX()
    update!(ctx, data)
    return digest!(ctx)
end

# AbstractStrings are a pretty handy thing to be able to crunch through
md5(str::AbstractString) = md5(codeunits(str))

# Convenience function for IO devices, allows for things like:
# open("test.txt") do f
#     sha256(f)
# done
function md5(io::IO, chunk_size=4*1024)
    ctx = MD5_CTX()
    buff = Vector{UInt8}(undef, chunk_size)
    while !eof(io)
        num_read = readbytes!(io, buff)
        update!(ctx, buff[1:num_read])
    end
    return digest!(ctx)
end

end # module


text:
using Images, Base.Test, Color, FixedPointNumbers

macro chk(a, b)
    :(@test $a == $b && typeof($a) == typeof($b))
end

# MapNone
mapi = MapNone{Int}()
@chk map(mapi, 7) 7
@chk map(mapi, 0x07) 7
a7 = Int[7]
@chk map(mapi, [0x07]) a7
@test map(mapi, a7) === a7
mapi = MapNone{RGB24}()
g = Ufixed8(0.1)
@chk Images.map1(mapi, 0.1) g
@chk map(mapi, 0.1) RGB24(g,g,g)
@chk map(mapi, Gray(0.1)) RGB24(g,g,g)
@chk map(mapi, g) RGB24(g,g,g)
@chk map(mapi, true) RGB24(0xff,0xff,0xff)
@chk map(mapi, false) RGB24(0x00,0x00,0x00)
mapi = MapNone{RGB{Float32}}()
g = 0.1f0
@chk Images.map1(mapi, 0.1) g
@chk map(mapi, 0.1) RGB(g,g,g)

c = RGB(0.1,0.2,0.3)
mapi = MapNone{HSV{Float64}}()
@chk map(mapi, c) convert(HSV, c)

# issue #200
c = RGBA{Ufixed8}(1,0.5,0.25,0.8)
mapi = MapNone{Images.ColorTypes.BGRA{Ufixed8}}()
@chk map(mapi, c) convert(Images.ColorTypes.BGRA{Ufixed8}, c)

# BitShift
mapi = BitShift{Uint8,7}()
@chk map(mapi, 0xff) 0x01
@chk map(mapi, 0xf0ff) 0xff
@chk map(mapi, [0xff]) Uint8[0x01]
mapi = BitShift{Ufixed8,7}()
@chk map(mapi, 0xffuf8) 0x01uf8
@chk map(mapi, 0xf0ffuf16) 0xffuf8
mapi = BitShift{ARGB32,4}()
@chk map(mapi, 0xffuf8) ARGB32(0xff0f0f0f)
mapi = BitShift{RGB24,2}()
@chk map(mapi, Gray(0xffuf8)) RGB24(0x003f3f3f)
mapi = BitShift{ARGB32,2}()
@chk map(mapi, Gray(0xffuf8)) ARGB32(0xff3f3f3f)
@chk map(mapi, GrayAlpha{Ufixed8}(Gray(0xffuf8),0x3fuf8)) ARGB32(0x0f3f3f3f)
mapi = BitShift{RGB{Ufixed8},2}()
@chk map(mapi, Gray(0xffuf8)) RGB(0x3fuf8, 0x3fuf8, 0x3fuf8)
mapi = BitShift{ARGB{Ufixed8},2}()
@chk map(mapi, Gray(0xffuf8)) ARGB{Ufixed8}(0x3fuf8, 0x3fuf8, 0x3fuf8, 0xffuf8)
@chk map(mapi, GrayAlpha{Ufixed8}(Gray(0xffuf8),0x3fuf8)) ARGB{Ufixed8}(0x3fuf8, 0x3fuf8, 0x3fuf8, 0x0fuf8)
mapi = BitShift{RGBA{Ufixed8},2}()
@chk map(mapi, Gray(0xffuf8)) RGBA{Ufixed8}(0x3fuf8, 0x3fuf8, 0x3fuf8, 0xffuf8)
@chk map(mapi, GrayAlpha{Ufixed8}(Gray(0xffuf8),0x3fuf8)) RGBA{Ufixed8}(0x3fuf8, 0x3fuf8, 0x3fuf8, 0x0fuf8)
mapi = BitShift(ARGB{Ufixed8}, 8)
@chk map(mapi, RGB{Ufixed16}(1.0,0.8,0.6)) ARGB{Ufixed8}(1.0,0.8,0.6,1.0)
mapi = BitShift(RGBA{Ufixed8}, 8)
@chk map(mapi, RGB{Ufixed16}(1.0,0.8,0.6)) RGBA{Ufixed8}(1.0,0.8,0.6,1.0)
# Issue, #269, IJulia issue #294
bs = BitShift(Gray{Ufixed8}, 8)
v = Gray(ufixed16(0.8))
@chk map(bs, v) Gray{Ufixed8}(0.8)

# Clamp
mapi = ClampMin(Float32, 0.0)
@chk map(mapi,  1.2) 1.2f0
@chk map(mapi, -1.2) 0.0f0
mapi = ClampMin(RGB24, 0.0f0)
@chk map(mapi, RGB{Float32}(-5.3,0.4,0.8)) RGB24(0x000066cc)
mapi = ClampMax(Float32, 1.0)
@chk map(mapi,  1.2)  1.0f0
@chk map(mapi, -1.2) -1.2f0
mapi = ClampMax(RGB24, 1.0f0)
@chk map(mapi, RGB{Float32}(0.2,1.3,0.8)) RGB24(0x0033ffcc)
mapi = ClampMinMax(Float32, 0.0, 1.0)
@chk map(mapi,  1.2) 1.0f0
@chk map(mapi, -1.2) 0.0f0
mapi = ClampMinMax(ARGB32, 0.0f0, 1.0f0)
@chk map(mapi, RGB{Float32}(-0.2,1.3,0.8)) ARGB32(0xff00ffcc)
@chk map(mapi, ARGB{Float32}(-0.2,1.3,0.8,0.6)) ARGB32(0x9900ffcc)
mapi = Clamp(Float32)
@chk map(mapi,  1.2) 1.0f0
@chk map(mapi, -1.2) 0.0f0
mapi = Clamp(Ufixed12)
@chk map(mapi, Ufixed12(1.2)) one(Ufixed12)
mapi = Clamp(Gray{Ufixed12})
@chk map(mapi, Gray(Ufixed12(1.2))) Gray(one(Ufixed12))
mapi = ClampMinMax(RGB24, 0.0, 1.0)
@chk map(mapi, 1.2) RGB24(0x00ffffff)
@chk map(mapi, 0.5) RGB24(0x00808080)
@chk map(mapi, -.3) RGB24(0x00000000)
mapi = ClampMinMax(RGB{Ufixed8}, 0.0, 1.0)
@chk map(mapi, 1.2) RGB{Ufixed8}(1,1,1)
@chk map(mapi, 0.5) RGB{Ufixed8}(0.5,0.5,0.5)
@chk map(mapi, -.3) RGB{Ufixed8}(0,0,0)
mapi = Clamp(RGB{Ufixed8})
@chk map(mapi, RGB(1.2,0.5,-.3)) RGB{Ufixed8}(1,0.5,0)
mapi = Clamp(ARGB{Ufixed8})
@chk map(mapi, ARGB{Float64}(1.2,0.5,-.3,0.2)) ARGB{Ufixed8}(1.0,0.5,0.0,0.2)
@chk map(mapi, RGBA{Float64}(1.2,0.5,-.3,0.2)) ARGB{Ufixed8}(1.0,0.5,0.0,0.2)
@chk map(mapi, 0.2) ARGB{Ufixed8}(0.2,0.2,0.2,1.0)
@chk map(mapi, GrayAlpha{Float32}(Gray(0.2),1.2)) ARGB{Ufixed8}(0.2,0.2,0.2,1.0)
@chk map(mapi, GrayAlpha{Float32}(Gray(-.4),0.8)) ARGB{Ufixed8}(0.0,0.0,0.0,0.8)
mapi = Clamp(RGBA{Ufixed8})
@chk map(mapi, ARGB{Float64}(1.2,0.5,-.3,0.2)) RGBA{Ufixed8}(1.0,0.5,0.0,0.2)
@chk map(mapi, RGBA{Float64}(1.2,0.5,-.3,0.2)) RGBA{Ufixed8}(1.0,0.5,0.0,0.2)
@chk map(mapi, 0.2) RGBA{Ufixed8}(0.2,0.2,0.2,1.0)
@chk map(mapi, GrayAlpha{Float32}(Gray(0.2),1.2)) RGBA{Ufixed8}(0.2,0.2,0.2,1.0)
@chk map(mapi, GrayAlpha{Float32}(Gray(-.4),0.8)) RGBA{Ufixed8}(0.0,0.0,0.0,0.8)
# Issue #253
mapi = Clamp(BGRA{Ufixed8})
@chk map(mapi, RGBA{Float32}(1.2,0.5,-.3,0.2)) BGRA{Ufixed8}(1.0,0.5,0.0,0.2)

@chk clamp(RGB{Float32}(-0.2,0.5,1.8)) RGB{Float32}(0.0,0.5,1.0)
@chk clamp(ARGB{Float64}(1.2,0.5,-.3,0.2)) ARGB{Float64}(1.0,0.5,0.0,0.2)
@chk clamp(RGBA{Float64}(1.2,0.5,-.3,0.2)) RGBA{Float64}(1.0,0.5,0.0,0.2)

# ScaleMinMax
mapi = ScaleMinMax(Ufixed8, 100, 1000)
@chk map(mapi, 100) Ufixed8(0.0)
@chk map(mapi, 1000) Ufixed8(1.0)
@chk map(mapi, 10) Ufixed8(0.0)
@chk map(mapi, 2000) Ufixed8(1.0)
@chk map(mapi, 550) Ufixed8(0.5)
mapinew = ScaleMinMax(Ufixed8, [100,500,1000])
@test mapinew == mapi
mapinew = ScaleMinMax(Ufixed8, [0,500,2000], convert(Uint16, 100), convert(Uint16, 1000))
@test mapinew == mapi
mapi = ScaleMinMax(ARGB32, 100, 1000)
@chk map(mapi, 100) ARGB32(0x00,0x00,0x00,0xff)
@chk map(mapi, 550) ARGB32(0x80,0x80,0x80,0xff)
@chk map(mapi,2000) ARGB32(0xff,0xff,0xff,0xff)
mapi = ScaleMinMax(RGB{Float32}, 100, 1000)
@chk map(mapi,  50) RGB(0.0f0, 0.0f0, 0.0f0)
@chk map(mapi, 550) RGB{Float32}(0.5, 0.5, 0.5)
@chk map(mapi,2000) RGB(1.0f0, 1.0f0, 1.0f0)
A = Gray{Ufixed8}[Ufixed8(0.1), Ufixed8(0.9)]
@test mapinfo(RGB24, A) == MapNone{RGB24}()
mapi = ScaleMinMax(RGB24, A, zero(Gray{Ufixed8}), one(Gray{Ufixed8}))
@test map(mapi, A) == map(mapinfo(RGB24, A), A)
mapi = ScaleMinMax(Float32, [Gray(one(Ufixed8))], 0, 1) # issue #180
@chk map(mapi, Gray(Ufixed8(0.6))) 0.6f0
@test_throws ErrorException ScaleMinMax(Float32, 0, 0, 1.0) # issue #245
A = [Gray{Float64}(0.2)]
mapi = ScaleMinMax(RGB{Ufixed8}, A, 0.0, 0.2)
@test map(mapi, A) == [RGB{Ufixed8}(1,1,1)]

# ScaleSigned
mapi = ScaleSigned(Float32, 1/5)
@chk map(mapi, 7) 1.0f0
@chk map(mapi, 5) 1.0f0
@chk map(mapi, 3) convert(Float32, 3/5)
@chk map(mapi, -3) convert(Float32, -3/5)
@chk map(mapi, -6) -1.0f0
mapi = ScaleSigned(RGB24, 1.0f0/10)
@chk map(mapi, 12) RGB24(0x00ff00ff)
@chk map(mapi, -10.0) RGB24(0x0000ff00)
@chk map(mapi, 0) RGB24(0x00000000)

# ScaleAutoMinMax
mapi = ScaleAutoMinMax()
A = [100,550,1000]
@chk map(mapi, A) ufixed8([0.0,0.5,1.0])
mapi = ScaleAutoMinMax(RGB24)
@chk map(mapi, A) RGB24[0x00000000, 0x00808080, 0x00ffffff]

# scaling, ssd
img = Images.grayim(fill(typemax(Uint16), 3, 3))
mapi = Images.mapinfo(Ufixed8, img)
img8 = map(mapi, img)
@assert all(img8 .== typemax(Ufixed8))
A = 0
mnA, mxA = 1.0, -1.0
while mnA > 0 || mxA < 0
    A = randn(3,3)
    mnA, mxA = extrema(A)
end
offset = 30.0
img = convert(Images.Image, A .+ offset)
mapi = Images.ScaleMinMax(Ufixed8, offset, offset+mxA, 1/mxA)
imgs = map(mapi, img)
@assert minimum(imgs) == 0
@assert maximum(imgs) == 1
@assert eltype(imgs) == Ufixed8
imgs = Images.imadjustintensity(img, [])
@test_throws ErrorException Images.imadjustintensity(img, [1])
mnA = minimum(A)
@assert Images.ssd(imgs, (A.-mnA)/(mxA-mnA)) < eps()
A = reshape(1:9, 3, 3)
B = map(Images.ClampMin(Float32, 3), A)
@assert eltype(B) == Float32 && B == [3 4 7; 3 5 8; 3 6 9]
B = map(Images.ClampMax(Uint8, 7), A)
@assert eltype(B) == Uint8 && B == [1 4 7; 2 5 7; 3 6 7]

A = reinterpret(Ufixed8, [convert(UInt8,1):convert(UInt8,24);], (3, 2, 4))
img = reinterpret(RGB{Ufixed8}, A, (2,4))
@test separate(img) == permutedims(A, (2,3,1))

# sc
arr = zeros(4,4)
arr[2,2] = 0.5
@assert sc(arr)[2,2] == 0xffuf8
@assert sc(arr, 0.0, 0.75)[2,2] == 0xaauf8

# color conversion
gray = collect(linspace(0.0,1.0,5)) # a 1-dimensional image
gray8 = round(Uint8, 255*gray)
gray32 = Uint32[convert(UInt32, g)<<16 | convert(UInt32, g)<<8 | convert(UInt32, g) for g in gray8]
imgray = Images.Image(gray, Dict{ASCIIString,Any}([("colordim",0), ("colorspace","Gray")]))
buf = map(Images.mapinfo(Uint32, imgray), imgray) # Images.uint32color(imgray)
@test buf == gray32
rgb = RGB{Float64}[RGB(g, g, g) for g in gray]
buf = map(Images.mapinfo(Uint32, rgb), rgb) # Images.uint32color(rgb)
@test buf == gray32
r = red(rgb)
@test r == gray
img = Images.Image(reinterpret(RGB24, gray32)) # , ["colordim"=>0, "colorspace"=>"RGB24"])
buf = map(Images.mapinfo(Uint32, img), img) # Images.uint32color(img)
@test buf == gray32
rgb = repeat(gray, outer=[1,3])
img = Images.Image(rgb, Dict{ASCIIString,Any}([("colordim",2), ("colorspace","RGB"), ("spatialorder",["x"])]))
buf = map(Images.mapinfo(Uint32, img), img) # Images.uint32color(img)
@test buf == gray32
g = green(img)
@test g == gray
rgb = repeat(gray', outer=[3,1])
img = Images.Image(rgb, Dict{ASCIIString,Any}([("colordim",1), ("colorspace","RGB"), ("spatialorder",["x"])]))
buf = map(Images.mapinfo(Uint32, img), img) # Images.uint32color(img)
@test buf == gray32
b = blue(img)
@test b == gray

# map and indexed images
img = Images.ImageCmap([1 2 3; 3 2 1], [RGB{Ufixed16}(1.0,0.6,0.4), RGB{Ufixed16}(0.2, 0.4, 0.6), RGB{Ufixed16}(0.5,0.5,0.5)])
mapi = MapNone(RGB{Ufixed8})
imgd = map(mapi, img)
cmap = [RGB{Ufixed8}(1.0,0.6,0.4), RGB{Ufixed8}(0.2, 0.4, 0.6), RGB{Ufixed8}(0.5,0.5,0.5)]
@test imgd == reshape(cmap[[1,3,2,2,3,1]], (2,3))


text:
using Pkg, Documenter

try
    using JAXTAM
catch ArgumentError
    Pkg.activate(joinpath(pwd(), ".."))
    using JAXTAM
end

@info "Current dir: $(pwd())"

DOCUMENTER_DEBUG=true

makedocs(
    modules  = [JAXTAM],
    doctest  = false,
    clean    = false,
    sitename = "JAXTAM Documentation",
    format   = :html
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    repo  = "github.com/RobertRosca/JAXTAM.jl.git",
    julia = "0.7",
    deps  = nothing,
    make  = nothing,
)


text:
using RvSpectraKitLearn
using Base.Test

# write your own tests here
@test 1 == 1

include(joinpath(Pkg.dir("RvSpectraKitLearn"),"path_to_spectra.jl"))
cd(path_to_spectra)

(datafiles, phases) = get_filenames();
(lambda, obs) = read_filelist(datafiles);

# Test 1
using MultivariateStats
@test test_pca_eford(obs) < 1e-8

include("test2.jl")

max_stars_repo_name: eford/RVSpectraKitLearn

max_issues_repo_name: eford/RVSpectraKitLearn

max_forks_repo_name: eford/RVSpectraKitLearn


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".

function ismeuv(wave::T, hcol::T, he1col::T, he2col::T, fano::Bool) where {T<:AbstractFloat}
    r = wave/911.75
    if r > 1
        return zero(T)
    end
    # minexp = -708.39642  from alog((machar(double=double).xmin)) where double = 1b in IDL
    z = sqrt(r/(1 - r))
    denom = -expm1(-2*T(π)*z)
    tauH = ((3.44e-16)*(r^4)*exp(-4*z*atan(1/z))*hcol)/denom
    r *= 4
    tauHe2 = zero(T)
    if r < 1
        z = sqrt(r/(1 - r))
        denom = -4*expm1(-2*T(π)*z)
        tauHe2 = ((3.44e-16)*(r^4)*exp(-4*z*atan(1/z))*he2col)/denom
    end
    q  = (2.81, 2.51, 2.45, 2.44)
    fano_gamma = (2.64061e-03, 6.20116e-04, 2.56061e-04, 1.320159e-04)
    esubi = (4.421529414644497, 4.679309217126802, 4.738680412951545, 4.764345016358231)
    tauHe1 = zero(T)
    if wave < 503.97
        x = log10(wave)
        if wave < 46
            y = @evalpoly(x, -2.465188e+01, 4.354679, -3.553024, 5.573040,
                          -5.872938, 3.720797, -1.226919, 1.576657e-01)
        else
            y = @evalpoly(x, -2.953607e+01, 7.083061, 8.678646e-01, -1.221932,
                          4.052997e-02, 1.317109e-01, -3.265795e-02,
                          2.500933e-03)
            if fano
                episilon = 911.2671/wave
                for i = 1:4
                    x = 2*(episilon - esubi[i])/fano_gamma[i]
                    y += log10(((x - q[i])^2)/(1 + x^2))
                end
            end
        end
        tauHe1 = exp10(y)*he1col
    end
    return tauH + tauHe1 + tauHe2
end

"""
    ismeuv(wave, hcol[, he1col=hcol*0.1, he2col=0, fano=false]) -> tau

### Purpose ###

Compute the continuum interstellar EUV optical depth

### Explanation ###

The EUV optical depth is computed from the photoionization of hydrogen and helium.

### Arguments ###

* `wave`: wavelength value (in Angstroms). Useful range is 40 - 912 A;
  at shorter wavelength metal opacity should be considered, at longer wavelengths
  there is no photoionization.
* `hcol`: interstellar hydrogen column density in cm-2.
* `he1col` (optional): neutral helium column density in cm-2.
  Default is 0.1*hcol (10% of hydrogen column)
* `he2col` (optional): ionized helium column density in cm-2
  Default is 0.
* `fano` (optional boolean keyword): If this keyword is true, then the 4 strongest
  auto-ionizing resonances of He I are included. The shape of these resonances
  is given by a Fano profile - see Rumph, Bowyer, & Vennes 1994, AJ, 107, 2108.
  If these resonances are included then the input wavelength vector should have
  a fine (>~0.01 A) grid between 190 A and 210 A, since the resonances are very narrow.

### Output ###

* `tau`: Vector giving resulting optical depth, non-negative values.

### Example ###

One has a model EUV spectrum with wavelength, w (in Angstroms).
Find the EUV optical depth by 1e18 cm-2 of HI, with N(HeI)/N(HI) = N(HeII)/N(HI) = 0.05.

```jldoctest
julia> using AstroLib

julia> ismeuv.([670, 910], 1e19, 5e17, 5e17)
2-element Array{Float64,1}:
 27.35393320556168
 62.683796028917286
```

### Notes ###

Code of this function is based on IDL Astronomy User's Library.
"""
ismeuv(wave::Real, hcol::Real, he1col::Real=hcol/10, he2col::Real=0,
       fano::Bool=false) =
           ismeuv(promote(float(wave), float(hcol), float(he1col),
                          float(he2col))..., fano)


text:
using LinearAlgebra
function loss(w, s=4.0)
	return (w[1]*w[2]*w[3] - s)^2/2
end
function grad_loss(w, s=4.0)
	term = w[1]*w[2]*w[3] - s
	return term*[w[2]*w[3],
				w[1]*w[3],
				w[1]*w[2]]
end
function hess_loss(w, s=4.0)
	term = w[1]*w[2]*w[3] - s
	dterm = [w[2]*w[3] w[1]*w[3] w[1]*w[2]]

	return term*[0 w[3] w[2];
				 w[3] 0 w[1];
				 w[2] w[1] 0] .+
	[w[2]*w[3]*dterm; w[1]*w[3]*dterm;
	w[1]*w[2]*dterm]
end
next(w, η, s=4.0) = w - η*grad_loss(w, s)
dnext(w, η, s=4.0) = I(3) - η*hess_loss(w, s)
function d2next(w, η, s=4.0)
	term = w[1]*w[2]*w[3] - s
	dterm = [w[2]*w[3] w[1]*w[3] w[1]*w[2]]
	d2term = [0 w[3] w[2];w[3] 0 w[1];
			 w[2] w[1] 0]
    dfirst = term*[0 0 0;0 0 1;0 1 0;
				   0 0 1;0 0 0;1 0 0;
				   0 1 0;1 0 0;0 0 0] .+
	[0 0 0;w[3]*dterm;w[2]*dterm;
	 w[3]*dterm;0 0 0;w[1]*dterm;
	w[2]*dterm;w[1]*dterm;0 0 0]

	dsecond = [0 2*w[2]*w[3]*w[3] 2*w[2]*w[2]*w[3];
			   w[2]*w[3]*w[3] w[1]*w[3]*w[3] 2*w[1]*w[2]*w[3];
			   w[2]*w[2]*w[3] 2*w[1]*w[2]*w[3] w[1]*w[2]*w[2];
			   w[2]*w[3]*w[3] w[1]*w[3]*w[3] 2*w[1]*w[2]*w[3];
			   2*w[3]*w[3]*w[1] 0 2*w[3]*w[1]*w[1];
			   2*w[1]*w[2]*w[3] w[1]*w[1]*w[3] w[1]*w[1]*w[2];
			   w[2]*w[2]*w[3] 2*w[2]*w[1]*w[3] w[1]*w[2]*w[2];
			   2*w[1]*w[2]*w[3] w[1]*w[1]*w[3] w[1]*w[1]*w[2];
			   2*w[1]*w[2]*w[2] 2*w[1]*w[1]*w[2] 0]

    return -η*(dfirst .+ dsecond)
end

max_stars_repo_name: nishaChandramoorthy/neuralOMET

max_issues_repo_name: nishaChandramoorthy/neuralOMET

max_forks_repo_name: nishaChandramoorthy/neuralOMET


text:
using HierarchicalGeometry
using Test

@testset "HierarchicalGeometry.jl" begin
    # Write your tests here.
end


text:
# Spatial calculations

export
    vector1D,
    vector,
    adjust_bounds

"""
    vector1D(c1, c2, box_size)

Displacement between two 1D coordinate values, accounting for the bounding box.
The minimum image convention is used, so the displacement is to the closest
version of the coordinate accounting for the periodic boundaries.
"""
function vector1D(c1::Real, c2::Real, box_size::Real)
    if c1 < c2
        return (c2 - c1) < (c1 - c2 + box_size) ? (c2 - c1) : (c2 - c1 - box_size)
    else
        return (c1 - c2) < (c2 - c1 + box_size) ? (c2 - c1) : (c2 - c1 + box_size)
    end
end

"""
    vector(c1, c2, box_size)

Displacement between two coordinate values, accounting for the bounding box.
The minimum image convention is used, so the displacement is to the closest
version of the coordinates accounting for the periodic boundaries.
"""
vector(c1, c2, box_size::Real) = vector1D.(c1, c2, box_size)

"""
    adjust_bounds(c, box_size)

Ensure a coordinate is within the simulation box and return the coordinate.
"""
function adjust_bounds(c::Real, box_size::Real)
    while c >= box_size
        c -= box_size
    end
    while c < zero(c)
        c += box_size
    end
    return c
end


text:
immutable Arcsine <: ContinuousUnivariateDistribution
    a::Float64
    b::Float64

    function Arcsine(a::Float64, b::Float64)
        a < b || error("a must be less than b.")
        new(a, b)
    end

    Arcsine(a::Real, b::Real) = Arcsine(float64(a), float64(b))
    Arcsine(b::Real) = Arcsine(0.0, float64(b))
    Arcsine() = new(0.0, 1.0)
end

@distr_support Arcsine d.a d.b

### Parameters

params(d::Arcsine) = (d.a, d.b)
location(d::Arcsine) = d.a
scale(d::Arcsine) = d.b - d.a

### Statistics

mean(d::Arcsine) = (d.a + d.b) * 0.5
median(d::Arcsine) = mean(d)
mode(d::Arcsine) = d.a
modes(d::Arcsine) = [d.a, d.b]

var(d::Arcsine) = 0.125 * abs2(d.b - d.a)
skewness(d::Arcsine) = 0.0
kurtosis(d::Arcsine) = -1.5

entropy(d::Arcsine) = -0.24156447527049044469 + log(scale(d))

### Evaluation

pdf(d::Arcsine, x::Float64) = insupport(d, x) ? 1.0 / (π * sqrt((x - d.a) * (d.b - x))) : 0.0

logpdf(d::Arcsine, x::Float64) = insupport(d, x) ? -(logπ + 0.5 * log((x - d.a) * (d.b - x))) : -Inf

cdf(d::Arcsine, x::Float64) = x < d.a ? 0.0 :
                              x > d.b ? 1.0 :
                              0.636619772367581343 * asin(sqrt((x - d.a) / (d.b - d.a)))

quantile(d::Arcsine, p::Float64) = location(d) + abs2(sin(halfπ * p)) * scale(d)

### Sampling

rand(d::Arcsine) = quantile(d, rand())


text:
@testset "approximators" begin
    @testset "TabularQ" begin
        A = TabularQ([1 2; 3 4])
        @test A(1) == [1., 2.]
        @test A(2) == [3., 4.]
        @test A(2, 2) == 4
        @test A(2,Val(:dist)) == [3., 4.]
        @test A(2,Val(:max)) == 4.
        @test A(2,Val(:argmax)) == 2

        update!(A, 2, 2, -1.)
        @test A(2) == [3., 3.]
        @test A(2, 2) == 3

        update!(A, 1, [-1., -1.])
        @test A(1) == [0., 1.]
        @test A(1, 1) == 0
        @test A(1, 2) == 1
    end

    @testset "TabularV" begin
        V = TabularV([1., 2., 3., 2., 1.])
        @test V(1) == 1.
        @test V(3) == 3.

        update!(V, 3, -1.)
        @test V(3) == 2.
    end

    @testset "AggregationV" begin
        table = [-1., 1.]
        V = AggregationV(table, s -> s < 0 ? 1 : 2)

        @test V(-5) == -1.
        @test V(5) == 1.

        update!(V, -1, 0.5)
        update!(V, 1, -0.5)

        @test V(-5) == -.5
        @test V(5) == 0.5
    end

    @testset "FourierV" begin
        V = FourierV([1, 1, 1])
        @test V(1) ≈ 1.0
        @test V(2) ≈ 3.0

        update!(V, 1, 0.5)
        @test V(1) ≈ 2.5
        @test V(2) ≈ 3.5
    end

    @testset "LinearV" begin
        features = [1 0 1; 0 1 1]
        weights = [1., 2., 3.]
        V = LinearV(features, weights)

        @test V(1) ≈ 1. + 3.
        @test V(2) ≈ 2. + 3.

        update!(V, 1, 0.5)
        @test V(1) ≈ 1.5 + 3.5
        @test V(2) ≈ 2. + 3.5
    end

    @testset "PolynomialV" begin
        V = PolynomialV([1., 2., 3.])

        @test V(2) ≈ 1 + 2. * 2  + 3. * 2^2

        update!(V, 2, 1.)
        @test V(2) ≈ 2 + 4. * 2 + 7. * 2^2
    end

    @testset "tilings" begin
        init_tiling = Tiling((0:5:50, 0:10:100))
        tilings = [init_tiling, init_tiling - [1, 2]]
        @testset "TilingsV" begin
            V = TilingsV(tilings)
            update!(V, [0, 0], 1.)
            @test V([0,0]) == 2.
            @test V([4,8]) == 1.
            @test V([5,10]) == 0.
        end
        @testset "TilingsQ" begin
            Q = TilingsQ(tilings, 2)
            update!(Q, [0, 0], 1, 1.)

            @test Q([0, 0]) ≈ [2., 0.]
            @test Q([0, 0], 1) ≈ 2.
            @test Q([0, 0], Val(:max)) ≈ 2.
            @test Q([0, 0], Val(:argmax)) == 1
        end
    end
end


text:
#!/usr/bin/env julia
# A wrapper script to launch one or several GeMM simulation runs, parallelising
# them where possible. *Warning:* Because this script includes its own commandline
# parsing, no commandline args are passed on to GeMM!

using # TODO: automatic package installation - ask user?
    Distributed,
    ArgParse,
    Random

@everywhere thisDir = joinpath(pwd(), "src")
@everywhere any(path -> path == thisDir, LOAD_PATH) || push!(LOAD_PATH, thisDir)
@everywhere using Pkg
@everywhere Pkg.activate(".")
@everywhere using GeMM

"""
    parseexperiments()

Get configuration file(s), random seed and number of replicates from the commandline.
"""
function parseexperiments()
    s = ArgParseSettings()
    @add_arg_table! s begin
        "--config", "-c"
            help = "list of configuration files, comma separated"
            arg_type = String
            required = true
        "--seed", "-s"
            help = "initial random seed"
            arg_type = Int
            required = false
            default = 0
        "--replicates", "-n"
            help = "number of replicates to run"
            arg_type = Int
            range_tester = x -> x > 0
            required = false
            default = 1
    end
    args = parse_args(s)
    for a in keys(args)
        (args[a] == nothing) && delete!(args, a)
    end
    args
end

args = parseexperiments()
configs = split(args["config"], ',')

simulations = []

for i in args["seed"]:(args["seed"] + args["replicates"] - 1)
    args["seed"] == 0 ? rep = abs(rand(RandomDevice(), Int32)) : rep = i
    for config in configs
        push!(simulations, [String(config), rep])
    end
end

# Only use `pmap` when we actually have multiple simulations.
# Otherwise, we can save on the overhead (plus, `pmap` obscures
# error messages).
if length(simulations) == 1
    sim = simulations[1]
    rungemm(sim[1], sim[2])
else
    pmap(x -> rungemm(x[1], x[2]), simulations)
end

max_stars_repo_name: CCTB-Ecomods/gemm

max_issues_repo_name: CCTB-Ecomods/gemm

max_forks_repo_name: CCTB-Ecomods/gemm


text:
# test basic manifold product behaviour

using ApproxManifoldProducts
using KernelDensityEstimate
using Test

# # 2D
#
# pts1 = [0.05*randn(100)'; 0.5*randn(100)']
# pts2 = [0.05*randn(100)'; 0.5*randn(100)']
#
# P1 = manikde!(pts1, (:Euclid, :Circular))
# P2 = manikde!(pts2, (:Euclid, :Circular))
#
# P12 = P1*P2
#
#
#
# plotKDE([P1;P2;P12], c=["red";"blue";"magenta"],levels=1)
#

# 3D

@testset "Basic 3D multiply and cross dimension covariance test..." begin

N = 100

pts1 = [0.05*randn(N)'; 0.05*randn(N)'; 0.75*randn(N)']
pts2 = [0.05*randn(N)'; 0.05*randn(N)'; 0.75*randn(N)']

# P1 = manikde!(pts1, (:Euclid,:Euclid,:Circular))
# P2 = manikde!(pts2, (:Euclid,:Euclid,:Circular))

# P1 = manikde!(pts1, (:Euclid,:Euclid,:Euclid))
# P2 = manikde!(pts2, (:Euclid,:Euclid,:Euclid))

P2 = kde!(pts2)
P1 = kde!(pts1)

P12 = P1*P2

pts = getPoints(P12)

@test 0.8*N < sum(abs.(pts[1,:]) .< 0.1)
@test 0.8*N < sum(abs.(pts[2,:]) .< 0.1)
@test 0.8*N < sum(abs.(pts[3,:]) .< 2.0)

P1 = manikde!(pts1, (:Euclid,:Euclid,:Euclid))
P2 = manikde!(pts2, (:Euclid,:Euclid,:Euclid))

P12 = P1*P2

pts = getPoints(P12)

@test 0.8*N < sum(abs.(pts[1,:]) .< 0.1)
@test 0.8*N < sum(abs.(pts[2,:]) .< 0.1)
@test 0.8*N < sum(abs.(pts[3,:]) .< 2.0)

P1 = manikde!(pts1, (:Euclid,:Euclid,:Circular))
P2 = manikde!(pts2, (:Euclid,:Euclid,:Circular))

P12 = P1*P2

pts = getPoints(P12)

@test 0.8*N < sum(abs.(pts[1,:]) .< 0.1)
@test 0.8*N < sum(abs.(pts[2,:]) .< 0.1)
@test 0.8*N < sum(abs.(pts[3,:]) .< 2.0)

# using KernelDensityEstimatePlotting
# using Gadfly
# Gadfly.set_default_plot_size(35cm,25cm)
#
# plotKDE([P1;P2;P12], c=["red";"blue";"magenta"],levels=1) |> PDF("/tmp/test.pdf")
# run(`evince /tmp/test.pdf`)

end

#


text:
# Properties of toric varieties
# Properties of toric varieties

@doc Markdown.doc"""
    is_isomorphic_to_projective_space( v::AbstractNormalToricVariety )

Checks if the normal toric variety `v` is isomorphic to projective space.

# Examples
```jldoctest
julia> is_isomorphic_to_projective_space( projective_space( 2 ) )
true
```
"""
function is_isomorphic_to_projective_space( v::AbstractNormalToricVariety )
    return GAP.Globals.IsIsomorphicToProjectiveSpace( v.GapNTV )::Bool
end
export is_isomorphic_to_projective_space

@doc Markdown.doc"""
    is_direct_product_of_projective_spaces( v::AbstractNormalToricVariety )

Checks if the normal toric variety `v` is isomorphic to a direct product of projective space.

# Examples
```jldoctest
julia> is_direct_product_of_projective_spaces( projective_space( 2 ) )
true
```
"""
function is_direct_product_of_projective_spaces( v::AbstractNormalToricVariety )
    return GAP.Globals.IsDirectProductOfPNs( v.GapNTV )::Bool
end
export is_direct_product_of_projective_spaces

# Attributes of toric varieties
# Attributes of toric varieties

@doc Markdown.doc"""
    affine_open_covering( v::AbstractNormalToricVariety )

Computes an affine open cover of the normal toric variety `v`, i.e. returns a list of affine toric varieties.
"""
function affine_open_covering( v::AbstractNormalToricVariety )
    gap_cover = GAP.Globals.AffineOpenCovering( v.GapNTV )
    return [ NormalToricVariety( v ) for v in gap_cover ]
end
export affine_open_covering

struct CoxRing
           GapCoxRing::GapObj
end
export CoxRing

@doc Markdown.doc"""
    cox_ring( v::AbstractNormalToricVariety )

Computes the Cox ring of the normal toric variety `v`.
"""
function cox_ring( v::AbstractNormalToricVariety )
    gap_ring = GAP.Globals.CoxRing( v.GapNTV )
    return CoxRing( gap_ring )
end
export cox_ring

@doc Markdown.doc"""
    list_of_variables_of_cox_ring( v::AbstractNormalToricVariety )

Computes the list of homogeneous variables of the Cox ring of the normal toric variety `v`.
"""
function list_of_variables_of_cox_ring(v::AbstractNormalToricVariety)
    vars = GAP.Globals.ListOfVariablesOfCoxRing(v.GapNTV)
    return Vector{String}(vars)
end
export list_of_variables_of_cox_ring

struct ClassGroup
           GapClassGroup::GapObj
end
export ClassGroup

@doc Markdown.doc"""
    class_group( v::AbstractNormalToricVariety )

Computes the class group of the normal toric variety `v`.
"""
function class_group( v::AbstractNormalToricVariety )
    gap_class_group = GAP.Globals.ClassGroup( v.GapNTV )
    return ClassGroup( gap_class_group )
end
export class_group

struct TorusInvariantDivisorGroup
           GapTorusInvariantDivisorGroup::GapObj
end
export TorusInvariantDivisorGroup

@doc Markdown.doc"""
    torus_invariant_divisor_group( v::AbstractNormalToricVariety )

Computes the torus invariant divisor class group of the normal toric variety `v`.
"""
function torus_invariant_divisor_group( v::AbstractNormalToricVariety )
    gap_TorusInvariantDivisorGroup = GAP.Globals.TorusInvariantDivisorGroup( v.GapNTV )
    return TorusInvariantDivisorGroup( gap_TorusInvariantDivisorGroup )
end
export torus_invariant_divisor_group

struct MapFromCharacterToPrincipalDivisor
           GapMapFromCharacterToPrincipalDivisor::GapObj
end
export MapFromCharacterToPrincipalDivisor

@doc Markdown.doc"""
    map_from_character_to_principal_divisor( v::AbstractNormalToricVariety )

Computes the map from the character lattice to the principal divisors of the normal toric variety `v`.
"""
function map_from_character_to_principal_divisor( v::AbstractNormalToricVariety )
    gap_MapFromCharacterToPrincipalDivisor = GAP.Globals.MapFromCharacterToPrincipalDivisor( v.GapNTV )
    return MapFromCharacterToPrincipalDivisor( gap_MapFromCharacterToPrincipalDivisor )
end
export map_from_character_to_principal_divisor

struct MapFromWeilDivisorsToClassGroup
           GapMapFromWeilDivisorsToClassGroup::GapObj
end
export MapFromWeilDivisorsToClassGroup

@doc Markdown.doc"""
    map_from_weil_divisors_to_class_group( v::AbstractNormalToricVariety )

Computes the map from the Weil divisors to the Class group of the normal toric variety `v`.
"""
function map_from_weil_divisors_to_class_group( v::AbstractNormalToricVariety )
    gap_MapFromWeilDivisorsToClassGroup = GAP.Globals.MapFromWeilDivisorsToClassGroup( v.GapNTV )
    return MapFromWeilDivisorsToClassGroup( gap_MapFromWeilDivisorsToClassGroup )
end
export map_from_weil_divisors_to_class_group

struct CoordinateRingOfTorus
           GapCoordinateRingOfTorus::GapObj
end
export CoordinateRingOfTorus

@doc Markdown.doc"""
    coordinate_ring_of_torus( v::AbstractNormalToricVariety )

Computes the coordinate ring of the torus of the normal toric variety `v`.
"""
function coordinate_ring_of_torus( v::AbstractNormalToricVariety )
    gap_CoordinateRingOfTorus = GAP.Globals.CoordinateRingOfTorus( v.GapNTV )
    return CoordinateRingOfTorus( gap_CoordinateRingOfTorus )
end
export coordinate_ring_of_torus

@doc Markdown.doc"""
    list_of_variables_of_coordinate_ring_of_torus( v::AbstractNormalToricVariety )

Computes the list of homogeneous coordinates of the coordinate ring of the torus of the normal toric variety `v`.
"""
function list_of_variables_of_coordinate_ring_of_torus(v::AbstractNormalToricVariety)
    vars = GAP.Globals.ListOfVariablesOfCoordinateRingOfTorus(v.GapNTV)
    return Vector{String}(vars)
end
export list_of_variables_of_coordinate_ring_of_torus

@doc Markdown.doc"""
    is_product_of( v::AbstractNormalToricVariety )

Identifies the factors from which the normal toric variety `v` has been constructed in GAP.
"""
function is_product_of(v::AbstractNormalToricVariety)
    factors = GAP.Globals.IsProductOf(v.GapNTV)
    return [ NormalToricVariety( f ) for f in factors ]
end
export is_product_of

@doc Markdown.doc"""
    factors( v::AbstractNormalToricVariety )

Identifies the factors from which the normal toric variety `v` has been constructed in GAP.
"""
function factors( v::AbstractNormalToricVariety )
    gap_factors = GAP.Globals.Factors( v.GapNTV )
    return [ NormalToricVariety( f ) for f in gap_factors ]
end
export factors

struct CharacterLattice
           GapCharacterLattice::GapObj
end
export CharacterLattice

@doc Markdown.doc"""
    character_lattice( v::AbstractNormalToricVariety )

Computes the character lattice of the normal toric variety `v`.
"""
function character_lattice( v::AbstractNormalToricVariety )
    gap_CharacterLattice = GAP.Globals.CharacterLattice( v.GapNTV )
    return CharacterLattice( gap_CharacterLattice )
end
export character_lattice

@doc Markdown.doc"""
    torus_invariant_prime_divisors( v::AbstractNormalToricVariety )

Computes the torus invariant prime divisors of the normal toric variety `v`.
"""
function torus_invariant_prime_divisors( v::AbstractNormalToricVariety )
    divisors = GAP.Globals.TorusInvariantPrimeDivisors( v.GapNTV )
    return [ ToricDivisor(extract_gap_divisor_coeffs(d), v) for d in divisors ]
end
export torus_invariant_prime_divisors

struct IrrelevantIdeal
           GapIrrelevantIdeal::GapObj
end
export IrrelevantIdeal

@doc Markdown.doc"""
    irrelevant_ideal( v::AbstractNormalToricVariety )

Computes the irrelevant ideal of the normal toric variety `v`.
"""
function irrelevant_ideal( v::AbstractNormalToricVariety )
    gap_IrrelevantIdeal = GAP.Globals.IrrelevantIdeal( v.GapNTV )
    return IrrelevantIdeal( gap_IrrelevantIdeal )
end
export irrelevant_ideal

struct StanleyReisnerIdeal
           GapSRIdeal::GapObj
end
export StanleyReisnerIdeal

@doc Markdown.doc"""
    stanley_reisner_ideal( v::AbstractNormalToricVariety )

Computes the Stanley-Reisner ideal of the normal toric variety `v`.
"""
function stanley_reisner_ideal( v::AbstractNormalToricVariety )
    gap_SRIdeal = GAP.Globals.SRIdeal( v.GapNTV )
    return StanleyReisnerIdeal( gap_SRIdeal )
end
export stanley_reisner_ideal

struct MorphismFromCoxVariety
           GapMorphismFromCoxVariety::GapObj
end
export MorphismFromCoxVariety

@doc Markdown.doc"""
    morphism_from_cox_variety( v::AbstractNormalToricVariety )

Computes the morphism from the Cox variety of the normal toric variety `v`.
"""
function morphism_from_cox_variety( v::AbstractNormalToricVariety )
    gap_MorphismFromCoxVariety = GAP.Globals.MorphismFromCoxVariety( v.GapNTV )
    return MorphismFromCoxVariety( gap_MorphismFromCoxVariety )
end
export morphism_from_cox_variety

@doc Markdown.doc"""
    cox_variety( v::AbstractNormalToricVariety )

Computes the Cox variety of the normal toric variety `v`.
"""
function cox_variety( v::AbstractNormalToricVariety )
    gap_CoxVariety = GAP.Globals.CoxVariety( v.GapNTV )
    return NormalToricVariety( gap_CoxVariety )
end
export cox_variety

struct Fan
           gap_fan::GapObj
           rays::Vector{Vector{Int}}
           cones::Vector{Vector{Int}}
end
export Fan

@doc Markdown.doc"""
    fan_of_variety( v::AbstractNormalToricVariety )

Computes the fan of the normal toric variety `v`.
"""
function fan_of_variety( v::AbstractNormalToricVariety )
    # collect data
    gap_fan = GAP.Globals.FanOfVariety( v.GapNTV )
    rays = Vector{Vector{Int}}( GAP.Globals.RayGenerators( gap_fan ) )
    cones = Vector{Vector{Int}}( GAP.Globals.RaysInMaximalCones( gap_fan ) )
    cones = [ findall( x -> x == 1, c ) for c in cones ]

    # return the fan
    return Fan( gap_fan, rays, cones )
end
export fan_of_variety

@doc Markdown.doc"""
    fan( v::AbstractNormalToricVariety )

A convenience method for `fan_of_variety( v::AbstractNormalToricVariety )`.
"""
function fan( v::AbstractNormalToricVariety )
    return fan_of_variety( v::AbstractNormalToricVariety )
end
export fan

struct CartierTorusInvariantDivisorGroup
           GapCartierTorusInvariantDivisorGroup::GapObj
end
export CartierTorusInvariantDivisorGroup

@doc Markdown.doc"""
    cartier_torus_invariant_divisor_group( v::AbstractNormalToricVariety )

Computes the group of Cartier and torus invariant divisors of the normal toric variety `v`.
"""
function cartier_torus_invariant_divisor_group( v::AbstractNormalToricVariety )
    gap_CartierTorusInvariantDivisorGroup = GAP.Globals.CartierTorusInvariantDivisorGroup( v.GapNTV )
    return CartierTorusInvariantDivisorGroup( gap_CartierTorusInvariantDivisorGroup )
end
export cartier_torus_invariant_divisor_group

struct PicardGroup
           GapPicardGroup::GapObj
end
export PicardGroup

@doc Markdown.doc"""
    picard_group( v::AbstractNormalToricVariety )

Computes the Picard group of the normal toric variety `v`.
"""
function picard_group( v::AbstractNormalToricVariety )
    gap_PicardGroup = GAP.Globals.PicardGroup( v.GapNTV )
    return PicardGroup( gap_PicardGroup )
end
export picard_group

@doc Markdown.doc"""
    name_of_variety( v::AbstractNormalToricVariety )

Returns the name of the normal toric variety `v`, if set. Otherwise returns "No name set for this variety".
"""
function name_of_variety( v::AbstractNormalToricVariety )
    if ! GAP.Globals.HasNameOfVariety( v.GapNTV )
            return "No name set for this variety"
    end

    return String( GAP.Globals.NameOfVariety( v.GapNTV ) )
end
export name_of_variety

struct ZariskiCotangentSheaf
           GapZariskiCotangentSheaf::GapObj
end
export ZariskiCotangentSheaf

@doc Markdown.doc"""
    zariski_cotangent_sheaf( v::AbstractNormalToricVariety )

Returns the Zariski cotangent sheaf of the normal toric variety `v`.
"""
function zariski_cotangent_sheaf( v::AbstractNormalToricVariety )
    gap_ZariskiCotangentSheaf = GAP.Globals.ZariskiCotangentSheaf( v.GapNTV )
    return ZariskiCotangentSheaf( gap_ZariskiCotangentSheaf )
end
export zariski_cotangent_sheaf

struct CotangentSheaf
           GapCotangentSheaf::GapObj
end
export CotangentSheaf

@doc Markdown.doc"""
    cotangent_sheaf( v::AbstractNormalToricVariety )

Returns the cotangent sheaf of the normal toric variety `v`.
"""
function cotangent_sheaf( v::AbstractNormalToricVariety )
    gap_CotangentSheaf = GAP.Globals.CotangentSheaf( v.GapNTV )
    return CotangentSheaf( gap_CotangentSheaf )
end
export cotangent_sheaf

@doc Markdown.doc"""
    weil_divisors_of_variety( v::AbstractNormalToricVariety )

Compute the Weil divisors of the normal toric variety `v`.
"""
function weil_divisors_of_variety( v::AbstractNormalToricVariety )
    gap_divisors = GAP.Globals.WeilDivisorsOfVariety( v.GapNTV )
    return [ToricDivisor(extract_gap_divisor_coeffs(d), v) for d in gap_divisors ]
end
export weil_divisors_of_variety

struct ZariskiCotangentSheafViaEulerSequence
           GapZariskiCotangentSheafViaEulerSequence::GapObj
end
export ZariskiCotangentSheafViaEulerSequence

@doc Markdown.doc"""
    zariski_cotangent_sheaf_via_euler_sequence( v::AbstractNormalToricVariety )

Computes the Zariski cotangent sheaf of the normal toric variety `v` via the Euler sequence.
"""
function zariski_cotangent_sheaf_via_euler_sequence( v::AbstractNormalToricVariety )
    gap_ZariskiCotangentSheafViaEulerSequence = GAP.Globals.ZariskiCotangentSheafViaEulerSequence( v.GapNTV )
    return ZariskiCotangentSheafViaEulerSequence( gap_ZariskiCotangentSheafViaEulerSequence )
end
export zariski_cotangent_sheaf_via_euler_sequence

struct ZariskiCotangentSheafViaPoincareResidueMap
           GapZariskiCotangentSheafViaPoincareResidueMap::GapObj
end
export ZariskiCotangentSheafViaPoincareResidueMap

@doc Markdown.doc"""
    zariski_cotangent_sheaf_via_poincare_residue_map( v::AbstractNormalToricVariety )

Computes the Zariski cotangent sheaf of the normal toric variety `v` via the Poincare residue map.
"""
function zariski_cotangent_sheaf_via_poincare_residue_map( v::AbstractNormalToricVariety )
    gap_ZariskiCotangentSheafViaPoincareResidueMap = GAP.Globals.ZariskiCotangentSheafViaPoincareResidueMap( v.GapNTV )
    return ZariskiCotangentSheafViaPoincareResidueMap( gap_ZariskiCotangentSheafViaPoincareResidueMap )
end
export zariski_cotangent_sheaf_via_poincare_residue_map

#struct UnderlyingSheaf
#          GapUnderlyingSheaf::GapObj
#end
#export UnderlyingSheaf

#function underlying_sheaf( v::AbstractNormalToricVariety )
#   gap_Underlying = GAP.Globals.UnderlyingSheaf( v.GapNTV )
#    return UnderlyingSheaf( gap_UnderlyingSheaf )
#end
#export underlying_sheaf

# Methods of toric varieties
# Methods of toric varieties

@doc Markdown.doc"""
    set_name_of_variety( v::AbstractNormalToricVariety, name::String )

Sets the name of the normal toric variety `v` to `name`.
"""
function set_name_of_variety( v::AbstractNormalToricVariety, s::String )
    GAP.Globals.SetNameOfVariety( v.GapNTV, GapObj( s ) )
    return true
end
export set_name_of_variety

@doc Markdown.doc"""
    coordinate_ring_of_torus( v::AbstractNormalToricVariety, names::Vector{String} )

Compute the coordinate ring of the torus factor of the normal toric variety `v`, using `names` as label for the homogeneous coordinates.
"""
function coordinate_ring_of_torus( v::AbstractNormalToricVariety, names::Vector{String} )
    gap_names = [ GapObj( names[ i ] ) for i in 1 : size(names)[1] ]
    gap_names = GapObj( gap_names )
    gap_CoordinateRingOfTorus = GAP.Globals.CoordinateRingOfTorus( v.GapNTV, gap_names )
    return CoordinateRingOfTorus( gap_CoordinateRingOfTorus )
end
export coordinate_ring_of_torus

@doc Markdown.doc"""
    cox_ring( v::AbstractNormalToricVariety, name::String )

Compute the Cox ring of the normal toric variety `v`, using `name` as label for the homogeneous coordinates.
"""
function cox_ring( v::AbstractNormalToricVariety, names::String )
    gap_names = GapObj( names )
    gap_CoxRing = GAP.Globals.CoxRing( v.GapNTV, gap_names )
    return CoxRing( gap_CoxRing )
end
export cox_ring

@doc Markdown.doc"""
    Base.:*( v::AbstractNormalToricVariety, w::AbstractNormalToricVariety )

Compute the direct product of the normal toric varieties `v` and `w`.

# Examples
```jldoctest
julia> projective_space( 2 ) * projective_space( 2 )
A normal toric variety corresponding to a polyhedral fan in ambient dimension 4
```
"""
function Base.:*( v::AbstractNormalToricVariety, w::AbstractNormalToricVariety )
    gap_NormalToricVariety = v.GapNTV * w.GapNTV
    return NormalToricVariety( gap_NormalToricVariety )
end
export *

struct CharacterToRationalFunction
           gap_CharacterToRationalFunction::GapObj
end
export CharacterToRationalFunction

@doc Markdown.doc"""
    character_to_rational_function( l::Vector{Int}, v::AbstractNormalToricVariety )

Turn the character `l` of the normal toric variety `v` into a rational function.
"""
function character_to_rational_function( l::Vector{Int}, v::AbstractNormalToricVariety )
    gap_CharacterToRationalFunction = GAP.Globals.CharacterToRationalFunction( GapObj( l ), v.GapNTV )
    return CharacterToRationalFunction( gap_CharacterToRationalFunction )
end
export character_to_rational_function

@doc Markdown.doc"""
    blowup_on_ith_minimal_torus_orbit( v::AbstractNormalToricVariety, i::Int )

Compute the blow up of the normal toric variety `v` on the i-th minimal torus orbit.
"""
function blowup_on_ith_minimal_torus_orbit( v::AbstractNormalToricVariety, i::Int )
    gap_blowup_variety = GAP.Globals.BlowUpOnIthMinimalTorusOrbit( v.GapNTV, GapObj( i ) )
    return NormalToricVariety( gap_blowup_variety )
end
export blowup_on_ith_minimal_torus_orbit

@doc Markdown.doc"""
    nr_of_q_rational_points( v::AbstractNormalToricVariety, i::Int )

Compute the number of q-rational points of the normal toric variety `v`.
"""
function nr_of_q_rational_points( v::AbstractNormalToricVariety, i::Int )
    return GAP.Globals.NrOfqRationalPoints( v.GapNTV, GapObj( i ) )::Int
end
export nr_of_q_rational_points

# ToricDivisors
# ToricDivisors

@doc Markdown.doc"""
    divisor_of_character( c::Vector{Int}, v::AbstractNormalToricVariety )

Construct the torus invariant divisor on the normal toric variety `v` corresponding to the character `c`.

# Examples
```jldoctest
julia> divisor_of_character( [1,2], projective_space( 2 ) )
A torus invariant divisor on a normal toric variety
```
"""
function divisor_of_character( character::Vector{Int}, v::AbstractNormalToricVariety )
    # create the divisor
    gap_character = GapObj( character )
    gap_divisor = GAP.Globals.DivisorOfCharacter( gap_character, v.GapNTV )

    # wrap and return
    return ToricDivisor(extract_gap_divisor_coeffs(gap_divisor), v)
end
export divisor_of_character

@doc Markdown.doc"""
    divisor_of_class(class::AbstractVector, v::AbstractNormalToricVariety)

Construct a torus invariant divisor on the normal toric variety `v` corresponding to the divisor class `c`.

# Examples
```jldoctest
julia> divisor_of_class( [1], projective_space( 2 ) )
A torus invariant divisor on a normal toric variety
```
"""
function divisor_of_class(class::AbstractVector, v::AbstractNormalToricVariety)
    # create the divisor
    gap_class = GapObj( class )
    gap_divisor = GAP.Globals.DivisorOfGivenClass( v.GapNTV, gap_class )

    # wrap and return
    return ToricDivisor(extract_gap_divisor_coeffs(gap_divisor), v)
end
export divisor_of_class

@doc Markdown.doc"""
    is_primedivisor( d::ToricDivisor )

Checks if the divisor `d` is prime.
"""
function is_primedivisor( d::ToricDivisor )
    return GAP.Globals.IsPrimedivisor( d.GapToricDivisor )::Bool
end
export is_primedivisor

# Disabled tests
# Disabled tests

@testset "Hirzebruch surface" begin
    @test JToric.is_isomorphic_to_projective_space( H5 ) == false
    @test JToric.is_direct_product_of_projective_spaces( H5 ) == false
    @test JToric.is_isomorphic_to_projective_space( P2 ) == true
    @test JToric.is_direct_product_of_projective_spaces( P2 ) == true
    @test JToric.is_isomorphic_to_projective_space( blowup_variety ) == false
    @test JToric.is_direct_product_of_projective_spaces( blowup_variety ) == false
    @test JToric.is_isomorphic_to_projective_space( v ) == false
    @test JToric.is_direct_product_of_projective_spaces( v ) == false
    @test size( affine_open_covering( H5 ) )[ 1 ] == 4
    cox_ring( H5 )
    cox_ring( H5, "u" )
    coordinate_ring_of_torus( H5, [ "u", "v", "w", "z" ] )
    @test size( list_of_variables_of_cox_ring( H5 ) )[ 1 ] == 4
    class_group( H5 )
    torus_invariant_divisor_group( H5 )
    map_from_character_to_principal_divisor( H5 )
    map_from_weil_divisors_to_class_group( H5 )
    coordinate_ring_of_torus( H5 )
    @test size( list_of_variables_of_coordinate_ring_of_torus( H5 ) )[ 1 ] == 4
    @test size( is_product_of( H5 ) )[ 1 ] == 1
    character_lattice( H5 )
    divisors = torus_invariant_prime_divisors( H5 )
    @test ( 0 in [ is_primedivisor( d ) for d in divisors ] ) == false
    irrelevant_ideal( H5 )
    stanley_reisner_ideal( H5 )
    morphism_from_cox_variety( H5 )
    cox_variety( H5 )
    fan_of_variety( H5 )
    fan( H5 )
    cartier_torus_invariant_divisor_group( H5 )
    picard_group( H5 )
    @test name_of_variety( H5 ) == "No name set for this variety"
    set_name_of_variety( H5, "Hirzebruch surface" )
    @test name_of_variety( H5 ) == "Hirzebruch surface"
    zariski_cotangent_sheaf( H5 )
    cotangent_sheaf( H5 )
    # UnderlyingSheaf( H5 ) <- Error in gap
    character_to_rational_function( [1,2,3,4], H5 )
    # @test size( weil_divisors_of_variety( H5 ) )[ 1 ] == 4
    @test size( factors( H5 ) )[ 1 ] == 1
    zariski_cotangent_sheaf_via_euler_sequence( H5 )
    zariski_cotangent_sheaf_via_poincare_residue_map( H5 )
    @test nr_of_q_rational_points( H5, 1 ) == 4
end

@testset "Projective space" begin
    @test nr_of_q_rational_points( P2, 2 ) == 4
end

@testset "Blowup of projective space" begin
    # Perform tests for blowup of projective space
    blowup_variety = blowup_on_ith_minimal_torus_orbit( P2, 1 )
    @test JToric.isnormal( blowup_variety ) == true
    @test JToric.isaffine( blowup_variety ) == false
    @test JToric.isprojective( blowup_variety ) == true
    @test JToric.issmooth( blowup_variety ) == true
    @test JToric.iscomplete( blowup_variety ) == true
    @test JToric.has_torusfactor( blowup_variety ) == false
    @test JToric.is_orbifold( blowup_variety ) == true
    @test JToric.issimplicial( blowup_variety ) == true
    @test ith_betti_number( blowup_variety, 0 ) == 1
    @test ith_betti_number( blowup_variety, 1 ) == 0
    @test ith_betti_number( blowup_variety, 2 ) == 2
    @test ith_betti_number( blowup_variety, 3 ) == 0
    @test ith_betti_number( blowup_variety, 4 ) == 1
    @test nr_of_q_rational_points( blowup_variety, 4 ) == 13
end

@testset "Direct products" begin
    # Perform tests for direct product of Hirzebruch surface and projective space
    v = H5 * P2;
    @test JToric.isnormal( v ) == true
    @test JToric.isaffine( v ) == false
    @test JToric.isprojective( v ) == true
    @test JToric.issmooth( v ) == true
    @test JToric.iscomplete( v ) == true
    @test JToric.has_torusfactor( v ) == false
    @test JToric.is_orbifold( v ) == true
    @test JToric.issimplicial( v ) == true
    @test size( factors( v ) )[ 1 ] == 2
    @test ith_betti_number( v, 0 ) == 1
    @test ith_betti_number( v, 1 ) == 0
    @test ith_betti_number( v, 2 ) == 3
    @test ith_betti_number( v, 3 ) == 0
    @test ith_betti_number( v, 4 ) == 4
    @test ith_betti_number( v, 5 ) == 0
    @test ith_betti_number( v, 6 ) == 3
    @test ith_betti_number( v, 7 ) == 0
    @test ith_betti_number( v, 8 ) == 1
    @test nr_of_q_rational_points( v, 3 ) == 106
end

@testset "Divisors" begin
    @test JToric.is_primedivisor( D ) == false
    D2 = divisor_of_character( [ 1,2 ], H5 )
    @test JToric.iscartier( D2 ) == true
    @test JToric.isprincipal( D2 ) == true
    @test JToric.isbasepoint_free( D2 ) == true
    @test JToric.isample( D2 ) == false
    @test JToric.isvery_ample( D2 ) == false
    @test JToric.isnef( D2 ) == true
    D3 = divisor_of_class( [ 1,2 ], H5 )
    @test JToric.iscartier( D3 ) == true
    @test JToric.isprincipal( D3 ) == false
    @test JToric.isbasepoint_free( D3 ) == true
    @test JToric.isample( D3 ) == true
    @test JToric.isvery_ample( D3 ) == true
    @test JToric.isnef( D3 ) == true
    @test JToric.is_primedivisor( D2 ) == false
    @test JToric.is_primedivisor( D3 ) == false
end


text:
#####
##### Resnet 50 from Metalhead
#####
struct ResidualBlock
  conv_layers
  norm_layers
  shortcut
end

Flux.@treelike ResidualBlock

function ResidualBlock(
        filters,
        kernels::Array{Tuple{Int,Int}},
        pads::Array{Tuple{Int,Int}},
        strides::Array{Tuple{Int,Int}},
        shortcut = identity
    )
    conv_layers = [
        Conv(kernels[i-1], filters[i-1] => filters[i], pad = pads[i-1], stride = strides[i-1])
        for i in 2:length(filters)
    ]

    ResidualBlock(Tuple(conv_layers), Tuple(norm_layers), shortcut)
end

function ResidualBlock(
        filters,
        kernels::Array{Int},
        pads::Array{Int},
        strides::Array{Int},
        shortcut = identity
    )
    ResidualBlock(
        filters,
        [(i,i) for i in kernels],
        [(i,i) for i in pads],
        [(i,i) for i in strides],
        shortcut
    )
end

function (block::ResidualBlock)(input)
    value = input
    for i in 1:length(block.conv_layers)-1
        value = relu.(block.norm_layers[i]((block.conv_layers[i])(value)))
    end
    relu.(((block.norm_layers[end])((block.conv_layers[end])(value))) + block.shortcut(input))
end

function Bottleneck(filters::Int, downsample::Bool = false, res_top::Bool = false)
    if(!downsample && !res_top)
        return ResidualBlock(
            [4 * filters, filters, filters, 4 * filters],
            [1,3,1],
            [0,1,0],
            [1,1,1]
        )
    elseif(downsample && res_top)
        return ResidualBlock(
            [filters, filters, filters, 4 * filters],
            [1,3,1],
            [0,1,0],
            [1,1,1],
            Chain(
                Conv(
                    (1,1),
                    filters=>4 * filters,
                    pad = (0,0),
                    stride = (1,1)
                ),
                BatchNorm(4 * filters)
           )
        )
    else
        shortcut = Chain(
            Conv((1,1), 2 * filters=>4 * filters, pad = (0,0), stride = (2,2)),
            BatchNorm(4 * filters)
        )
        return ResidualBlock(
            [2 * filters, filters, filters, 4 * filters],
            [1,3,1],
            [0,1,0],
            [1,1,2],
            shortcut
        )
    end
end

function resnet50()
    local layers = [3, 4, 6, 3]
    local layer_arr = []

    push!(layer_arr, Conv((7,7), 3=>64, pad = (3,3), stride = (2,2)))
    push!(layer_arr, x -> maxpool(x, (3,3), pad = (1,1), stride = (2,2)))

    initial_filters = 64
    for i in 1:length(layers)
        push!(layer_arr, Bottleneck(initial_filters, true, i==1))
        for j in 2:layers[i]
            push!(layer_arr, Bottleneck(initial_filters))
        end
        initial_filters *= 2
    end

    push!(layer_arr, x -> meanpool(x, (7,7)))
    push!(layer_arr, x -> reshape(x, :, size(x,4)))
    push!(layer_arr, (Dense(2048, 1000)))
    push!(layer_arr, softmax)

    #Chain(layer_arr...)
    function f(x)
        for l in layer_arr
            x = l(x)
        end
        return x
    end
end

@testset "Testing Resnet" begin
    f = resnet50()
    x = rand(Float32, 224, 224, 3, 16)
    z = f(x)

    backend = nGraph.Backend()
    X = nGraph.Tensor(backend, x)
    ex = nGraph.compile(backend, f, X)

    @test isapprox(read(ex(X)), z)
    @time ex(X)
    @time ex(X)

    y = similar(z)
    loss(x,y) = Flux.crossentropy(f(x), y)
    Y = nGraph.Tensor(backend, y)

    G = nGraph.compile(backend, loss, X, Y)
    G(X,Y)
end


text:
#module Cabida_nlo_evol

##############################################
# PARTE "1": CARGA CABIDA                    #
##############################################

using LotMassing, .poly2D, .polyShape, Random, CSV

#Random.seed!(1236)
#Random.seed!(1230)

##############################################
# PARTE "2": GENERACIÓN DE PARÁMETROS        #
##############################################

idPredio = 8 #8 predio = 1,2,3,4,5,6,75,8
conjuntoTemplates = [4] #4 [1:L, 2:C, 4:V]

fpe = FlagPlotEdif3D(true,  # predio
                     true,  # volTeor
                     true,  # restSombra
                     true,  # edif
                     true,  # sombraVolTeor_p
                     true,  # sombraVolTeor_o
                     true,  # sombraVolTeor_s
                     true,  # sombraEdif_p
                     true,  # sombraEdif_o
                     true)  # sombraEdif_s

dcn = datosCabidaNormativa(5, # SEPMIN (m): max(4, separación mínima deslindes) OGUC 2.6.3
                     5, # ANTEJARDIN (m)
                     2.75, # RASANTE (m de altura / m de separación)
                     60, # 24, #ALTURAMAX (m)
                     30, # 9, #MAXPISOS (unidades)
                     .7, # COEFOCUPACION (m2 / m2 de terreno)
                     1000, # SUBPREDIALMIN (m2)
                     2500, # DENSIDADMAX (Habitantes / 10000 m2 de terreno bruto)
                     true, # FLAGDENSIDADBRUTA
                     4, # COEFCONSTRUCTIBILIDAD (m2 / m2 de terreno)
                     [.5, 1, 1, 1], # ESTACIONAMIENTOSPORVIV (unidades / departamento)
                     .15, # PORCADICESTACVISITAS (unidades / estacionamiento vendible)
                     34, # SUPPORESTACIONAMIENTO (m2 / Estacionamiento)
                     .25, # ESTBICICLETAPORVIV (unidades / estacionamientos totales)
                     3, # BICICLETASPOREST
                     true, # FLAGCAMBIOESTPORBICICLETA
                     10, # MAXSUBTE (unidades)
                     .8, # COEFOCUPACIONEST (m2 / superficieTerreno)
                     7, # SEPESTMIN (m)
                     0 # REDUCCIONESTPORDISTMETRO
                     );

dca = datosCabidaArquitectura(2.55, # 2.625, # ALTURAPISO (m / piso)
                         .2, # PORCSUPCOMUN (m2 / m2 útil)
                         .05, # PORCTERRAZA (m2 / m2 útil)
                         8,#12, # ANCHOMAX (m)
                         [],
                         []
                         );

dcc = datosCabidaComercial([20, 30, 38, 45], # SUPDEPTOUTIL (m2)
                     [.3, .4, .4, 1], # MAXPORCTIPODEPTO
                     [56, 55, 54, 52], # PRECIOVENTA (UF / m2 vendible)
                     200 # PRECIOVENTAEST (UF / unidad)
                     );

dcu = datosCabidaUnit(30, # duracionProyecto
             20, # CostoTerreno
             2 * 1.19 / 100, # ComisionCorredor
             0.60, # Demolicion
             20, # LosaSNT
             10, # LosaBNT
             5, # EstacVisita
             3 / 100, # ExtrasPostVenta
             1.75 / 100, # Arquitectura
             2.75 / 100, # GestionAdministracion
             30, # ITO
             30, # Contabilidad
             0.85 / 100, # Legales
             0.145, # Calculo
             250, # MecanicaSuelo
             50, # Topografia
             .45 / 100, # ProyectoElectrico
             0, # ProyectoSanitario
             0, # ProyectoSeguridad
             50, # ProyectoPavimentacion
             25, # ProyectoBasura
             0, # ProyectoCalefaccion
             1.2 / 100, # Marketing
             0.5 / 100, # PilotSalaVentas
             0, # NotarialesLegales
             0, # CBR
             0, # EmpalmesAportes
             600, # DerechosPermisosMunicipales
             15, # InspeccionCertificacionTecnica
             0, # Contribuciones
             120, # CopiaPlanosArquitectura
             60, # CopiaPlanosObra
             300, # HallAcceso
             400, # JardinesPaisajismo
             0, # HabilitacionSalaUsosMultiples
             0, # HabilitacionPiscinaEnfermeria
             0, # CuentasServiciosDeptos
             0, # GastosComunes
             0, # ContribucionesDeptos
             0, # Seguros
             1 / 100, # ResponsabilidadCivil
             1 / 100, # PostVentaInmobiliaria
             0.1 / 100, # SeguroVentaEnVerde
             1 / 100 # Imprevistos
             );

dcf = datosCabidaFlujo([0.05, 0.05, 0.05, 0.05, 0.10, 0.70], # IngresosVentas
             [1, 0, 0, 0, 0, 0], # CostoTerreno
             [0.05, 0.10, 0.35, 0.35, .15, 0], # CostoConstruccion
             [0.40, 0.15, 0.15, 0.15, .15, 0], # CostoHonorariosProyectos
             [0, 0.10, 0.20, 0.20, 0.25, 0.25], # CostoVenta
             [0.1, 0.40, 0.25, 0.15, 0.05, 0.05], # CostoInmobiliarioObra
             [0, 0, 0, 0, 0.50, 0.50], # CostosHabilitacion
             [0, 0, 0, 0, 0.25, 0.75], # CostosPuestaEnMarcha
             [0.05, 0.05, 0.05, 0.05, 0.05, 0.75], # CostosAtencionCliente
             [0, 0.20, 0.20, 0.20, 0.20, 0.20], # Imprevistos
             [0.05, 0.10, 0.35, 0.35, .15, 0], # ingresoLineaCredito
             [0, 0, 0.10, 0.35, 0.35, .2], # pagoLineaCredito
             .27, # TasaImpuestoRenta
             .05, # TasaInteresLineaCredito
             0.8 * dcu.duracionProyecto, # DuracionLineaCredito
             15 # Tasaciones
             );

dcr = datosCabidaRentabilidad(1.20) # RetornoExigido

if idPredio == 1
    factorCorreccion = 2;
    x = factorCorreccion * [0 30 40 45 10]';
    y = factorCorreccion * [10 0 15 35 30]';
# Calle: [4]
    dcp = datosCabidaPredio(x, y, [1 4], [20 20], 1, 200);

elseif idPredio == 2
    factorCorreccion = 5;
    x = factorCorreccion * [0 10 15 20 25 13 12 4 9]';
    y = factorCorreccion * [10 0 12 7 20 25 22 28 14]';
# Calle: [1 5]
    dcp = datosCabidaPredio(x, y, [1 5], [20 20], 1, 200);

elseif idPredio == 3
# EJEMPLO 1 (Córdova y Figueroa con Vargas Fontecilla, Quinta Normal)
    factorCorreccion = 0.8345;
    x = factorCorreccion * [-7869557.5 -7869559.7 -7869535.4 -7869535.5 -7869490.4 -7869488.3]';
    y = factorCorreccion * [-3952906.8 -3952965.7 -3952966.5 -3952970.9 -3952972.4 -3952909.6]';
# Calle: [1 6]
    dcp = datosCabidaPredio(x, y, [1 6], [18 20], 0, 200);

elseif idPredio == 4
# EJEMPLO 2 (Catedral llegando a General Velásquez, Quinta Normal)
    factorCorreccion = 0.8547;
    x = factorCorreccion * [-7869442.1 -7869413.9 -7869409.5 -7869438.5]';
    y = factorCorreccion * [-3953939.3 -3953935.5 -3953861.6 -3953864.5]';
# Calle: [1]
    dcp = datosCabidaPredio(x, y, [1], [20], 1, 200);

elseif idPredio == 5
    factorCorreccion = 1;
    x = factorCorreccion * [0 150 150 100 100 50 50 0]';
    y = factorCorreccion * [0 0 100 100 50 50 100 100]';
# Calle: [4]
    dcp = datosCabidaPredio(x, y, [1 8], [15 15], 1, 200);

elseif idPredio == 6
    # EJEMPLO 6 (Augusto Leguía)
    factorCorreccion = 0.4;
    x = factorCorreccion * [568 571 612 613 646 648 683 681]';
    y = factorCorreccion * [-339 -405 -403 -418 -419 -452 -442 -355]';
    # Calle: [1 6]
    dcp = datosCabidaPredio(x, y, [1 6 8], [15 10 20], 0, 200);

elseif idPredio == 7
    # EJEMPLO 7 (Independencia)

    dcn.SEPMIN = 5 # SEPMIN (m): max(4, separación mínima deslindes) OGUC 2.6.3
    dcn.ANTEJARDIN = 5 # ANTEJARDIN (m)
    dcn.ALTURAMAX = 30 # 24, #ALTURAMAX (m)
    dcn.MAXPISOS = 20 # 9, #MAXPISOS (unidades)
    dcn.COEFOCUPACION = .5 # COEFOCUPACION (m2 / m2 de terreno)
    dcn.SUBPREDIALMIN = 300 # SUBPREDIALMIN (m2)
    dcn.DENSIDADMAX = 2000 # DENSIDADMAX (Habitantes / 10000 m2 de terreno bruto)
    dcn.COEFCONSTRUCTIBILIDAD = 3 # COEFCONSTRUCTIBILIDAD (m2 / m2 de terreno)

    dca.ANCHOMAX = 10 # ANCHOMAX (m)

    #dcc.SUPDEPTOUTIL = [20, 35, 50, 60] # SUPDEPTOUTIL (m2)
    #dcc.PRECIOVENTA = [65, 58, 53, 50] # PRECIOVENTA (UF / m2 vendible)

    dcc.SUPDEPTOUTIL = [30, 40, 50, 65] # SUPDEPTOUTIL (m2)
    dcc.PRECIOVENTA = [65, 58, 53, 50] # PRECIOVENTA (UF / m2 vendible)
    dcc.MAXPORCTIPODEPTO = [.50, .10, .33, .18] # MAXPORCTIPODEPTO
    dcr = datosCabidaRentabilidad(1.15) # RetornoExigido

    nombreArchivo = "Independencia.csv"
    loadData = CSV.File(string("C:/Users/rjvia/Downloads/", nombreArchivo); header=false)
    numDatos = length(loadData)
    x = zeros(1,numDatos)
    y = zeros(1,numDatos)
    for i = 1:numDatos
        x[i] = loadData[i].Column1
        y[i] = loadData[i].Column2
    end
    areaSup = x[1]
    x = x[2:end]
    y = y[2:end]
    V = [x y]
    factorCorreccion = ajusteArea(V, areaSup)
    dcp = datosCabidaPredio(factorCorreccion * x, factorCorreccion * y, [1 4], [16.5 13.7], 1, 200);

elseif idPredio == 8
    # EJEMPLO 7 (La Florida)

    dcn.SEPMIN = 5 # SEPMIN (m): max(4, separación mínima deslindes) OGUC 2.6.3
    dcn.ANTEJARDIN = 5 # ANTEJARDIN (m)
    dcn.ALTURAMAX = 60 #47 # 24, #ALTURAMAX (m)
    dcn.MAXPISOS = 40 # 9, #MAXPISOS (unidades)
    dcn.COEFOCUPACION = .4 # COEFOCUPACION (m2 / m2 de terreno)
    dcn.SUBPREDIALMIN = 1000 # SUBPREDIALMIN (m2)
    dcn.DENSIDADMAX = 20000 # DENSIDADMAX (Habitantes / 10000 m2 de terreno bruto)
    dcn.COEFCONSTRUCTIBILIDAD = 3.2 # COEFCONSTRUCTIBILIDAD (m2 / m2 de terreno)

    dcc.SUPDEPTOUTIL = [20, 35, 45, 55] # SUPDEPTOUTIL (m2)
    dcc.PRECIOVENTA = [65, 63, 60, 57] # PRECIOVENTA (UF / m2 vendible)
    dcc.MAXPORCTIPODEPTO = [0, 1, 1, 1];

    dca.ANCHOMAX = 8 # Ancho Crujía (m)

    nombreArchivo = "predioLaFlorida2.csv"
    loadData = CSV.File(string("C:/Users/rjvia/Downloads/", nombreArchivo); header=false)
    numDatos = length(loadData)
    x = zeros(1,numDatos)
    y = zeros(1,numDatos)
    for i = 1:numDatos
        x[i] = loadData[i].Column1
        y[i] = loadData[i].Column2
    end
    areaSup = x[1]
    x = x[2:end]
    y = y[2:end]
    V = [x y]
    V = polyShape.reversePath(V);
    x = V[1:end,1]
    y = V[1:end,2]
    factorCorreccion = ajusteArea(V, areaSup)
    dcp = datosCabidaPredio(factorCorreccion * x, factorCorreccion * y, [1 3], [17.5 55], 0, 200);

elseif idPredio == 9
    # EJEMPLO 9 (El Dante)

    dcn.SEPMIN = 8 # SEPMIN (m): max(4, separación mínima deslindes) OGUC 2.6.3
    dcn.ANTEJARDIN = 7 # ANTEJARDIN (m)
    dcn.ALTURAMAX = 52.5 #47 # 24, #ALTURAMAX (m)
    dcn.MAXPISOS = 15 # 9, #MAXPISOS (unidades)
    dcn.COEFOCUPACION = .3 # COEFOCUPACION (m2 / m2 de terreno)
    dcn.SUBPREDIALMIN = 1500 # SUBPREDIALMIN (m2)
    dcn.DENSIDADMAX = 880 # DENSIDADMAX (Habitantes / 10000 m2 de terreno bruto)
    dcn.FLAGDENSIDADBRUTA = false # FLAGDENSIDADBRUTA
    dcn.COEFCONSTRUCTIBILIDAD = 2.8 # COEFCONSTRUCTIBILIDAD (m2 / m2 de terreno)
    dcn.ESTACIONAMIENTOSPORVIV = [1, 2, 2, 2] # ESTACIONAMIENTOSPORVIV
    dcn.FLAGCAMBIOESTPORBICICLETA = false # FLAGCAMBIOESTPORBICICLETA

    dcc.SUPDEPTOUTIL = [50, 90, 110, 140] # SUPDEPTOUTIL (m2)
    dcc.PRECIOVENTA = [95, 91, 89, 87] # PRECIOVENTA (UF / m2 vendible)
    dcc.MAXPORCTIPODEPTO = [1, 1, 1, 1];
    dcc.PRECIOVENTAEST = 350 # PRECIOVENTAEST (UF / unidad)

    dcu.LosaSNT = 30 # LosaSNT
    dcu.LosaBNT = 10 # LosaBNT

    dca.ANCHOMAX = 8 # Ancho Crujía (m)
    dca.ALTURAPISO = 2.7 # 2.625, # ALTURAPISO (m / piso)
    dca.PORCSUPCOMUN = .2 # PORCSUPCOMUN (m2 / m2 útil)

    #nombreArchivo = "el_dante_2.csv"
    nombreArchivo = "predio_ElDante.csv"
    loadData = CSV.File(string("C:/Users/rjvia/Downloads/", nombreArchivo); header=false)
    numDatos = length(loadData)
    x = zeros(1,numDatos)
    y = zeros(1,numDatos)
    for i = 1:numDatos
        x[i] = loadData[i].Column1
        y[i] = loadData[i].Column2
    end
    areaSup = x[1]
    x = x[2:end]
    y = y[2:end]
    V = [x y]
    factorCorreccion = ajusteArea(V, areaSup)
    x = factorCorreccion * x
    y = factorCorreccion * y

#    dcp = datosCabidaPredio(x, y, [1,2,3], [15,15,15], 0, 200);
    dcp = datosCabidaPredio(x, y, [1,2,3,4,5,6], [15,15,15,15,15,20], 0, 200);

end

resultados, ps_predio, ps_base, xopt_cs, fopt_cs, polyCorte,
ps_SombraVolTeor_p, ps_sombraEdif_p,
ps_SombraVolTeor_o, ps_sombraEdif_o,
ps_SombraVolTeor_s, ps_sombraEdif_s = executaCalculoCabidas(dcp, dcn, dca, dcc, dcu, dcf, dcr, fpe, conjuntoTemplates);


text:
# ---
# title: 1550. Three Consecutive Odds
# id: problem1550
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Array
# link: <https://leetcode.com/problems/three-consecutive-odds/description/>
# hidden: true
# ---
#
# Given an integer array `arr`, return `true` if there are three consecutive odd
# numbers in the array. Otherwise, return `false`.
#
#
#
# **Example 1:**
#
#
#
#     Input: arr = [2,6,4,1]
#     Output: false
#     Explanation: There are no three consecutive odds.
#
#
# **Example 2:**
#
#
#
#     Input: arr = [1,2,34,3,4,5,7,23,12]
#     Output: true
#     Explanation: [5,7,23] are three consecutive odds.
#
#
#
#
# **Constraints:**
#
#   * `1 <= arr.length <= 1000`
#   * `1 <= arr[i] <= 1000`
#
#
## @lc code=start
using LeetCode

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


text:
using .Threads, BenchmarkTools

function machin_series(n::Int; handover=false)
    qpi = 0.0
    for i in 1:n
        qpi += (-1)^(i+1)/(2i-1)
        handover && yield()  # we yield here !!!!
    end
    qpi*4
end

function startonthread(id::Int, f::F) where {F<:Function}
    t = Task(nothing)
    @threads for i in 1:nthreads()
        if i == id
            t = @async f()
        end
    end
    fetch(t)
end

println("Benchmark results on thread 1-4 without yielding:")
@btime startonthread(1, ()->machin_series(10_000))
@btime startonthread(2, ()->machin_series(10_000))
@btime startonthread(3, ()->machin_series(10_000))
@btime startonthread(4, ()->machin_series(10_000))

println("Benchmark results on thread 1-4 with yielding:")
@btime startonthread(1, ()->machin_series(10_000, handover=true))
@btime startonthread(2, ()->machin_series(10_000, handover=true))
@btime startonthread(3, ()->machin_series(10_000, handover=true))
@btime startonthread(4, ()->machin_series(10_000, handover=true))


text:
# define tgrid while storing cross-correlations
xcorr_range(tgrid)=range(-1*abs(tgrid[end]-tgrid[1]),stop=abs(tgrid[end]-tgrid[1]),length=2*length(tgrid)-1)

function save(pa::P_fbd, folder; tgridg=nothing, tgrid=nothing, tgrids=nothing)
	!(isdir(folder)) && error("invalid directory")

	save(pa.pfibd,folder,tgridg=tgridg, tgrid=tgrid,tgrids=tgrids)
	save(pa.pfpr,folder,tgridg=tgridg)
	#save(pa.lsbd,tgridg=tgridg, tgrid=tgrid)
end

"""
Save BD
* `tgridg` : 1D grid for Green's functions
* `tgrid` : 1D grid for the data
"""
function save(pa::BD, folder; tgridg=nothing, tgrid=nothing)
	!(isdir(folder)) && error("invalid directory")

	save(pa.om, folder, tgridg=tgridg, tgrid=tgrid,prefix="BD")
	save(pa.optm, folder, tgridg=tgridg, tgrid=tgrid,prefix="BD")
	#savex(pa.optm, folder, tgridg=tgridg, tgrid=tgrid,prefix="BD")

	# finally save err
	err!(pa) # compute err in cal
	file=joinpath(folder, "BDerr.csv")
	CSV.write(file, pa.err)
end

"""
Save IBD
"""
function save(pa::IBD, folder; tgridg=nothing, tgrid=nothing, tgrids=nothing)
	!(isdir(folder)) && error("invalid directory")

	@info "saving results of fibd"
	(tgridg===nothing) && (tgridg = range(0.0, stop=(pa.om.ntg-1)*1.0, length=pa.om.ntg))
	(tgrid===nothing) && (tgrid = range(0.0, stop=(pa.om.nt-1)*1.0, length=pa.om.nt))
	(tgrids===nothing) && (tgrids = range(0.0, stop=(pa.om.nts-1)*1.0, length=pa.om.nts))

	save(pa.om, folder, tgridg=tgridg, tgrid=tgrid, tgrids=tgrids, prefix="IBD")
	save(pa.optm, folder, tgridg=xcorr_range(tgridg), tgrid=xcorr_range(tgrid), prefix="IBD")

	# finally save err
	err!(pa) # compute err in cal
	file=joinpath(folder, "IBDerr.csv")
	CSV.write(file, pa.err)
end

function save(pa::FPR, folder; tgridg=nothing, gobs=nothing)
	!(isdir(folder)) && error("invalid directory")

	@info "saving results of fpr"
	ntg=size(pa.g,1)
	nr=size(pa.g,2)
	(tgridg===nothing) && (tgridg = range(0.0, stop=(ntg-1)*1.0, length=ntg))

	# save observed vs modelled data
#	save_cross(hcat(pa.p_misfit_xcorr.cy...),hcat(pa.p_misfit_xcorr.pxcorr.cg...),
#	    "dfpr", folder)

	file=joinpath(folder, "FPRgfpr.jld")
	#CSV.write(file,DataFrame(hcat(tgridg, pa.g)))
	JLD.save(file, "tgridg", tgridg, "g", pa.g)

	file=joinpath(folder, "FPRimgfpr.csv")
	#CSV.write(file,DataFrame(hcat(repeat(tgridg,outer=nr),repeat(1:nr,inner=ntg),vec(pa.g))),)

	if(!(gobs===nothing))
		file=joinpath(folder, "FPRgfprobs.jld")
		JLD.save(file, "tgridg", tgridg, "g", gobs)
		#CSV.write(file,DataFrame(hcat(tgridg, gobs)))

		file=joinpath(folder, "FPRimgfprobs.csv")
		#CSV.write(file,DataFrame(hcat(repeat(tgridg,outer=nr),repeat(1:nr,inner=ntg),vec(gobs))),)
	end

end

function save(pa::ObsModel, folder; tgridg=nothing, tgrid=nothing, tgrids=nothing, prefix="")
	!(isdir(folder)) && error("invalid directory")

	(tgridg===nothing) && (tgridg = range(0.0, stop=(pa.ntg-1)*1.0, length=pa.ntg))
	(tgrid===nothing) && (tgrid = range(0.0, stop=(pa.nt-1)*1.0, length=pa.nt))
	(tgrids===nothing) && (tgrids = range(0.0, stop=(pa.nts-1)*1.0, length=pa.nts))

	# save original g
	file=joinpath(folder, string(prefix, "gobs.jld"))
	#CSV.write(file,DataFrame(hcat(tgridg, pa.g)))
	JLD.save(file, "tgridg", tgridg, "g", pa.g)

	# save g for imagesc plots in TikZ
	file=joinpath(folder, string(prefix, "imgobs.csv"))
	#CSV.write(file,DataFrame(hcat(repeat(tgridg,outer=pa.nr),repeat(1:pa.nr,inner=pa.ntg),vec(pa.g))),)

	# save original s
	file=joinpath(folder, string(prefix, "sobs.jld"))
	JLD.save(file, "tgrids", tgrids, "s", pa.s)
	CSV.write(file,DataFrame(hcat(tgrids,pa.s)))

	# save data
	file=joinpath(folder, string(prefix, "dobs.jld"))
	JLD.save(file, "tgrid", tgrid, "d", pa.d)
	#CSV.write(file,DataFrame(hcat(tgrid, pa.d)))
end

function save(pa::OptimModel, folder; tgridg=nothing, tgrid=nothing,tgrids=nothing, prefix="")
	!(isdir(folder)) && error("invalid directory")

	(tgridg===nothing) && (tgridg = range(0.0, stop=(pa.ntg-1)*1.0, length=pa.ntg))
	(tgrids===nothing) && (tgrids = range(0.0, stop=(pa.nts-1)*1.0, length=pa.nts))
	(tgrid===nothing) && (tgrid = range(0.0, stop=(pa.nt-1)*1.0, length=pa.nt))

	save_obscal(pa.obs.g, pa.cal.g, tgridg, "goptm", folder, prefix=prefix)
	save_obscal(pa.obs.d, pa.cal.d, tgrid, "doptm", folder,prefix=prefix)
	save_obscal(pa.obs.s, pa.cal.s, tgrids, "soptm", folder,prefix=prefix)

	save_cross(pa.obs.d,pa.cal.d, "doptm", folder,prefix=prefix)
	save_cross(pa.obs.g,pa.cal.g, "goptm", folder,prefix=prefix)
	save_cross(pa.obs.s,pa.cal.s, "soptm", folder,prefix=prefix)
end

function savex(pa::OptimModel, folder; tgridg=nothing, tgrid=nothing, prefix="")
	!(isdir(folder)) && error("invalid directory")

	(tgridg===nothing) && (tgridg = range(0.0, stop=(pa.ntg-1)*1.0, length=pa.ntg))
	(tgrid===nothing) && (tgrid = range(0.0, stop=(pa.nt-1)*1.0, length=pa.nt))

	xtgridg=xcorr_range(tgridg);
	xtgrid=xcorr_range(tgrid)

	xgobs=Conv.xcorr(pa.obs.g)[1]
	xgcal=Conv.xcorr(pa.cal.g)[1]

	save_obscal(xgobs, xgcal, xtgridg, "xg", folder, prefix=prefix)

	# compute cross-correlations without blind Decon
	xdobs=Conv.xcorr(pa.obs.d,Conv.P_xcorr(pa.nt, pa.nr, cglags=[pa.ntg-1, pa.ntg-1]))[1]
	xdcal=Conv.xcorr(pa.cal.d,Conv.P_xcorr(pa.nt, pa.nr, cglags=[pa.ntg-1, pa.ntg-1]))[1]
	xsobs=(Conv.xcorr(pa.obs.s,Conv.P_xcorr(pa.nt, 1, cglags=[pa.ntg-1, pa.ntg-1]))[1])
	xscal=(Conv.xcorr(pa.cal.s,Conv.P_xcorr(pa.nt, 1, cglags=[pa.ntg-1, pa.ntg-1]))[1])

	save_obscal(xdobs, xdcal, xtgridg, "xdat", folder, prefix=prefix)
	save_obscal(xsobs, xscal, xtgridg, "xs", folder, prefix=prefix)

	save_cross(xdobs,xdcal, "xdat", folder, prefix=prefix)
	save_cross(xsobs,xscal, "xs", folder, prefix=prefix)
	save_cross(xgobs,xgcal, "xg", folder, prefix=prefix)
end

function save_cross(a,b, name, folder; num=1000, prefix="")
	(length(a) ≠ length(b)) && error("size mismatch")

	# x plots of data
	nd=length(b);
	fact=(nd>num) ? round(Int,nd/num) : 1
	b=b[1:fact:nd]
	a=a[1:fact:nd]
	datmat=hcat(vec(b), vec(a))
	rmul!(datmat, inv(maximum(abs, datmat)))

	file=joinpath(folder, string(prefix, name, "cross.jld"))
	#CSV.write(file,DataFrame( datmat))
	JLD.save(file, "d", datmat)
end

function save_obscal(a,b, x, name, folder; prefix="")
	(size(a) ≠ size(b)) && error("size mismatch")
	nr=size(a,2)
	nt=size(a,1)

	# save obs
	file=joinpath(folder, string(prefix, name,"obs.jld"))
	JLD.save(file, "t", x, "d", a)
	# save cal
	file=joinpath(folder, string(prefix, name, "cal.jld"))
	JLD.save(file, "t", x, "d", b)

	# save g for imagesc plots in TikZ
	file=joinpath(folder, string(prefix, "im", name, "obs.csv"))
	#CSV.write(file,DataFrame(hcat(repeat(x,outer=nr),repeat(1:nr,inner=nt),vec(a))))

	file=joinpath(folder, string(prefix, "im", name, "cal.csv"))
	#CSV.write(file,DataFrame(hcat(repeat(x,outer=nr),repeat(1:nr,inner=nt),vec(b))))

end