text:
function  magnetic_field(B,x)

# evaluates the radial magnetic field as a function of x = cos(theta);
#
# input
# B   - amplitude of magnetic field (dimensionless)
# x   - cos(theta)
#
# output
# br     - radial field as a function of x

#  size of x grid
   nx = length(x);
   br = zeros(nx);

#  constant radial field
#   br[1:nx] .= B;

#  dipole radial field
#  br = B * x;
#
#  increase rms field toward pole
   Gamma = 2.5;
   B02 = B^2 / (1.0  + Gamma/3.0);
   br2 = B02 *(Gamma*x.*x .+ 1);
   br = sqrt.(br2);

   return br

end

max_stars_repo_name: bbuffett/MacWaves

max_issues_repo_name: bbuffett/MacWaves

max_forks_repo_name: bbuffett/MacWaves


text:
using DifferentiableNAS
using Flux
using Flux: throttle, logitcrossentropy, onecold, onehotbatch, Optimiser
using Zygote: @nograd
using StatsBase: mean
using Parameters
using CUDA
using Distributions
using BSON
using Dates
using Random
include("../CIFAR10.jl")
include("../training_utils.jl")

@show beginscript = now()

@show args = parse_commandline()

if args["random_seed"] > -1
    Random.seed!(args["random_seed"])
end

num_ops = length(PRIMITIVES)

optimiser_α = Optimiser(WeightDecay(1f-3),ADAM(3f-4,(0.5f0,0.999f0)))
optimiser_w = Optimiser(WeightDecay(3f-4),CosineAnnealing(args["epochs"]),Momentum(0.025f0, 0.9f0))
#optimiser_w = Optimiser(WeightDecay(3e-4),Momentum(0.025, 0.9))

train, val = get_processed_data(args["val_split"], args["batchsize"], args["trainval_fraction"], args["random_seed"])
test = get_test_data(args["test_fraction"], args["random_seed"])

histepoch = historiessml()
histbatch = historiessml()
losses = [0f0, 0f0]

base_folder = prepare_folder("admm", args)

cbepoch = CbAll(histepoch, save_progress)
cbbatch = CbAll(histbatch, GC.gc)

function (hist::historiessml)()
    @show losses
    push!(hist.normal_αs_sm, [tanh.(relu.(a)) for a in m.normal_αs |> cpu])
    push!(hist.reduce_αs_sm, [tanh.(relu.(a)) for a in m.reduce_αs |> cpu])
    #push!(hist.activations, copy(m.activations.currentacts) |> cpu)
    push!(hist.train_losses, losses[1])
    push!(hist.val_losses, losses[2])
    #CUDA.reclaim()
    #GC.gc()
end

m = DARTSModelBN(num_cells = args["num_cells"], channels = args["channels"])
zu = ADMMaux(0f0*vcat(m.normal_αs, m.reduce_αs), 0f0*vcat(m.normal_αs, m.reduce_αs))
disc = 7
for epoch in 1:args["epochs"]
    @show epoch
    display(Dates.format(convert(DateTime,now()-beginscript), "HH:MM:SS"))
    if epoch <= 5
        @time DARTStrain1st!(loss, m, train, val, optimiser_α, optimiser_w, losses, epoch; cbepoch = cbepoch, cbbatch = cbbatch)
    else
        @time ADMMtrain1st!(loss, m, train, val, optimiser_w, optimiser_α, zu, args["rho"], losses, epoch, args["epochs"], disc; cbepoch = cbepoch, cbbatch = cbbatch)
        if epoch%5 == 0
            global zu = ADMMaux(0f0*vcat(m.normal_αs, m.reduce_αs), 0f0*vcat(m.normal_αs, m.reduce_αs))
            global disc -= 1
        end
    end
end
display(("done", Dates.format(convert(DateTime,now()-beginscript), "HH:MM:SS")))

ext: jl

max_stars_repo_name: huangyxi/DifferentiableNAS

max_issues_repo_name: huangyxi/DifferentiableNAS

max_forks_repo_name: huangyxi/DifferentiableNAS


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

using WriteVTK
using FileIO
using DataFrames
using JLD2

using CommonUtils
using Basis1D
using Basis2DQuad
using UniformQuadMesh
using UnPack
using Polyester
using SparseArrays

using SetupDG

N = 2
K1D = 40
CFL = 0.9
Np = (N+1)*(N+1)

"Mesh related variables"
VX, VY, EToV = uniform_quad_mesh(2*K1D,K1D)
# @. VX = sin(.5*pi*(VX+1)/2)
# @. VY = ((VY+1)/2)^2/2
@. VX = (VX+1)/2
@. VY = (VY+1)/4
# @. VX = ((VX+1)/2)^(.25)
# @. VY = ((VY+1)/2)^(2)/2

rd = init_reference_quad(N,gauss_lobatto_quad(0,0,N))
md = init_mesh((VX,VY),EToV,rd)
@unpack x,y,mapP,rxJ,syJ,J,sJ,nxJ,nyJ = md;

@unpack VDM,Vf,LIFT,M,Dr,Ds = rd
droptol!(sparse(M),1e-16)
rp,sp = equi_nodes_2D(10)
Vp = vandermonde_2D(N,rp,sp)/VDM

x_vtk = []
y_vtk = []
for k = 1:Int(round(2*K1D))
    push!(x_vtk,x[1,k])
    push!(x_vtk,x[2,k])
    push!(x_vtk,x[3,k])
    # push!(x_vtk,x[4,k])
    # push!(x_vtk,x[5,k])
end
for k = 1:K1D
    # push!(y_vtk,y[1,1+(k-1)*2*K1D])
    # push!(y_vtk,y[5,1+(k-1)*2*K1D])
    # push!(y_vtk,y[9,1+(k-1)*2*K1D])
    # push!(y_vtk,y[13,1+(k-1)*2*K1D])
    # push!(y_vtk,y[1,1+(k-1)*2*K1D])
    # push!(y_vtk,y[6,1+(k-1)*2*K1D])
    # push!(y_vtk,y[11,1+(k-1)*2*K1D])
    # push!(y_vtk,y[16,1+(k-1)*2*K1D])
    # push!(y_vtk,y[21,1+(k-1)*2*K1D])
    push!(y_vtk,y[1,1+(k-1)*2*K1D])
    push!(y_vtk,y[4,1+(k-1)*2*K1D])
    push!(y_vtk,y[7,1+(k-1)*2*K1D])
end
x_vtk = x_vtk .* ones((N+1)*K1D,1)'
y_vtk = y_vtk' .* ones((N+1)*Int(round(2*K1D)),1)

rho_vtk,rhou_vtk,rhov_vtk,E_vtk,gradg_vtk = (zeros(size(x_vtk)) for i = 1:5)

df = load("dg2D_CNS_quad_shocktube.jld2","data")
Uhist = df[end,:Uhist]
thist = df[end,:thist]
# TODO: add savethist to df

pvd = paraview_collection("data/shocktube/N=$N,K1D=$K1D,CFL=$CFL,shocktube", append=true)

for i in 1:length(Uhist)
    t = thist[100*(i-1)+1]
    println("Current time $t")
    flush(stdout)
    vtk_grid("data/shocktube/N=$N,K1D=$K1D,t=$t,CFL=$CFL",x_vtk,y_vtk) do vtk
        rho  = Uhist[i][1,:,:]
        rhou = Uhist[i][2,:,:]
        rhov = Uhist[i][3,:,:]
        E    = Uhist[i][4,:,:]

        # TODO: hardcode rx,ry, put in dataframe
        u     = rhou./rho
        v     = rhov./rho
        vmag  = @. u^2 + v^2
        g     = sqrt.(((rxJ./J).*(Dr*rho)).^2 .+ ((syJ./J).*(Ds*rho)).^2)
        g_min = minimum(g)
        g_max = maximum(g)
        vv    = @. exp(-10*(g-g_min)/(g_max-g_min))

        for k = 1:Int(round(2*K1D))*K1D
            kx = mod1(k,Int(round(2*K1D)))
            ky = div(k-1,Int(round(2*K1D)))+1

            rho_vtk[(kx-1)*(N+1)+1:(kx-1)*(N+1)+N+1, (ky-1)*(N+1)+1:(ky-1)*(N+1)+N+1] = rho[:,k]
            rhou_vtk[(kx-1)*(N+1)+1:(kx-1)*(N+1)+N+1, (ky-1)*(N+1)+1:(ky-1)*(N+1)+N+1] = rhou[:,k]
            rhov_vtk[(kx-1)*(N+1)+1:(kx-1)*(N+1)+N+1, (ky-1)*(N+1)+1:(ky-1)*(N+1)+N+1] = rhov[:,k]
            E_vtk[(kx-1)*(N+1)+1:(kx-1)*(N+1)+N+1, (ky-1)*(N+1)+1:(ky-1)*(N+1)+N+1] = E[:,k]
            gradg_vtk[(kx-1)*(N+1)+1:(kx-1)*(N+1)+N+1, (ky-1)*(N+1)+1:(ky-1)*(N+1)+N+1] = vv[:,k]
        end

        vtk["rho"]   = rho_vtk
        vtk["rhou"]  = rhou_vtk
        vtk["rhov"]  = rhov_vtk
        vtk["E"]     = E_vtk
        vtk["gradg"] = gradg_vtk
        pvd[t] = vtk
    end
end

vtk_save(pvd)

max_stars_repo_name: yiminllin/ESDG-PosLimit

max_issues_repo_name: yiminllin/ESDG-PosLimit

max_forks_repo_name: yiminllin/ESDG-PosLimit


text:
module PowerSystemCaseBuilder

# exports
export SystemCategory
export SystemBuildStats
export SystemDescriptor

export SystemCatalog
export PSYTestSystems
export PSITestSystems
export SIIPExampleSystems
export PSIDTestSystems
export PSSETestSystems
export MatPowerTestSystems

export build_system
export list_categories
export show_categories
export list_systems
export show_systems

export SYSTEM_CATALOG

# imports
import InfrastructureSystems
import InfrastructureSystems: InfrastructureSystemsType
import PowerSystems
import DataStructures: SortedDict
import DataFrames
import PrettyTables

#TimeStamp Management Imports
import TimeSeries
import Dates
import Dates: DateTime, Hour, Minute
import CSV
const PSY = PowerSystems
const IS = InfrastructureSystems

# includes
abstract type PowerSystemCaseBuilderType <: IS.InfrastructureSystemsType end

include("definitions.jl")
include("system_library.jl")

include("system_build_stats.jl")
include("system_descriptor.jl")
include("system_catalog.jl")

include("system_descriptor_data.jl")
include("utils/print.jl")
include("utils/utils.jl")
include("build_system.jl")

end # module

max_stars_repo_name: NREL-SIIP/PowerSystemBuilder

max_issues_repo_name: NREL-SIIP/PowerSystemBuilder

max_forks_repo_name: NREL-SIIP/PowerSystemBuilder


text:
# Julia has the easiest way of doing hello world
println("Hello World")

max_stars_repo_name: kenanchristian/hacktoberfest

max_issues_repo_name: kenanchristian/hacktoberfest

max_forks_repo_name: kenanchristian/hacktoberfest


text:
# This file is a part of EmpiricalDistributions.jl, licensed under the MIT License (MIT).

using EmpiricalDistributions
using Test

using LinearAlgebra

@testset "utils" begin
    @testset "_bin_widths" begin
        @test @inferred(EmpiricalDistributions._bin_widths([2, 3, 5, 6, 10])) == [1, 2, 1, 4]
        @test @inferred(EmpiricalDistributions._bin_widths(2:2:10)) == fill(2, 4)
        @test @inferred(EmpiricalDistributions._bin_widths(2:1.5:10)) == fill(1.5, 5)
    end
end


text:
@noinline function old_cfunction(f, r, a)
  ccall(:jl_function_ptr, Ptr{Cvoid}, (Any, Any, Any), f, r, a)
end


text:
const LT = Union{Diagonal, LowerTriangular}

abstract type AbstractPSDMatrix{T<:Real} <: AbstractMatrix{T} end
struct PSDMatrix{T<:Real} <: AbstractPSDMatrix{T}
    L::LowerTriangular{T,Matrix{T}}
    mat::Matrix{T}
end
function PSDMatrix(mat::Matrix)
    if !(mat ≈ mat') error("Matrix not symmetric") end
    mat = Symmetric(mat)
    vals, U = eigen(mat)
    if any(vals .< 0) error("Matrix not positive semi definite") end
    D = Diagonal(vals)
    # A ≈ U * D * U'
    Q, R = qr(sqrt.(D)*U')
    # A ≈ L*L'
    return PSDMatrix(LowerTriangular(collect(R')))
end
PSDMatrix(L::LowerTriangular) = (nonnegativediag!(L); PSDMatrix(L, L * L'))

Base.size(a::PSDMatrix) = size(a.mat)
Base.getindex(a::PSDMatrix, I::Vararg{Int, N}) where {N} = getindex(a.mat, I...)
Base.copy(A::PSDMatrix) = PSDMatrix(copy(A.L), copy(A.mat))
Base.copy!(dst::PSDMatrix, src::PSDMatrix) = (dst.L .= src.L; dst.mat .= src.mat)

Base.Matrix(a::PSDMatrix) = copy(a.mat)
diag(a::PSDMatrix) = diag(a.mat)
det(Σ::PSDMatrix) = prod(diag(Σ.L).^2)
logdet(Σ::PSDMatrix) = sum(2*log.(diag(Σ.L)))

cholesky(A::PSDMatrix) = cholesky(A.mat)

inv(a::PSDMatrix) = (Li = inv(a.L); Li'Li)
\(a::PSDMatrix, x::AbstractVecOrMat) = a.L' \ (a.L \ x)

function *(c::T, A::PSDMatrix{S}) where {S<:Real, T<:Real}
    if c < 0
        return Matrix(A)*c
    else
        return PSDMatrix(A.L * sqrt(c))
    end
end
*(A::PSDMatrix{S}, c::T) where {S<:Real, T<:Real} = c*A

function X_A_Xt(A::PSDMatrix, X::AbstractMatrix)
    _, R = qr((X*A.L)')
    return PSDMatrix(LowerTriangular(collect(R')))
end
X_A_Xt(A::PSDMatrix, X::LT) = PSDMatrix(X*A.L)

function +(A::PSDMatrix, B::PSDMatrix)
    Q, R = qr([A.L B.L]')
    return PSDMatrix(LowerTriangular(collect(R')))
end

function nonnegativediag!(L::LowerTriangular)
    signs = signbit.(diag(L))
    if !any(signs) return end
    L.data .*= (1 .- 2 .* signs)'
end


text:
using Combinatorics: combinations
using Random

"""
    argmin_rand(rng, v)

Given a vector, find the index of the smallest value. If more than
one value is the smallest, then randomly choose among the smallest
values.
"""
function argmin_rand(rng, v)
    small = typemax(v[1])
    small_extra_cnt = 0
    small_idx = -1
    for i in 1:length(v)
        if v[i] < small
            small = v[i]
            small_extra_cnt = 0
            small_idx = i
        elseif v[i] == small
            small_extra_cnt += 1
        # else not the smallest.
        end
    end
    if small_extra_cnt == 0
        return small_idx
    else
        which = rand(rng, 1:(small_extra_cnt + 1))
        which_cnt = 1
        for s_idx in small_idx:length(v)
            if v[s_idx] == small
                if which_cnt == which
                    return s_idx
                end
                which_cnt += 1
            end
        end
    end
    return 0
end

function n_way_coverage(arity, n_way, M, rng)
    param_cnt = length(arity)
    allc = all_combinations_matrix(arity, n_way)
    maximum_match_score = combination_number(param_cnt, n_way)
    param_cnt = length(arity)

    # Array of arrays.
    coverage = Array{Array{Int64,1},1}()

    trials = zeros(Int, M, param_cnt)
    trial_scores = zeros(Int, M)
    params = zeros(Int, param_cnt)
    entry = zeros(Int, param_cnt)

    loop_idx = 1
    while remaining_uncovered(allc) > 0
        params[:] = 1:param_cnt
        param_coverage = coverage_by_parameter(allc)
        params[1] = argmin_rand(rng, -param_coverage)
        params[params[1]] = 1
        for trial_idx in 1:M
            params[2:end] = shuffle(rng, params[2:end])
            entry[:] .= 0
            candidate_params = coverage_by_value(allc, params[1])
            entry[params[1]] = argmin_rand(rng, -candidate_params)
            for p_idx in 2:param_cnt
                candidate_values = most_matches_existing(
                        allc, entry, params[p_idx]
                        )
                entry[params[p_idx]] = argmin_rand(rng, -candidate_values)
            end
            score = match_score(allc, entry)
            trial_scores[trial_idx] = score
            trials[trial_idx, :] = entry
        end
        chosen_idx = argmin_rand(rng, -trial_scores)
        maximum_score = trial_scores[chosen_idx]
        if maximum_score > 0
            chosen_trial = trials[chosen_idx, :]
            remain = add_coverage!(allc, chosen_trial)
            push!(coverage, chosen_trial)
        end  # nothing was covered, but this can happen
        loop_idx += 1
    end
    coverage
end

function n_way_coverage_init(arity, n_way, seed, M, rng)
    param_cnt = length(arity)
    allc = all_combinations_matrix(arity, n_way)
    # Every combination is nonzero.
    maximum_match_score = combination_number(param_cnt, n_way)
    param_cnt = length(arity)

    # Array of arrays.
    coverage = Array{Array{Int64,1},1}()

    for seed_idx in 1:size(seed, 2)
        push!(coverage, seed[:, seed_idx])
        remain = add_coverage!(allc, seed[:, seed_idx])
    end

    trials = zeros(Int, M, param_cnt)
    trial_scores = zeros(Int, M)
    params = zeros(Int, param_cnt)
    entry = zeros(Int, param_cnt)

    exhausted_by_filter = false
    loop_idx = 1
    while remaining_uncovered(allc) > 0 || exhausted_by_filter
        params[:] = 1:param_cnt
        param_coverage = coverage_by_parameter(allc)
        params[1] = argmin_rand(rng, -param_coverage)
        params[params[1]] = 1
        for trial_idx in 1:M
            params[2:end] = shuffle(rng, params[2:end])
            entry[:] .= 0
            candidate_params = coverage_by_value(allc, params[1])
            entry[params[1]] = argmin_rand(rng, -candidate_params)
            for p_idx in 2:param_cnt
                candidate_values = most_matches_existing(allc, entry, params[p_idx])
                entry[params[p_idx]] = argmin_rand(rng, -candidate_values)
            end
            score = match_score(allc, entry)
            trial_scores[trial_idx] = score
            trials[trial_idx, :] = entry
        end
        chosen_idx = argmin_rand(rng, -trial_scores)
        maximum_score = trial_scores[chosen_idx]
        if maximum_score > 0
            chosen_trial = trials[chosen_idx, :]
            remain = add_coverage!(allc, chosen_trial)
            push!(coverage, chosen_trial)
        else
            exhausted_by_filter = true
        end
        loop_idx += 1
    end
    coverage
end

function arguments_to_arity(arguments)
    [length(x) for x in arguments]
end

function indices_to_arguments(arguments, indices)
    Tuple(arguments[i[1]][i[2]] for i in zip(1:length(indices), indices))
end

"""
    n_way_coverage_filter(arity, n_way, disallow, seed, M, rng)

Greedy coverage generator that includes a filter for tuples
to exclude from consideration.

# Arguments

- `arity`: An array of the number of values for each parameter, in order.
- `n_way`: Integer coverage level > 0.
- `disallow`: A function that returns `true` for values to avoid.
- `seed`: A 2d array of initial values with parameters along the column.
- `M`: The number of tries to get a good test case.
- `rng`: A random number generator.
"""
function n_way_coverage_filter(arity, n_way, disallow, seed, M, rng)
    param_cnt = length(arity)
    allc = all_combinations_matrix(arity, n_way)
    # Every combination is nonzero.
    @assert sum(sum(allc.allc, dims = 2) == 0) == 0
    before_disallow = size(allc.allc, 1)
    remove_combinations!(allc, disallow)
    maximum_match_score = combination_number(param_cnt, n_way)
    param_cnt = length(arity)

    # Array of arrays.
    coverage = Array{Array{Int64,1},1}()

    for seed_idx in 1:size(seed, 2)
        push!(coverage, seed[:, seed_idx])
        remain = add_coverage!(allc, seed[:, seed_idx])
    end

    trials = zeros(Int, M, param_cnt)
    trial_scores = zeros(Int, M)
    params = zeros(Int, param_cnt)
    entry = zeros(Int, param_cnt)

    loop_idx = 1
    while remaining_uncovered(allc) > 0
        params[:] = 1:param_cnt
        param_coverage = coverage_by_parameter(allc)
        params[1] = argmin_rand(rng, -param_coverage)
        params[params[1]] = 1
        trial_cnt = 0
        while trial_cnt < M
            params[2:end] = shuffle(rng, params[2:end])
            entry[:] .= 0
            candidate_params = coverage_by_value(allc, params[1])
            entry[params[1]] = argmin_rand(rng, -candidate_params)
            for p_idx in 2:param_cnt
                candidate_values = most_matches_existing(allc, entry, params[p_idx])
                entry[params[p_idx]] = argmin_rand(rng, -candidate_values)
            end
            if !disallow(entry)
                trial_cnt += 1
                score = match_score(allc, entry)
                trial_scores[trial_cnt] = score
                trials[trial_cnt, :] = entry
            end
        end
        chosen_idx = argmin_rand(rng, -trial_scores[1:trial_cnt])
        chosen_trial = trials[chosen_idx, :]
        maximum_score = trial_scores[chosen_idx]
        if 0 < maximum_score
            remain = add_coverage!(allc, chosen_trial)
            push!(coverage, chosen_trial)
        # else Failure to cover can happen for a bad draw in the shuffle.
        end
        loop_idx += 1
    end
    coverage
end

function n_way_coverage_multi(allc, disallow, seed, M, rng)
    param_cnt = parameter_cnt(allc)
    remove_combinations!(allc, disallow)

    # Array of arrays.
    coverage = Array{Array{Int64,1},1}()

    for seed_idx in 1:size(seed, 2)
        push!(coverage, seed[:, seed_idx])
        add_coverage!(allc, seed[:, seed_idx])
    end

    trials = zeros(Int, M, param_cnt)
    trial_scores = zeros(Int, M)
    params = zeros(Int, param_cnt)
    entry = zeros(Int, param_cnt)

    loop_idx = 1
    while remaining_uncovered(allc) > 0
        params[:] = 1:param_cnt
        param_coverage = coverage_by_parameter(allc)
        params[1] = argmin_rand(rng, -param_coverage)
        params[params[1]] = 1
        trial_cnt = 0
        while trial_cnt < M
            params[2:end] = shuffle(rng, params[2:end])
            entry[:] .= 0
            candidate_params = coverage_by_value(allc, params[1])
            entry[params[1]] = argmin_rand(rng, -candidate_params)
            for p_idx in 2:param_cnt
                candidate_values = most_matches_existing(allc, entry, params[p_idx])
                entry[params[p_idx]] = argmin_rand(rng, -candidate_values)
            end
            if !disallow(entry)
                trial_cnt += 1
                score = match_score(allc, entry)
                trial_scores[trial_cnt] = score
                trials[trial_cnt, :] = entry
            end
        end
        chosen_idx = argmin_rand(rng, -trial_scores[1:trial_cnt])
        chosen_trial = trials[chosen_idx, :]
        maximum_score = trial_scores[chosen_idx]
        if 0 < maximum_score
            add_coverage!(allc, chosen_trial)
            push!(coverage, chosen_trial)
        # else Failure to cover can happen for a bad draw in the shuffle.
        end
        loop_idx += 1
    end
    coverage
end


text:
"""
    cas_parse(str)

given a CAS number string, returns a tuple of 3 Int32 containing the CAS numbers, it performs no validation on the data

"""
function cas_parse(str)
    a1,a2,a3 = split(str,'-')
    n1 = parse(Int32,a1)
    n2 = parse(Int16,a2)
    n3 = parse(Int16,a3)
    return (n1,n2,n3)
end

function unique_idxs_sorted(vec)
    n = length(vec)
    res = zeros(Int,n)
    k = 1
    res[1] = k
    @inbounds for i = 2:n
        if !isequal(vec[i],vec[i-1])
            k+=1
            res[k] = i
        end
    end
    resize!(res,k)
    return res
end
"""
    load_db!(dbtype::Symbol)
    downloads, processes and stores a database corresponding to the one with the same key stored in DATA_INFO

"""

function load_db!(dbtype::Symbol)
    data = DATA_INFO[dbtype]
    if !isfile(data.textdb)
        @info ":" * string(dbtype) * " database file not found, downloading from " * data.url
        url  = data.url
        fname = data.textdb
        path = Downloads.download(url,fname)
        @info ":" * string(dbtype) *  " database file downloaded."
    end

    path = data.textdb
    if !isfile(data.db)
        @info ":" * string(dbtype)  * " arrow file not generated, processing..."
        arrow_db,arrow_synonym_db,arrow_sort_db =parse_and_write_db!(dbtype)
    else
        arrow_db = Arrow.Table(data.db)
        arrow_synonym_db = Arrow.Table(data.symsdb)
        arrow_sort_db = Arrow.Table(data.sorteddb)

    end
    DATA_DB[dbtype] = (arrow_db,arrow_synonym_db,arrow_sort_db)
    return arrow_db,arrow_synonym_db,arrow_sort_db
end

function parse_and_write_db!(dbtype::Symbol)
    data = DATA_INFO[dbtype]
    path = data.textdb
    i = 0
    for _ in eachline(path)
        i +=1
    end

    pubchemid = zeros(Int64,i)
    CAS = Vector{Tuple{Int32,Int16,Int16}}(undef,i)
    formula = Vector{String}(undef,i)
    MW = Vector{Float64}(undef,i)
    smiles = Vector{String}(undef,i)
    InChI = Vector{String}(undef,i)
    InChI_key = Vector{String}(undef,i)
    iupac_name = Vector{String}(undef,i)
    common_name = Vector{String}(undef,i)
    _synonyms  = Vector{Vector{String}}(undef,i)

    i = 0
    for line in eachline(path)
        i += 1
        strs = line |> z->rstrip(z) |> z->split(z,'\t',limit=10)
        pubchemid[i] = parse(Int64,strs[1])
        CAS[i] = cas_parse(strs[2])
        formula[i] = strs[3]
        MW[i] = parse(Float64,strs[4])
        smiles[i] = strs[5]
        InChI[i] = strs[6]
        InChI_key[i] = strs[7]
        iupac_name[i] = strs[8]
        common_name[i] = strs[9]
        if length(strs) >= 10 #not any synonyms
            sym_i = split(strs[10],('\t',';'))
            push!(sym_i,strs[8])
            _synonyms[i]  = sym_i
        else
            _synonyms[i] = String[strs[8]]
        end
    end

    syms_i = mapreduce(length,+,_synonyms)
    synonyms_list = Vector{String}(undef,syms_i)
    synonyms_index = Vector{Int}(undef,syms_i)

    #for some reason,some empty strings are generated as synonyms.
    #those are eliminated here.

    k = 0
    for (ii,sym_vec) in pairs(_synonyms)

        for (jj,sym) in pairs(sym_vec)
            if !isempty(sym)
                k+=1
                synonyms_list[k] = sym
                synonyms_index[k] = ii
            end
        end
    end

    resize!(synonyms_list,k)
    resize!(synonyms_index,k)

    pubchemid_sort =sortperm(pubchemid)
    CAS_sort = sortperm(CAS)
    formula_sort =sortperm(formula)
    MW_sort =sortperm(MW)
    smiles_sort =sortperm(smiles)
    InChI_sort =sortperm(InChI)
    InChI_key_sort =sortperm(InChI_key)
    iupac_name_sort =sortperm(iupac_name)
    common_name_sort =sortperm(common_name)
    synonyms_sort = sortperm(synonyms_list)

    list = synonyms_list[synonyms_sort]
    index = synonyms_index[synonyms_sort]
    #there is the posibility of repeated elements.
    #this is not present in CalebBell/chemicals because it uses a dict, so equal keys
    #store the same value
    list_unique_idx = unique_idxs_sorted(list)
    list = list[list_unique_idx]
    index = index[list_unique_idx]

    db = (;pubchemid, CAS, formula, MW, smiles, InChI, InChI_key, iupac_name, common_name)
    synonym_db = (;list,index)
    sort_db = (;pubchemid_sort, CAS_sort, formula_sort, MW_sort, smiles_sort, InChI_sort, InChI_key_sort, iupac_name_sort, common_name_sort)

    Arrow.write(data.db,db)
    Arrow.write(data.symsdb,synonym_db)
    Arrow.write(data.sorteddb,sort_db)

    arrow_db = Arrow.Table(data.db)
    arrow_synonym_db = Arrow.Table(data.symsdb)
    arrow_sort_db = Arrow.Table(data.sorteddb)
    return arrow_db,arrow_synonym_db,arrow_sort_db
end

"""
    load_data!(key::Symbol;url=nothing,file=nothing)

generates and adds to the global DATA_INFO dict a new database. download and process this database with ´load_db(key)´
"""
function load_data!(key::Symbol;url=nothing,file=nothing)
    if url == file == nothing
        throw(ArgumentError("a file or a url must be provided."))
    elseif (url !== nothing) & (file !== nothing)
        throw(ArgumentError("a file or a url must be provided."))
    elseif (url !== nothing)
        fname = joinpath(download_cache, "pubchem_" * string(key))
        farrow = fname * ".arrow"
        fsyms = fname * "_synonyms.arrow"
        sorted = fname * "_sorted.arrow"
        data = (url=url,textdb=fname,db=farrow,symsdb=fsyms,sorteddb=sorted)

    else #file
        fname = file
        url = fname
        farrow = fname * ".arrow"
        fsyms = fname * "_synonyms.arrow"
        sorted = fname * "_sorted.arrow"
        data = (url=url,textdb=fname,db=farrow,symsdb=fsyms,sorteddb=sorted)
    end
    DATA_INFO[key] = data
end

#3 * x^0.7 - 2 * x + 1
#a,b,c,a
#3,2,1,0.7
#=
ax = bx-c
c = (b-a)x
=#


text:
# Base type from which all solvers inherit from
abstract type AbstractSolver{NDIMS} end

# "eval is evil"
# This is a temporary hack until we have switched to a library based approach
# with pure Julia code instead of parameter files.
# This hack enables users to define functions outside of Trixi.jl and pass
# them as keyword arguments to Trixi.run etc.
eval_if_not_function(func) = eval(Symbol(func))
eval_if_not_function(func::Function) = func

# Create an instance of a solver based on a given name
function make_solver(name::String, equations::AbstractEquation, mesh::TreeMesh;
                     surface_flux_function=nothing, volume_flux_function=nothing)
  if name == "dg"
    # "eval is evil"
    # This is a temporary hack until we have switched to a library based approach
    # with pure Julia code instead of parameter files.
    if isnothing(surface_flux_function)
      surface_flux_type = parameter("surface_flux", "flux_lax_friedrichs")
      surface_flux_function = eval_if_not_function(surface_flux_type)
    end
    if isnothing(volume_flux_function)
      volume_flux_type = parameter("volume_flux", "flux_central")
      volume_flux_function = eval_if_not_function(volume_flux_type)
    end

    initial_conditions_type = parameter("initial_conditions")
    initial_conditions = eval_if_not_function(initial_conditions_type)

    source_terms_type = parameter("source_terms", "nothing")
    source_terms = eval_if_not_function(source_terms_type)

    if ndims(equations) == 2
      return Dg2D(equations, surface_flux_function, volume_flux_function, initial_conditions, source_terms, mesh, parameter("polydeg"))
    elseif ndims(equations) == 3
      return Dg3D(equations, surface_flux_function, volume_flux_function, initial_conditions, source_terms, mesh, parameter("polydeg"))
    else
      error("Unsupported number of spatial dimensions: ", ndims(equations))
    end
  else
    error("'$name' does not name a valid solver")
  end
end

"""
    calc_error_norms([func=(u,equation)->u,] solver, t)

Calculate the discrete L2 and L∞ errors of `func` applied to the conservative variables of
the problem encapsulated by `solver` at time `t`, where `func` is called as `func(u, equation)`.
"""
function calc_error_norms end

@inline calc_error_norms(solver::AbstractSolver, t) = calc_error_norms(cons2cons, solver, t)

####################################################################################################
# Include files with actual implementations for different systems of equations.

# DG
include("dg/dg.jl")


text:
function forward(backend::GPUBackend, state::SubMeanState, input::Blob)
  fea_dim = get_fea_size(input)
  num     = get_num(input)
  CuBLAS.gemm(backend.cublas_ctx, CuBLAS.OP_N, CuBLAS.OP_N, fea_dim, num, 1, convert(eltype(input), -1),
      state.mean_blob.ptr, fea_dim, state.multiplier.ptr, 1, convert(eltype(input), 1),
      input.ptr, fea_dim)
end

function forward(backend::GPUBackend, state::ScaleState, input::Blob)
  CuVec.mul_scal!(backend, input, state.scale)
end


text:
#*************************************************************************#
#
# AbstractType
#
#*************************************************************************#
abstract type ElementwiseOperator <: AparametricModel end

#*************************************************************************#
#
# Generate Buffers
#
#*************************************************************************#
function dx_buffer(f::ElementwiseOperator)
    return Diagonal{Float64}(zeros(f.noutputs))
end

#*************************************************************************#
#
# LogisticFunction
#
#*************************************************************************#
"""
    LogisticFunction(n)

Initialize `ParametricModel` representing logistic function. It applies
`f(x) = 1 / (1+exp(-x))` elementwise to a given input vector.
`n` is the input vector dimension.
"""
struct LogisticFunction <: ElementwiseOperator
    ninputs::Int
    noutputs::Int
end

function LogisticFunction(n::Int)
    LogisticFunction(n, n)
end

function _evaluate!(f::LogisticFunction,
                    x::AbstractVector{<:Real},
                    z::AbstractVector{Float64},
                    dx::AbstractMatrix{Float64})

    z .= 1 ./ (1+exp.(-x))

    if !isempty(dx)
        dx.diag .=  z .* (1-z)
    end
end

#*************************************************************************#
#
# Hyperbolic Tangent
#
#*************************************************************************#
"""
    HyperbolicTangent(n)

Initialize `ParametricModel` representing hyperbolic tangent. It
applies `f(x) = tanh(x)` elementwise to a given input vector.
 `n` is the input vector dimension.
"""
struct HyperbolicTangent <: ElementwiseOperator
    ninputs::Int
    noutputs::Int
end

function HyperbolicTangent(n::Int)
    HyperbolicTangent(n, n)
end

function _evaluate!(f::HyperbolicTangent,
                    x::AbstractVector{<:Real},
                    z::AbstractVector{Float64},
                    dx::AbstractMatrix{Float64})
    z .= tanh.(x)

    if !isempty(dx)
        dx.diag .=  1 - z.^2
    end
end

#*************************************************************************#
#
# AffineMap
#
#*************************************************************************#
"""
    AffineMap(a[, b])

Initialize `ParametricModel` representing affine function. It applies
`f(x) = a*x + b` elementwise to a given input vector. `a` and
`b` are vectors.
"""
struct AffineMap <: ElementwiseOperator
    ninputs::Int
    noutputs::Int
    a::Vector{Float64}
    b::Vector{Float64}
end

function AffineMap{R<:Real, S<:Real}(a::Vector{R}, b::Vector{S}=zeros(length(a)))
    if length(a) != length(b)
        throw(ArgumentError("Vector `a` and `b` should have the same size."))
    end

    n = length(a)
    AffineMap(n, n, a, b)
end

function _evaluate!(f::AffineMap,
                    x::AbstractVector{<:Real},
                    z::AbstractVector{Float64},
                    dx::AbstractMatrix{Float64})
    z .= f.a.*x + f.b

    if !isempty(dx)
        dx.diag .=  f.a
    end
end

#*************************************************************************#
#
# Identity
#
#*************************************************************************#
"""
    Identity(n)

Initialize `ParametricModel` representing identity function.  It applies
`f(x) = x` elementwise to a given input vector. `n` is the input vector
dimension.
"""
struct Identity <: ElementwiseOperator
    ninputs::Int
    noutputs::Int
end

function Identity(n::Int)
    Identity(n, n)
end

function _evaluate!(f::Identity,
                    x::AbstractVector{<:Real},
                    z::AbstractVector{Float64},
                    dx::AbstractMatrix{Float64})
    z .= x

    if !isempty(dx)
        dx.diag .=  1
    end
end


text:
#################### Approximate Bayesian Computation ####################

#################### Types and Constructors ####################

type ABCTune
  datakeys::Vector{Symbol}
  Tsim::Vector{Vector{Float64}}
  epsilon::Vector{Float64}
  epsilonprime::Vector{Float64}

  function ABCTune()
    new(
      Vector{Symbol}(),
      Vector{Vector{Float64}}(),
      Vector{Float64}()
    )
  end
end

#################### Sampler Constructor ####################

function ABC{T<:Real}(params::ElementOrVector{Symbol},
                      scale::ElementOrVector{T}, summary::Function,
                      epsilon::Real; kernel::KernelDensityType=SymUniform,
                      dist::Function=(Tsim, Tobs) -> sqrt(sumabs2(Tsim - Tobs)),
                      proposal::SymDistributionType=Normal, maxdraw::Integer=1,
                      nsim::Integer=1, decay::Real=1.0, randeps::Bool=false,
                      args...)
  0 <= decay <= 1 || throw(ArgumentError("decay is not in [0, 1]"))

  params = asvec(params)
  kernelpdf = (epsilon, d) -> pdf(kernel(0.0, epsilon), d)
  local obsdata, simdata, summarizenodes

  samplerfx = function(model::Model, block::Integer)
    tune = gettune(model, block)

    ## current parameter and density values
    theta0 = unlist(model, block, true)
    logprior0 = logpdf(model, params, true)
    pi_epsilon0 = 0.0
    pi_error0 = 1.0

    ## initialize tuning parameters
    local Tobs
    if model.iter == 1
      ## data node symbols
      targets = keys(model, :target, params)
      stochastics = keys(model, :stochastic)
      tune.datakeys = intersect(setdiff(targets, params), stochastics)

      ## local functions to get and summarize data nodes
      obsdata = key -> unlist(model[key])
      simdata = key -> unlist(model[key], rand(model[key]))
      summarizenodes = length(tune.datakeys) > 1 ?
        data -> vcat(map(key -> summary(data(key)), tune.datakeys)...) :
        data -> asvec(summary(data(tune.datakeys[1])))

      ## observed data summary statistics
      Tobs = summarizenodes(obsdata)

      tune.Tsim = Array{Vector{Float64}}(nsim)
      tune.epsilon = Array{Float64}(nsim)
      tune.epsilonprime = Array{Float64}(nsim)
      for i in 1:nsim
        ## simulated data summary statistics for current parameter values
        tune.Tsim[i] = summarizenodes(simdata)
        d = dist(tune.Tsim[i], Tobs; args...)

        ## starting tolerance
        tune.epsilon[i] = decay > 0 ? max(epsilon, d) : epsilon
        if randeps
          dexp = Exponential(tune.epsilon[i])
          tune.epsilonprime[i] = rand(dexp)
          pi_error0 = pdf(dexp, tune.epsilonprime[i])
        else
          tune.epsilonprime[i] = tune.epsilon[i]
        end

        ## kernel density evaluation
        pi_epsilon0 += kernelpdf(tune.epsilonprime[i], d) * pi_error0
      end
    else
      Tobs = summarizenodes(obsdata)
      for i in 1:nsim
        d = dist(tune.Tsim[i], Tobs; args...)
        if randeps
          dexp = Exponential(tune.epsilon[i])
          pi_error0 = pdf(dexp, tune.epsilonprime[i])
        end
        pi_epsilon0 += kernelpdf(tune.epsilonprime[i], d) * pi_error0
      end
    end

    Tsim1 = similar(tune.Tsim)
    epsilon1 = similar(tune.epsilon)
    epsilonprime1 = similar(tune.epsilonprime)

    for k in 1:maxdraw
      ## candidate draw and prior density value
      theta1 = theta0 + scale .* rand(proposal(0.0, 1.0), length(theta0))
      logprior1 = mapreduce(keyval -> logpdf(model[keyval[1]], keyval[2]), +,
                            relist(model, theta1, block, true))
      if logprior1 == -Inf
        continue
      end
      relist!(model, theta1, block, true)

      ## tolerances and kernel density
      pi_epsilon1 = 0.0
      pi_error1 = 1.0
      for i in 1:nsim
        ## simulated data summary statistics for candidate draw
        Tsim1[i] = summarizenodes(simdata)
        d = dist(Tsim1[i], Tobs; args...)

        ## monotonically decrease tolerance to target
        epsilon1[i] = (1 - decay) * tune.epsilon[i] +
                      decay * max(epsilon, min(d, tune.epsilon[i]))
        if randeps
          dexp = Exponential(epsilon1[i])
          epsilonprime1[i] = rand(dexp)
          pi_error1 = pdf(dexp, epsilonprime1[i])
        else
          epsilonprime1[i] = epsilon1[i]
        end

        ## kernel density evaluation
        pi_epsilon1 += kernelpdf(epsilonprime1[i], d) * pi_error1
      end

      ## accept/reject the candidate draw
      if rand() < pi_epsilon1 / pi_epsilon0 * exp(logprior1 - logprior0)
        theta0 = theta1
        tune.Tsim = Tsim1
        tune.epsilon = epsilon1
        tune.epsilonprime = epsilonprime1
        break
      end
    end

    relist(model, theta0, block, true)
  end

  Sampler(params, samplerfx, ABCTune())
end


text:
# Launch with `julia --project --track-allocation=user`
if !("." in LOAD_PATH)
    push!(LOAD_PATH, ".")
end
import Profile

include("common.jl")
constructor = ENV["ALLOCATION_CONSTRUCTOR"]
@info "Recording allocations for $constructor"

thermo_state_map = Dict(
    "ρeq" => thermo_state_ρeq,
    "pθq" => thermo_state_pθq,
    "pTq" => thermo_state_pTq,
)
thermo_state = thermo_state_map[constructor]

thermo_state() # compile first
Profile.clear_malloc_data()
thermo_state()

# Quit julia (which generates .mem files), then call
#=
import Coverage
allocs = Coverage.analyze_malloc("src")
=#


text:
include("libswresample_h.jl")

#include("swresample.jl")


text:
import AQFED.TermStructure: ConstantBlackModel
import AQFED.Black: blackScholesFormula
import Roots: find_zero, A42
export AndersenLakeNRepresentation, priceAmerican

#Andersen-Lake American option pricing under negative rates
struct AndersenLakeNRepresentation
    isCall::Bool
    model::ConstantBlackModel
    tauMax::Float64
    tauMaxOrig::Float64
    tauHat::Float64
    nC::Int
    nTS1::Int
    nTS2::Int
    capX::Float64
    capXD::Float64
    avec::Vector{Float64}
    avecD::Vector{Float64}
    qvec::Vector{Float64}
    qvecD::Vector{Float64}
    wvec::Vector{Float64}
    yvec::Vector{Float64}
end

function AndersenLakeNRepresentation(
    model::ConstantBlackModel,
    tauMax::Float64,
    atol::Float64,
    nC::Int,
    nIter::Int,
    nTS1::Int,
    nTS2::Int;
    isCall::Bool = false
)
    aUp = AndersenLakeRepresentation(model, tauMax, atol, nC, nIter, nTS1, nTS2, isCall = isCall, isLower = false)
    if (model.r < 0) && (model.q < model.r)
        aDown = AndersenLakeRepresentation(model, tauMax, atol, nC, nIter, nTS1, nTS2, isCall = isCall, isLower = true)
        tauHat = aUp.tauHat
        tauStar = tauHat
        #calculate intersection tauStar
        logCapXD = log(aDown.capX)
        logCapX = log(aUp.capX)
        logBdown = logCapXD + sqrt(aDown.qvec[1])
        logBup = logCapX - sqrt(aUp.qvec[1])
        if logBdown > logBup
            obj = function (τ)
                z = 2 * sqrt((τ) / tauHat) - 1
                qck = chebQck(aUp.avec, z)
                lnBUp = logCapX - sqrt(qck)
                qck = chebQck(aDown.avec, z)
                lnBDown = logCapXD + sqrt(qck)
                return lnBUp - lnBDown
            end
            tauStar = find_zero(obj, (0, tauHat), A42())
            #  println("tauStar ", tauStar)
        end
        return AndersenLakeNRepresentation(isCall, model, tauStar, tauMax, tauHat, nC, nTS1, nTS2,
            aUp.capX, aDown.capX, aUp.avec, aDown.avec, aUp.qvec, aDown.qvec, aUp.wvec, aUp.yvec)
    else
        return AndersenLakeNRepresentation(isCall, model, tauMax, tauMax, tauMax, nC, nTS1, nTS2,
            aUp.capX, Float64[], aUp.avec, Float64[], aUp.qvec, Float64[], aUp.wvec, aUp.yvec)
    end
end

function priceAmerican(p::AndersenLakeNRepresentation, K::Float64, S::Float64)::Float64
    if isempty(p.qvecD)
        return priceAmerican(AndersenLakeRepresentation(p.isCall, p.model, p.tauMax, p.tauMax, p.nC, p.nTS1, p.nTS2, p.capX, p.avec, p.qvec, p.wvec, p.yvec), K, S)
    end
    vol, r, q = p.model.vol, p.model.r, p.model.q
    if p.isCall #use McDonald and Schroder symmetry
        K, S = S, K
        r, q = q, r
    end
    capX, capXD = p.capX * K, p.capXD * K

    f0 = exp(-sqrt(p.qvec[1])) * capX
    f0D = exp(sqrt(p.qvecD[1])) * capXD
    if S <= f0 && S >= f0D && p.tauMax == p.tauMaxOrig
        # println(f0, " ",f0D)
        return max(K - S, 0.0)
    end

    tauMax, tauMaxOrig, tauHat, nTS2 = p.tauMax, p.tauMaxOrig, p.tauHat, p.nTS2
    wvec, yvec, avec, avecD = p.wvec, p.yvec, p.avec, p.avecD
    nC, rK, qS = p.nC, r * K, q * S

    uMax = tauMax
    uMin = 0.0
    uScale = (uMax - uMin) / 2
    uShift = (uMax + uMin) / 2
    sum4k = 0.0
    isCrossed = false
    euro = blackScholesFormula(
        false,
        K,
        S,
        vol * vol * tauMaxOrig,
        exp(-(r - q) * tauMaxOrig),
        exp(-r * tauMaxOrig),
    )
    if tauMax == tauMaxOrig
        #only 1 integral to compute
        isLower = S < f0D
        if isLower
            avec = avecD
            capX = capXD
        end
        for sk2 = nTS2:-1:1
            wk = wvec[sk2]
            yk = yvec[sk2]
            uk = uScale * yk + uShift
            if abs(yk) != 1
                zck = 2 * sqrt(uk / tauHat) - 1 #cheb from tauHat.
                qck = chebQck(avec, zck)
                Bzk = isLower ? capX * exp(sqrt(qck)) : capX * exp(-sqrt(qck))
                tauk = uMax - uk
                d1k, d2k = vaGBMd1d2(S, Bzk, r, q, tauk, vol)
                sum4k += wk * rK * exp(-r * tauk) * normcdf(-d2k)
                sum4k += -wk * qS * exp(-q * tauk) * normcdf(-d1k)
            end
        end
        if isLower
            sum4k = -sum4k
            euro = K-S
        end
    else
        # two integrals
        for sk2 = nTS2:-1:1
            wk = wvec[sk2]
            yk = yvec[sk2]
            uk = uScale * yk + uShift
            if abs(yk) != 1
                zck = 2 * sqrt(uk / tauHat) - 1 #cheb from tauHat.
                qck = chebQck(avec, zck)
                qckD = chebQck(avecD, zck)
                Bzk = capX * exp(-sqrt(qck))
                BzkD = capXD * exp(sqrt(qckD))
                if Bzk <= BzkD
                    isCrossed = true
                end
                if !isCrossed
                    tauk = uMax - uk + tauMaxOrig - tauMax
                    d1k, d2k = vaGBMd1d2(S, Bzk, r, q, tauk, vol)
                    sum4k += wk * rK * exp(-r * tauk) * normcdf(-d2k)
                    sum4k += -wk * qS * exp(-q * tauk) * normcdf(-d1k)
                    d1k, d2k = vaGBMd1d2(S, BzkD, r, q, tauk, vol)
                    sum4k -= wk * rK * exp(-r * tauk) * normcdf(-d2k)
                    sum4k += wk * qS * exp(-q * tauk) * normcdf(-d1k)
                end
            end
        end
    end

    price = euro + uScale * sum4k
    price = max(K - S, price)
    return price
end


text:
__precompile__()

module ObjGeom

import Base: isvalid
import Statistics: mean
import LinearAlgebra: cross, norm, dot, normalize!

export load_obj, get_indexed, isvalid, add_normals, add_values, delete_values_at, delete_unused_values,
       regularpoly, prism, pyramid, sphere,
       ObjModel, SMOOTHING

const defaultIds = [:position, :texCoord, :normal]

const ITuple{N} = NTuple{N, Int}

mutable struct ObjModel
  values::Vector{Matrix} # vector of data streams, each containing a type of values (positions, texture coordinates, normals)
  faces::Vector{Matrix{Int}} # vector of face data, each face containing indices inside values for each vertex
  valueIds::Vector{Symbol}

  ObjModel(ids, values, faces) = new(values, faces, copy(ids))
  ObjModel() = new([:position], [Matrix{Float32}(undef, 3, 0)], Matrix{Int}[])
end

matrixcols(v::Vector, dim::Int) = div(length(v), dim)
asmatrix(v::Vector, dim::Int) = reshape(v, dim, matrixcols(v, dim))

function value_matrix(values::Vector{Float32}, default::Vector{Float32})
  mat = asmatrix(values, length(default))
  rows = size(mat, 1)
  for r = rows:-1:1
    any(mat[r, i] != default[r] for i = 1:size(mat, 2)) && break
    rows = r
  end
  if rows < length(default)
    mat = mat[1:rows, :]
  end
  mat
end

load_obj(file::String) = open(io->load_obj(io), file)

const defaultVecs = [ [0f0, 0f0, 0f0, 1f0], [0f0, 0f0, 0f0], [0f0, 0f0, 0f0] ]

function load_obj(io::IOStream)
  values = [Float32[] for i = 1:3]
  smoothGroups = Dict{Int, Vector{Int}}()
  smoothGroupId = 0
  faces = Vector{Matrix{Int}}()

  for l in eachline(io)
    tokens = split(strip(l))
    comment = findfirst(t->t[1] == '#', tokens)
    if comment != nothing
      tokens = tokens[1:comment-1]
    end
    isempty(tokens) && continue
    if tokens[1] == "s"
      smoothGroupId = tokens[2] == "off" ? 0 : parse(Int, tokens[2])
    elseif tokens[1] == "v" || tokens[1] == "vt" || tokens[1] == "vn"
      valInd = tokens[1] == "v" ? 1 : (tokens[1] == "vt" ? 2 : 3)
      vals = map(v->parse(Float32, v), tokens[2:end])
      append!(vals, defaultVecs[valInd][length(vals)+1:end])
      append!(values[valInd], vals)
    elseif tokens[1] == "f"
      vals=map(tokens[2:end]) do t
        indices = split(t, '/')
        indVals = Matrix{Int}(undef, 3, 1)
        for i = 1:3
          ind = (i > length(indices) || isempty(indices[i])) ? 0 : parse(Int, indices[i])
          if ind < 0
            ind = matrixcols(values[i], length(defaultVecs[i])) + 1 + ind
          end
          indVals[i, 1] = ind
        end
        indVals
      end
      face = hcat(vals...)
      push!(faces, face)
      if smoothGroupId != 0
        smoothGroup = get!(smoothGroups, smoothGroupId) do; Int[] end
        push!(smoothGroup, length(faces))
      end
    end
  end
  objValues = Matrix[value_matrix(values[1], [NaN32, NaN32, NaN32, 1.0f0]),
                     value_matrix(values[2], [NaN32, 0f0, 0f0]),
                     value_matrix(values[3], [NaN32, NaN32, NaN32])]
  model = ObjModel(defaultIds, objValues, faces)
  add_normals(model, Base.values(smoothGroups))
  delete_unused_values(model)
  return model
end

# ObjModel

function isvalid(model::ObjModel, checkIndices::Bool = false)
  dim = length(model.values)
  if length(model.valueIds) != dim || any(face->size(face, 1) != dim, model.faces)
    return false
  end
  return !checkIndices ||
         all(model.faces) do face
           all(1 <= face[i, j] <= length(model.values[i]) for j = 1:size(face, 2), i = 1:size(face, 1))
         end
end

function facenormal(model::ObjModel, faceInd::Int)
  posInd = findfirst(isequal(:position), model.valueIds)
  positions = model.values[posInd]
  face = model.faces[faceInd]
  n = zeros(Float32, 3)
  for i = 1:size(face, 2)-2
    v1 = positions[1:3, face[posInd, i]]
    v2 = positions[1:3, face[posInd, i+1]]
    v3 = positions[1:3, face[posInd, i+2]]
    n = cross(v1-v2, v2-v3)
    len = norm(n)
    if len > eps(len)
      n /= len
      break
    end
  end
  n
end

function facenormals(model::ObjModel)
  normals = Matrix{Float32}(undef, 3, length(model.faces))
  for i = 1:length(model.faces)
    normals[:, i] = facenormal(model, i)
  end
  normals
end

averagenormals(normals::Matrix) = mean(normals; dims = 2) # maybe change to calculate axis of bounding cone?

function build_vertex_faces(model::ObjModel)
  posInd = findfirst(isequal(:position), model.valueIds)
  normalsInd = findfirst(isequal(:normal), model.valueIds)
  vertexFaces = Dict{Int, Vector{Int}}() # map between vertex indices and faces they are part of
  for faceInd = 1:length(model.faces)
    face = model.faces[faceInd]
    for i = 1:size(face, 2)
      if face[normalsInd, i] == 0 # face vertex has no specified normal
        faceArray = get!(vertexFaces, face[posInd, i]) do; Int[] end
        if isempty(faceArray) || faceArray[end] != faceInd # check if this vertex index wasn't specified already in this face's vertices
          push!(faceArray, faceInd)
        end
      end
    end
  end
  vertexFaces
end

function add_normals(model::ObjModel, smoothGroups = Vector{Int}[])
  posInd = findfirst(isequal(:position), model.valueIds)
  normalsInd = findfirst(isequal(:normal), model.valueIds)
  if normalsInd == nothing
    normalsInd = add_values(model, :normal, 3)
  end
  vertexFaces = build_vertex_faces(model)
  faceNormals = facenormals(model)
  normals = model.values[normalsInd]
  normals = reshape(normals, length(normals))
  addedNormalIndices = Dict{ITuple, Int}()
  # add normals to vertices on faces inside smoothing groups
  for smoothFaces in smoothGroups, faceInd in smoothFaces
    face = model.faces[faceInd]
    for i = 1:size(face, 2)
      if face[normalsInd, i] == 0
        vertInd = face[posInd, i]
        faceArray = intersect(vertexFaces[vertInd], smoothFaces) # all faces in the current smoothing group that this vertex is part of
        faceIndexTuple = (faceArray...,)
        face[normalsInd, i] = get!(addedNormalIndices, faceIndexTuple) do
          normalsToAverage = faceNormals[:, faceArray]
          normal = averagenormals(normalsToAverage)
          append!(normals, normal)
          return matrixcols(normals, 3)
        end
      end
    end
  end
  # add normals to faces outside smoothing groups
  for faceInd in 1:length(model.faces)
    face = model.faces[faceInd]
    for i = 1:size(face, 2)
      if face[normalsInd, i] == 0
        face[normalsInd, i] = get!(addedNormalIndices, (faceInd,)) do
          append!(normals, faceNormals[:, faceInd])
          return matrixcols(normals, 3)
        end
      end
    end
  end
  model.values[normalsInd] = asmatrix(normals, 3)
  nothing
end

function add_texcoord(model::ObjModel, direction::Vector{Float32} = [0f0, 0f0, 1f0], uniform::Bool = true)
  posInd = findfirst(isequal(:position), model.valueIds)
  tcInd = findfirst(isequal(:texCoord), model.valueIds)
  if tcInd == nothing
    tcInd = add_values(model, :texCoord, 2)
  end
  normalize!(direction)
  xDir = [1f0, 0f0, 0f0]
  zDir = [0f0, 0f0, 1f0]
  other = abs(dot(direction, xDir)) < abs(dot(direction, zDir)) ? xDir : zDir
  u = cross(direction, other)
  v = cross(u, direction)
  proj = vcat(u', v', direction')
  uvw = proj * model.values[posInd]
  ext = extrema(uvw, dims=2)
  if uniform
    mini = minimum(map(t->t[1], ext))
    maxi = maximum(map(t->t[2], ext))
    ext = fill((mini, maxi), 3)
  end
  baseInd = length(model.values[tcInd])
  newTexCoord = [(uvw[i, j] - ext[i][1]) / (ext[i][2] - ext[i][1]) for i = 1:min(size(uvw, 1), size(model.values[tcInd], 1)), j = 1:size(uvw, 2)]
  model.values[tcInd] = hcat(model.values[tcInd], newTexCoord)
  for face in model.faces, i = 1:size(face, 2)
    if face[tcInd, i] == 0
      face[tcInd, i] = face[posInd, i] + baseInd
    end
  end
  model
end

function append_face(indices::Vector{Int}, face::Vector{Int})
  # push the face as a fan in triangle list format
  for i = 1:length(face) - 2
    push!(indices, face[1], face[i+1], face[i+2])
  end
end

const defaultVals = [0, 0, 0, 1]

function get_indexed(model::ObjModel)
  @assert length(model.values) == length(model.valueIds)
  generatedVertices = Dict{ITuple, Int}()
  streams = map(model.values) do v
    Vector{eltype(v)}()
  end
  indices = Int[]
  for face in model.faces
    @assert size(face, 1) == length(model.values)
    faceIndices = Int[]
    for i = 1:size(face, 2)
      indTuple = (face[:, i]...,)
      vertInd = get!(generatedVertices, indTuple) do
        for s = 1:length(indTuple)
          value = indTuple[s] != 0 ?
                    model.values[s][:, indTuple[s]] :
                    defaultVals[1:size(model.values[s], 1)]
          append!(streams[s], value)
        end
        return matrixcols(streams[1], size(model.values[1], 1)) - 1
      end
      push!(faceIndices, vertInd)
    end
    append_face(indices, faceIndices)
  end
  vertices = Dict{Symbol, Array}(model.valueIds[s] => asmatrix(streams[s], size(model.values[s], 1)) for s = 1:length(model.values))
  return vertices, indices
end

function add_values(model::ObjModel, id::Symbol, dim::Int; elType::DataType = Float32)
  @assert isvalid(model)
  push!(model.values, Matrix{elType}(undef, dim, 0))
  push!(model.valueIds, id)
  map!(model.faces, model.faces) do face
    vcat(face, fill(0, 1, size(face, 2)))
  end
  return length(model.values)
end

function delete_values_at(model::ObjModel, ind::Int)
  @assert isvalid(model)
  deleteat!(model.values, ind)
  deleteat!(model.valueIds, ind)
  map!(model.faces, model.faces) do face
    face[[i!=ind for i = 1:size(face, 1)], :]
  end
  nothing
end

function delete_unused_values(model::ObjModel)
  i = 1
  while i <= length(model.values)
    if isempty(model.values[i])
      delete_values_at(model, i)
    else
      i += 1
    end
  end
  nothing
end

# Model generation

@enum SMOOTHING SmoothNone=0 SmoothSides=1 SmoothAll=2

function regularpoly(sides::Int, z::Float32 = 0f0; normals = true)
	points = Matrix{Float32}(undef, 3, sides)
	for i in 1:sides
		ang = (i - 1) * 2pi / sides
		s = sin(ang)
		c = cos(ang)
		points[:, i] = [cos(ang), sin(ang), z]
	end

  face = [j for i=1:1, j=sides:-1:1] # Matrix with a single row

  model = ObjModel([:position], [points], [face])
  if normals
    add_normals(model)
  end
  model
end

clampindex(i::Int, n::Int) = (i-1)%n+1

function prism(sides::Int, zMin::Float32 = -0.5f0, zMax::Float32 = 0.5f0; smooth::SMOOTHING = SmoothNone)
  model = regularpoly(sides, zMin, normals = false)
  upperPoints = model.values[1].+[0, 0, zMax - zMin]
  upperFace = model.faces[1][:, end:-1:1].+sides
  model.values[1] = hcat(model.values[1], upperPoints)
  push!(model.faces, upperFace)
  smoothGroup = Int[]
  if smooth == SmoothAll
    push!(smoothGroup, 1, 2)
  end
  for i = 1:sides
    iNext = clampindex(i+1, sides)
    side = [i iNext iNext+sides i+sides]
    push!(model.faces, side)
    if smooth != SmoothNone
      push!(smoothGroup, length(model.faces))
    end
  end

  add_normals(model, [smoothGroup])
  model
end

function pyramid(sides::Int, zMin::Float32 = 0f0, zMax::Float32 = 1f0; smooth::SMOOTHING = SmoothNone)
  model = regularpoly(sides, zMin, normals = false)
  model.values[1] = hcat(model.values[1], [0, 0, zMax])
  smoothGroup = Int[]
  if smooth == SmoothAll
    push!(smoothGroup, 1)
  end
  apexInd = size(model.values[1], 2)
  for i = 1:sides
    side = [i clampindex(i+1, sides) apexInd]
    push!(model.faces, side)
    if smooth != SmoothNone
      push!(smoothGroup, length(model.faces))
    end
  end

  add_normals(model, [smoothGroup])
  model
end

function sphere(sides::Int, rh::Float32 = 1f0, rv::Float32 = 1f0; smooth::SMOOTHING = SmoothAll)
  vsides = div(sides, 2) - 1
  points = Matrix{Float32}(undef, 3, sides*vsides+2)
  points[:, end] = [0, 0, rv]
  points[:, end-1] = [0, 0, -rv]
  faces = Matrix{Int}[]
  smoothGroup = Int[]
  ptsCount = size(points, 2)
  for x = 0:sides-1
    nextX = (x+1) % sides
    baseX = vsides*x
    nextBaseX = vsides*nextX
    angH = x*2pi/sides
    sinH = sin(angH)
    cosH = cos(angH)
    for y = 1:vsides
      angV = y*pi/(vsides+1)
      sinV = sin(angV)
      cosV = cos(angV)
      points[:, baseX+y] = Float32[rh*sinV*cosH, rh*sinV*sinH, rv*cosV]
    end

    for y = 1:vsides-1
      push!(faces, [nextBaseX+y baseX+y baseX+y+1 nextBaseX+y+1])
      if smooth != SmoothNone
        push!(smoothGroup, length(faces))
      end
    end
    push!(faces, [baseX+1 nextBaseX+1 ptsCount])
    push!(faces, [nextBaseX+vsides baseX+vsides ptsCount-1])
    if smooth != SmoothNone
      i = length(faces)
      push!(smoothGroup, i-1, i)
    end
  end

  model = ObjModel([:position], [points], faces)
  add_normals(model, [smoothGroup])
  model
end

end


text:
import LazySets: dim, inner
using LazySets: isapproxzero

"""
    AbstractDirections{N}

Abstract type for template direction representations.

### Notes

All subtypes should implement the standard iterator methods from `Base` and the
function `dim(d<:AbstractDirections)::Int`.
"""
abstract type AbstractDirections{N} end

"""
    dim(ad::AbstractDirections)::Int

Returns the dimension of the generated directions.

### Input

- `ad` -- box direction representation

### Output

The dimension of the generated directions.
"""
function dim(ad::AbstractDirections)::Int
    return ad.n
end

# ==================================================
# Box directions
# ==================================================

"""
    UnitVector{T} <: AbstractVector{T}

A lazy unit vector with arbitrary one-element.

### Fields

- `i` -- index of non-zero entry
- `n` -- vector length
- `v` -- non-zero entry
"""
struct UnitVector{T} <: AbstractVector{T}
    i::Int
    n::Int
    v::T
end

function Base.getindex(e::UnitVector{T}, i::Int) where T
    @boundscheck @assert 1 <= i <= e.n
    return i == e.i ? e.v : zero(T)
end

Base.size(e::UnitVector) = (e.n,)

function Base.:(*)(A::AbstractMatrix{N}, e::UnitVector{N}) where {N}
    return A[:, e.i] * e.v
end

function Base.:(*)(A::Transpose{N, <:AbstractMatrix{N}}, e::UnitVector{N}
                  ) where {N}
    return A[:, e.i] * e.v
end

# diagonal matrix times unit vector
function Base.:(*)(D::Diagonal{N, V}, e::UnitVector{N}) where {N, V<:AbstractVector{N}}
    return UnitVector(e.i, e.n, D.diag[e.i] * e.v)
end

function inner(e1::UnitVector{N}, A::AbstractMatrix{N}, e2::UnitVector{N}
              ) where {N}
    return A[e1.i, e2.i] * e1.v * e2.v
end

"""
    BoxDirections{N} <: AbstractDirections{N}

Box direction representation.

### Fields

- `n` -- dimension
"""
struct BoxDirections{N} <: AbstractDirections{N}
    n::Int
end

# constructor for type Float64
BoxDirections(n::Int) = BoxDirections{Float64}(n)

Base.eltype(::Type{BoxDirections{N}}) where {N} = AbstractVector{N}
Base.length(bd::BoxDirections) = 2 * bd.n

function Base.iterate(bd::BoxDirections{N}, state::Int=1) where {N}
    if state == 0
        return nothing
    end
    vec = UnitVector{N}(abs(state), bd.n, convert(N, sign(state)))
    state = (state == bd.n) ? -bd.n : state + 1
    return (vec, state)
end

# ==================================================
# Octagonal directions
# ==================================================

"""
    OctDirections{N} <: AbstractDirections{N}

Octagon direction representation.

### Fields

- `n` -- dimension

### Notes

Octagon directions consist of all vectors that are zero almost everywhere except
in two dimensions ``i``, ``j`` (possibly ``i = j``) where it is ``±1``.
"""
struct OctDirections{N} <: AbstractDirections{N}
    n::Int
end

# constructor for type Float64
OctDirections(n::Int) = OctDirections{Float64}(n)

Base.eltype(::Type{OctDirections{N}}) where {N} = AbstractVector{N}
Base.length(od::OctDirections) = 2 * od.n^2

function Base.iterate(od::OctDirections{N}) where {N}
    if od.n == 1
        # fall back to box directions in 1D case
        return iterate(od, 1)
    end
    vec = spzeros(N, od.n)
    vec[1] = one(N)
    vec[2] = one(N)
    return (copy(vec), (vec, 1, 2))
end

# basic idea: modify the vector from the previous iteration
# have two indices i = 1, j = i + 1
# - run j through all indices > i
# - i = i + 1, j = i + 1
# - for any pair (i, j) create four vectors [..., i: ±1, ..., j: ±1, ...]
# in the end continue with box directions
function Base.iterate(od::OctDirections{N}, state::Tuple) where {N}
    # continue with octagon directions
    vec = state[1]
    i = state[2]
    j = state[3]
    if vec[j] > 0
        # change sign at j
        vec[j] = -one(N)
    elseif vec[i] > 0
        # change sign at i and j
        vec[i] = -one(N)
        vec[j] = one(N)
    elseif j < od.n
        # advance j, reset i
        vec[i] = one(N)
        vec[j] = zero(N)
        j = j + 1
        vec[j] = one(N)
    elseif i < od.n - 1
        # advance i
        vec[i] = zero(N)
        vec[j] = zero(N)
        i = i + 1
        j = i + 1
        vec[i] = one(N)
        vec[j] = one(N)
    else
        # continue with box directions
        return iterate(od, 1)
    end
    return (copy(vec), (vec, i, j))
end

function Base.iterate(od::OctDirections{N}, state::Int) where {N}
    # continue with box directions
    return iterate(BoxDirections{N}(od.n), state)
end

# ==================================================
# Box-diagonal directions
# ==================================================

"""
    BoxDiagDirections{N} <: AbstractDirections{N}

Box-diagonal direction representation.

### Fields

- `n` -- dimension

### Notes

Box-diagonal directions can be seen as the union of diagonal directions (all
entries are ±1) and box directions (one entry is ±1, all other entries are 0).
The iterator first enumerates all diagonal directions, and then all box
directions.
"""
struct BoxDiagDirections{N} <: AbstractDirections{N}
    n::Int
end

# constructor for type Float64
BoxDiagDirections(n::Int) = BoxDiagDirections{Float64}(n)

Base.eltype(::Type{BoxDiagDirections{N}}) where {N} = AbstractVector{N}
Base.length(bdd::BoxDiagDirections) = bdd.n == 1 ? 2 : 2^bdd.n + 2 * bdd.n

Base.iterate(bdd::BoxDiagDirections{N}) where {N} =
    (ones(N, bdd.n), ones(N, bdd.n))

function Base.iterate(bdd::BoxDiagDirections{N},
                      state::AbstractVector) where {N}
    # continue with diagonal directions
    i = 1
    while i <= bdd.n && state[i] < 0
        state[i] = -state[i]
        i = i+1
    end
    if i > bdd.n
        if bdd.n == 1
            # finish here to avoid duplicates
            return nothing
        else
            # continue with box directions
            return iterate(bdd, 1)
        end
    else
        state[i] = -state[i]
        return (copy(state), state)
    end
end

function Base.iterate(bdd::BoxDiagDirections{N}, state::Int) where {N}
    # continue with box directions
    return iterate(BoxDirections{N}(bdd.n), state)
end

# ==================================================
# Polar directions
# ==================================================

"""
    PolarDirections{N<:AbstractFloat} <: AbstractDirections{N}

Polar directions representation.

### Fields

- `Nφ`    -- length of the partition of the polar angle

### Notes

The `PolarDirections` constructor provides a sample of the unit sphere
in ``\\mathbb{R}^2``, which is parameterized by the polar angles
``φ ∈ Dφ := [0, 2π]`` respectively; see the wikipedia entry
[Polar coordinate system](https://en.wikipedia.org/wiki/Polar_coordinate_system).
The domain ``Dφ`` is discretized in ``Nφ`` pieces.
Then the Cartesian components of each direction are obtained with

```math
[cos(φᵢ), sin(φᵢ)].
```

### Examples

The integer passed as an argument is used to discretize ``φ``:

```jldoctest spherical_directions
julia> using LazySets.Approximations: PolarDirections

julia> pd = PolarDirections(2)
PolarDirections{Float64}(2, Array{Float64,1}[[1.0, 0.0], [-1.0, 1.22465e-16]])

julia> pd.Nφ
2
```
"""
struct PolarDirections{N<:AbstractFloat} <: AbstractDirections{N}
    Nφ::Int
    stack::Vector{Vector{N}} # stores the polar directions

    function PolarDirections{N}(Nφ::Int) where {N<:AbstractFloat}
        if Nφ <= 0
            throw(ArgumentError("Nφ = $Nφ is invalid; it shoud be at least 1"))
        end
        stack = Vector{Vector{N}}()
        # discretization of the polar angle
        φ = range(N(0.0), N(2*pi), length=Nφ+1)

        for φᵢ in φ[1:Nφ]  # skip last (repeated) angle
            d = N[cos(φᵢ), sin(φᵢ)]
            push!(stack, d)
        end
        return new{N}(Nφ, stack)
    end
end

PolarDirections(Nφ::Int) = PolarDirections{Float64}(Nφ)

# common functions
Base.eltype(::Type{PolarDirections{N}}) where {N} = Vector{N}
Base.length(pd::PolarDirections) = length(pd.stack)
dim(::PolarDirections) = 2

function Base.iterate(pd::PolarDirections, state::Int=1)
    state == length(pd.stack)+1 && return nothing
    return (pd.stack[state], state + 1)
end

# ==================================================
# Spherical directions
# ==================================================

"""
    SphericalDirections{N<:AbstractFloat} <: AbstractDirections{N}

Spherical directions representation.

### Fields

- `Nθ`    -- length of the partition of the azimuthal angle
- `Nφ`    -- length of the partition of the polar angle
- `stack` -- list of computed directions

### Notes

The `SphericalDirections` constructor provides a sample of the unit sphere
in ``\\mathbb{R}^3``, which is parameterized by the azimuthal and polar angles
``θ ∈ Dθ := [0, π]`` and ``φ ∈ Dφ := [0, 2π]`` respectively, see the wikipedia
entry [Spherical coordinate system](https://en.wikipedia.org/wiki/Spherical_coordinate_system).
The domains ``Dθ`` and ``Dφ`` are discretized in ``Nθ`` and ``Nφ`` respectively.
Then the Cartesian componentes of each direction are obtained with

```math
[sin(θᵢ)*cos(φᵢ), sin(θᵢ)*sin(φᵢ), cos(θᵢ)].
```

The north and south poles are treated separately so that those points
are not considered more than once.

### Examples

A `SphericalDirections` can be built in different ways. If you pass only one integer,
it is used to discretize both ``θ`` and ``φ``:

```jldoctest spherical_directions
julia> using LazySets.Approximations: SphericalDirections

julia> sd = SphericalDirections(3)
SphericalDirections{Float64}(3, 3, Array{Float64,1}[[0.0, 0.0, 1.0], [0.0, 0.0, -1.0], [1.0, 0.0, 6.12323e-17], [-1.0, 1.22465e-16, 6.12323e-17]])

julia> sd.Nθ, sd.Nφ
(3, 3)
```

Pass two integers to control the discretization in ``θ`` and in ``φ`` separately:

```jldoctest spherical_directions
julia> sd_4_5 = SphericalDirections(4, 5);

julia> length(sd_4_5)
10

julia> sd_4_8 = SphericalDirections(4, 8);

julia> length(sd_4_8)
16
```
"""
struct SphericalDirections{N<:AbstractFloat} <: AbstractDirections{N}
    Nθ::Int
    Nφ::Int
    stack::Vector{Vector{N}} # stores the spherical directions

    function SphericalDirections{N}(Nθ::Int, Nφ::Int) where {N<:AbstractFloat}
        if Nθ <= 1 || Nφ <= 1
            throw(ArgumentError("(Nθ, Nφ) = ($Nθ, $Nφ) is invalid; both shoud be at least 2"))
        end
        stack = Vector{Vector{N}}()
        θ = range(N(0), N(pi), length=Nθ)      # discretization of the azimuthal angle
        φ = range(N(0.0), N(2*pi), length=Nφ)  # discretization of the polar angle

        # add north pole (θ = 0)
        push!(stack, N[0, 0, 1])

        # add south pole (θ = pi)
        push!(stack, N[0, 0, -1])

        for φᵢ in φ[1:Nφ-1]  # delete repeated angle
            for θⱼ in θ[2:Nθ-1] # delete north and south poles
                d = N[sin(θⱼ)*cos(φᵢ), sin(θⱼ)*sin(φᵢ), cos(θⱼ)]
                push!(stack, d)
            end
        end
        return new{N}(Nθ, Nφ, stack)
    end
end

# convenience constructors
SphericalDirections(Nθ::Int) = SphericalDirections{Float64}(Nθ, Nθ)
SphericalDirections(Nθ::Int, Nφ::Int) = SphericalDirections{Float64}(Nθ, Nφ)

# common functions
Base.eltype(::Type{SphericalDirections{N}}) where {N} = Vector{N}
Base.length(sd::SphericalDirections) = length(sd.stack)
dim(::SphericalDirections) = 3

function Base.iterate(sd::SphericalDirections, state::Int=1)
    state == length(sd.stack)+1 && return nothing
    return (sd.stack[state], state + 1)
end


text:
@testset "Overorders" begin
  Qx,  x = FlintQQ["x"]

  f = x^2 + 100x + 100

  A = AlgAss(f * (f - 1))
  O = Order(A, basis(A))

  orders = @inferred poverorders(O, 2)
  @test length(orders) == 6

  orders = @inferred poverorders(O, fmpz(5))
  @test length(orders) == 2

  orders = @inferred poverorders(O, fmpz(7))
  @test length(orders) == 3

  orders = @inferred overorders(O)
  @test length(orders) == 36

  @test count(is_gorenstein, orders) == 36
  @test count(is_bass, orders) == 36
  @test count(is_maximal, orders) == 1

  orders = @inferred overorders(O, type = :bass)
  @test length(orders) == 36

  orders = @inferred overorders(O, type = :gorenstein)
  @test length(orders) == 36

  G = small_group(6, 1)
  A = AlgAss(group_algebra(FlintQQ, G))[1]
  O = Order(A, basis(A))
  orders = @inferred overorders(O)
  @test length(orders) == 12

  orders = @inferred overorders(O, type = :all)
  @test length(orders) == 12

  @test count(is_maximal, orders) == 2

  @test_throws ErrorException overorders(O, type =  :gorenstein)
  @test_throws ErrorException overorders(O, type =  :bass)
  @test_throws ErrorException overorders(O, type =  :blub)

  f = x^4-1680*x^3-25920*x^2-1175040*x+25214976
  O = EquationOrder(f)
  orders = @inferred overorders(O)
  @test length(orders) == 2535
  @test count(is_maximal, orders) == 1
  @test count(is_gorenstein, orders) == 657
  @test count(is_bass, orders) == 5

  @test length(overorders(O, type = :gorenstein)) == 657
  @test length(overorders(O, type = :bass)) == 5

  f = x^4-1680*x^3-25920*x^2-1175040*x+25214976
  A = AlgAss(f)
  O = Order(A, basis(A))
  @test length(overorders(O)) == 2535
end


text:
push!(LOAD_PATH, pwd())
include("startup.jl")

#using PyCall
#pygui(:tk) # :tk, :gtk3, :gtk, :qt5, :qt4, :qt, or :wx
import PyPlot
import .InterpolationFunctions
import  .MasslistFunctions
import .ResultFileFunctions

file = "ExampleFiles/TOFDATA/results/_result.hdf5"
availableTraces = ResultFileFunctions.loadResults(file,masslistOnly=true)
linewidth = 2

fig=PyPlot.figure()
concAx = PyPlot.subplot(1,1,1)
data =   ResultFileFunctions.getTraces(file,massIndices=[1])
concAx.semilogy(availableTraces.Times, data, linewidth=linewidth, label="m/z $(round(availableTraces.MasslistMasses[1],digits=3)) -  $(MasslistFunctions.sumFormulaStringFromCompositionArray(availableTraces.MasslistCompositions[:,1]))","b")
concAx.set_ylabel("CPS")
concAx.set_xlabel("Time")
concAx.grid()
concAx.legend(loc=1)
lastPlottedIndex = 1
activeIndices = []

function changeLastPlotTo(i)
    global lastPlottedIndex, data
    if length(concAx.lines) > 0
        concAx.lines[end].remove()
    end
    data = ResultFileFunctions.getTraces(file,massIndices=[i])
    concAx.semilogy(availableTraces.Times, data, linewidth=linewidth, label="m/z $(round(availableTraces.MasslistMasses[i],digits=3)) -  $(MasslistFunctions.sumFormulaStringFromCompositionArray(availableTraces.MasslistCompositions[:,i]))","b")
    concAx.legend(loc=1)
    global lastPlottedIndex = i
    println("Setting $lastPlottedIndex as last plot")
end

function addNewPlot()
    global lastPlottedIndex, data
    println("Adding $lastPlottedIndex to plots")
    if length(concAx.lines) > 0
        concAx.lines[end].remove()
    end
    data = ResultFileFunctions.getTraces(file,massIndices=[lastPlottedIndex])
    concAx.semilogy(availableTraces.Times, data, linewidth=linewidth, label="m/z $(round(availableTraces.MasslistMasses[lastPlottedIndex],digits=3)) -  $(MasslistFunctions.sumFormulaStringFromCompositionArray(availableTraces.MasslistCompositions[:,lastPlottedIndex]))")
    concAx.semilogy(availableTraces.Times, data, linewidth=linewidth, label="m/z $(round(availableTraces.MasslistMasses[lastPlottedIndex],digits=3)) -  $(MasslistFunctions.sumFormulaStringFromCompositionArray(availableTraces.MasslistCompositions[:,lastPlottedIndex]))","b")
    concAx.legend(loc=1)
end

function keyReleaseHandler(event)
    global activeIndices
    k=event.key
    if k == " "
        println("Adding Plot permanently")
        addNewPlot()
        push!(activeIndices,lastPlottedIndex)
    end
end
eventStepCumulator=0
function scrollHandler(event)
    global lastPlottedIndex, eventStepCumulator
    if event.step != 0
        println("Scrolled $(event[:step])")
    end
    eventStepCumulator += event.step # mouse pads scroll fractional
    if (eventStepCumulator >= 1) || (eventStepCumulator <=-0.6)
        plotIndex = clamp(lastPlottedIndex + eventStepCumulator, 1, length(availableTraces.MasslistMasses))
        changeLastPlotTo(Int64(floor(plotIndex)))
        eventStepCumulator = 0
    end
end

fig.canvas.mpl_connect("scroll_event", scrollHandler)
fig.canvas.mpl_connect("key_release_event", keyReleaseHandler)


text:
function histopryremapp(t,A::Matrix,SolutionSetV,solindex)
    #t [-τ,0]
    sum([SolutionSetV[i](t) * A[i,solindex] for i in 1:length(SolutionSetV)])
end

#function extract_solutionpoints(sol,time)
#    [solset[solind](t) for t in timelocations, solind in 1:eigvecnum];
#end


text:
# youtube-dl backend mode for unsupported video
function youtube_dl(v_id, opts)
    @info("Entered youtube-dl MODE")
    if !isavailable(YTDLbin)
        return @error("""

            Command `youtube-dl` is unavailable on your computer.
            Install it and try again!
        """)
    end
    u = "https://wwww.youtube.com/watch?v=" * v_id
    codecs = youtube_dl_available_codecs(u)
    c = sellect_format(codecs, opts)
    ext = sellect_ext(c.a, c.v)
    base = (
        "youtube-dl",
        "-o",
        joinpath(opts.dir, "%(title)s.%(ext)s"),
        "-f",
        isnothing(c.v) ? c.a : string(c.v, "+", c.a)
    )
    if isnothing(c.v)
        fmt = ("--extract-audio", "--audio-format", ext.a)
    else
        fmt = ("--merge-output-format", ext.v)
    end
    if haskey(opts, :backend_options)
        backend_opts = opts[:backend_options]
    else
        backend_opts = [
            "--no-continue",
            "--no-mtime",
            "-q",
            "--console-title",
            "--add-metadata"
        ]
    end
    t = Sys.get_process_title()
    run(`$base $backend_opts $fmt $u`)
    Sys.set_process_title(t)
    return nothing
end

function youtube_dl_available_codecs(url)
    lines = readlines(`youtube-dl -F $url`)
    Int[
        parse(Int, first(l, 3))
            for l ∈ lines
                if isdigit(l[1]) && isdigit(l[3])
    ]
end


text:
# when calculating a density product, there is a certain element selection of incoming elements
# wont yet allow too much hybridization in types but enough to get going
struct DensityProductElements{D,Vo<:AbstractVector,B,Vi<:AbstractVector}
  outElements::Vo
  outBW::Vector{B}
  outElemName::Symbol
  outdatedBW::Base.RefValue{Bool}
  lblCombinations::Vector{Vector{Int64}}
  inElements::NTuple{D,Vi}
  inBWs::NTuple{D,<:AbstractVector{B}}
  inElemNames::NTuple{D,Symbol}
  inElemFactorNames::NTuple{D,Symbol}
  # Pair{inDens, elemIdx}
  _sequencedUpdates::Channel{Pair{Int,Int}}
end

function Base.show(io::IO, x::DensityProductElements)
  # println(io, )
  printstyled(io, "DensityProductElements:", color=:blue)
  println(io)
  println(io, "  Out Name:     ", x.outElemName )
  println(io, "   in names:    ", x.inElemNames)
  println(io, "   # densities: ", length(x.inElements))
  println(io, "  # dims:       ", length(x.inElements[1][1]))
  println(io, "  # out pts:    ", length(x.outElements))
  println(io, "  .outBW[1]~:   ", x.outBW[1] .|> x->round(x,digits=5) )
  println(io, "  .outdatedBW:  ", x.outdatedBW[] )
  # println(io)
end

Base.show(io::IO, ::MIME"text/plain", x::DensityProductElements) = show(io, x)

# assume `mdp.inElements` is already set and corroborates with `lbsChs`
function _setLabelCombinations!(mdp::DensityProductElements,
                                lbsChs::Dict )
  #
  npts_ = length(lbsChs)
  ndens = length(lbsChs[1])

  resize!(mdp.outElements, npts_)
  resize!(mdp.lblCombinations, npts_)

  NptsD = length.(mdp.inElements)

  # how many points i in the product belief
  for i in 1:npts_
    # make sure memory allocations are correct
    resize!(mdp.lblCombinations[i], ndens)
    # how many densities to make the product
    for densid in 1:ndens
      # level in the kd tree
      level = length(lbsChs[i][densid])
      # get label selection for densid
      lbl = lbsChs[i][densid][level] - NptsD[densid]
      # which row in xarrTiles to mark
      mdp.lblCombinations[i][densid] = lbl
      # xarrTiles[densid][lbl,i] = 1
    end
  end

  nothing
end

# # can only do for Array, not view
# function _setProductElements!(mdp::DensityProductElements{D},
#                               prd::BallTreeDensity)
#   #
#   # also set the bandwidth
#   dim = Ndim(prd)
#   resize!(mdp.outBW, dim)

#   npts = Npts(prd)
#   for i in 1:
#     mdp.outBW[:,i] .= getBW(prd)[:,1] # fix for all elements
#   end

#   # set kernel center elements
#   resize!(mdp.outElements, )
#   for i in 1:Npts(prd)
#     resize!(mdp.outElements[i], dim)
#     mdp.outElements[i][:] .= getPoints(prd, i)
#   end

#   #
#   nothing
# end

##

"""
Likely new bug see KDE #70
"""
function _buildDensityProductElements(XX::AbstractVector{B};
                                      outName::Symbol=:product,
                                      inNames::Union{<:AbstractVector{Symbol},NTuple{D,Symbol}}=[Symbol("belief$i") for i in 1:length(XX)],
                                      inFctNames::Union{<:AbstractVector{Symbol},NTuple{D,Symbol}}=[Symbol("factor$i") for i in 1:length(XX)],
                                      _glbs=KDE.makeEmptyGbGlb(recordChoosen=true),
                                      product::B = *(XX, glbs=_glbs, addEntropy=false ) ) where {B<: BallTreeDensity, D}
  #

  npts = Npts(product)
  ndim = Ndim(product)
  ndens = length(XX)

  # pts = getPoints(product)
  # @cast outArr[j][i] := pts[i,j]
  outArr = [getPoints(product,i) for i in 1:npts]
  bw = [getBW(product)[:,i] for i in 1:npts]
  lblComb = [zeros(Int,ndens) for i in 1:npts]

  # restructure points for incoming densities
  # indens = (XX...,)
  XXarr = [([getPoints(x, i) for i in 1:Npts(x)]) for x in XX]
  BWarr = [([getBW(x)[:,i] for i in 1:Npts(x)]) for x in XX]

  # build the object
  mdp = DensityProductElements(outArr, bw, outName, Ref(true), lblComb, (XXarr...,), (BWarr...,), (inNames...,), (inFctNames...,), Channel{Pair{Int,Int}}(100))

  # set the labels selections used for current product
  if length(XX) == 1
    resize!(mdp.lblCombinations, npts)
    for i in 1:npts
      resize!(mdp.lblCombinations[i],1)
      mdp.lblCombinations[i][1] = i
    end
  else
    _setLabelCombinations!(mdp, _glbs.labelsChoosen)
  end

  return mdp
end

function _recalcProductKernel(mdp::DensityProductElements, idx::Int)
  #
  dim = length(mdp.inElements[1][1])
  selections = mdp.lblCombinations[idx]

  # TODO still have to do per dimension...
  infor = zeros(dim)
  informn = zeros(dim)
  for (idd, lbid) in enumerate(selections)
    # dens = mdp.inElements[idd]
    # mn_i = getPoints(dens, lbid)
    mn_i = mdp.inElements[idd][lbid]
    recpcov = 1/mdp.inBWs[idd][lbid].^2
    infor[:] .+= recpcov[:] # NOTE currently only diagonals
    informn += recpcov.*mn_i
  end
  cov = 1 ./ infor
  mn = informn./infor
  return mn, cov
end

##

function _updateProposalElement(mdp::DensityProductElements, elemval, densid::Int, elemidx::Int)
  # won't work for non-vector elements
  mdp.inElements[densid][elemidx] .= elemval
  mdp.outdatedBW[] = false
  nothing
end

function _updateOutElement(dpe::DensityProductElements, elemval, elemidx::Int)
  # won't work for non-vector elements
  dpe.outElements[elemidx] .= elemval
  dpe.outdatedBW[] = false
  nothing
end

function _sequenceUpdateProposalElement(mdp::DensityProductElements, elemval, densid::Int, elemidx::Int)
  _updateProposalElement(mdp, elemval, densid, elemidx)
  newpair = densid=>elemidx
  put!(mdp._sequencedUpdates, newpair)
  nothing
end

function _sequenceUpdateOutElement(dpe::DensityProductElements, elemval, elemidx::Int)
  _updateOutElement(dpe, elemval, elemidx)
  @warn("WIP")
end

_listOutElementLabelSelections(mdp::DensityProductElements, idx::Int) = mdp.lblCombinations[idx]

function _findOutElementConnectedLabel(mdp::DensityProductElements, densidx, idx::Int)
  outelems = Int[]
  # check for all elements
  for outi in 1:length(mdp.outElements)
    # if there happens to be a hit (out product relies on an in label selection)
    if _listOutElementLabelSelections(mdp, idx)[densidx] == outi
      push!(outelems, outi)
    end
  end

  # return list of outElements that depend on one particular inbound kernel `densidx=>idx`
  return outelems
end

# if something in the queue effects an outElement, then update values.
function _takeOutElementFromSequence(mdp::DensityProductElements)
  if !isready(mdp._sequencedUpdates)
    return nothing
  end
  nextpair = take!(mdp._sequencedUpdates)

  outElemList = _findOutElementConnectedLabel(mdp, nextpair[1], nextpair[2])
  println("Connected label $outElemList")

  return nothing
end


text:
# define IC function, BC functor here
function TestIC(mesh, sbp, eqn, opts, q_vec)
  println("hello world")
end

mutable struct TestBCType <: BCType
end

function (obj::TestBCType)()
  println("hello edge of the world")
end

function test_frontend()

  @testset "----- Testing PDESolver Module frontend -----" begin
    physics_name = "Advection"
    @test ( haskey(PDESolver.PhysicsModDict, physics_name) )== true
    mod, _createObjects, _checkOptions = retrieve_physics(physics_name)
    # FactCheck doesn't compare functions correctly, so use == instead
    @test ( mod == AdvectionEquationMod )== true
    @test ( _createObjects == AdvectionEquationMod.createObjects )== true
    @test ( _checkOptions == AdvectionEquationMod.checkOptions )== true

    # test run_solver()
    fname = "input_vals_channel.jl"
    mesh, sbp, eqn, opts = solvePDE(fname)
    mesh, sbp, eqn2, opts = run_solver(fname)

    @test isapprox( eqn.q_vec, eqn2.q_vec) atol=1e-13

    ic_name = "TestIC"
    registerIC(AdvectionEquationMod, ic_name, TestIC)
    @test ( haskey(AdvectionEquationMod.ICDict, ic_name) )== true
    @test ( AdvectionEquationMod.ICDict[ic_name] == TestIC )== true

    bc_name = "TestBC"
    registerBC(AdvectionEquationMod, bc_name, TestBCType())
    @test ( haskey(AdvectionEquationMod.BCDict, bc_name) )== true
    @test ( AdvectionEquationMod.BCDict[bc_name] )== TestBCType()

  end  # end facts block

  return nothing

end

add_func1!(AdvectionTests, test_frontend, [TAG_FRONTEND, TAG_SHORTTEST])


text:
using Combinatorics

function findMinWrapDist(f_true,f_hat)
# Assign each of the ground truth frequencies to closest
# (squared error) canditate of estimated frequencies.
# Metric is the wrap-aroud distance on the torus $\mathds{T}=[0,1]$.
#
# Inputs:
# f_true: ground truth frequencies
# f_hat: estimated frequencies
#
# Output:
# dists: wrap-around distances of the frequencies in f_true and f_hat

K = length(f_true)
WrapDist(x,y) = min(abs(x-y),abs(x-y+1))

dists = Array{Float64,1}(undef,K)
for ii = 1:K
    dists[ii] = minimum(WrapDist.(f_true[ii],f_hat))
end

return dists
end


text:
using FractionalDiffEq, Plots

fun(x) = 0.5*x+1
α=0.5;x0=1;
T=1; h=0.1
prob = FractionalDifferenceProblem(fun, α, x0)
t, y=solve(prob, T, h, PECEDifference())
plot(t, y, seriestype=:scatter, legend=:bottomright)


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

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

                     @variable(m, -23.82106505746807 <= q <= 20.87497636142748)

                     add_NL_constraint(m, :(softplus(-0.7432570628691102 + 0.9028765283094948*softplus(-0.11514299287103702 + 0.5608817277242992*$(x[1]) + 0.6700218006141032*$(x[2]) + -0.9693868545458724*$(x[3])) + -0.9177729976456397*softplus(-0.5568784315984976 + 0.281764667612443*$(x[1]) + 0.42150753955170206*$(x[2]) + -0.5829321104114529*$(x[3])) + -0.6091404444372248*softplus(-0.04388549286647825 + -0.9010617152865183*$(x[1]) + -0.7540101668519821*$(x[2]) + -0.4758905699629783*$(x[3])) + 0.7075296204914552*softplus(-0.511868403698823 + 0.5804725192481515*$(x[1]) + 0.7124288843879394*$(x[2]) + -0.06834010138653612*$(x[3])) + -0.9115635329002854*softplus(-0.9281357966735606 + -0.33421670762637223*$(x[1]) + -0.5412233562926474*$(x[2]) + -0.2936515077150612*$(x[3]))) + softplus(-0.2671327260424561 + -0.018499351907995898*softplus(-0.11514299287103702 + 0.5608817277242992*$(x[1]) + 0.6700218006141032*$(x[2]) + -0.9693868545458724*$(x[3])) + 0.12529725503678657*softplus(-0.5568784315984976 + 0.281764667612443*$(x[1]) + 0.42150753955170206*$(x[2]) + -0.5829321104114529*$(x[3])) + 0.9839396030143157*softplus(-0.04388549286647825 + -0.9010617152865183*$(x[1]) + -0.7540101668519821*$(x[2]) + -0.4758905699629783*$(x[3])) + 0.28500560463058955*softplus(-0.511868403698823 + 0.5804725192481515*$(x[1]) + 0.7124288843879394*$(x[2]) + -0.06834010138653612*$(x[3])) + 0.970334521625404*softplus(-0.9281357966735606 + -0.33421670762637223*$(x[1]) + -0.5412233562926474*$(x[2]) + -0.2936515077150612*$(x[3]))) + softplus(0.40040910594674495 + -0.7454849743358074*softplus(-0.11514299287103702 + 0.5608817277242992*$(x[1]) + 0.6700218006141032*$(x[2]) + -0.9693868545458724*$(x[3])) + 0.7108986625228169*softplus(-0.5568784315984976 + 0.281764667612443*$(x[1]) + 0.42150753955170206*$(x[2]) + -0.5829321104114529*$(x[3])) + -0.30349624151173904*softplus(-0.04388549286647825 + -0.9010617152865183*$(x[1]) + -0.7540101668519821*$(x[2]) + -0.4758905699629783*$(x[3])) + -0.26057109207128715*softplus(-0.511868403698823 + 0.5804725192481515*$(x[1]) + 0.7124288843879394*$(x[2]) + -0.06834010138653612*$(x[3])) + 0.7858005505148937*softplus(-0.9281357966735606 + -0.33421670762637223*$(x[1]) + -0.5412233562926474*$(x[2]) + -0.2936515077150612*$(x[3]))) + softplus(-0.7120753920031948 + -0.005433419997351141*softplus(-0.11514299287103702 + 0.5608817277242992*$(x[1]) + 0.6700218006141032*$(x[2]) + -0.9693868545458724*$(x[3])) + -0.8839514190404207*softplus(-0.5568784315984976 + 0.281764667612443*$(x[1]) + 0.42150753955170206*$(x[2]) + -0.5829321104114529*$(x[3])) + -0.2933551592941672*softplus(-0.04388549286647825 + -0.9010617152865183*$(x[1]) + -0.7540101668519821*$(x[2]) + -0.4758905699629783*$(x[3])) + -0.5418193939730438*softplus(-0.511868403698823 + 0.5804725192481515*$(x[1]) + 0.7124288843879394*$(x[2]) + -0.06834010138653612*$(x[3])) + 0.8978172795430157*softplus(-0.9281357966735606 + -0.33421670762637223*$(x[1]) + -0.5412233562926474*$(x[2]) + -0.2936515077150612*$(x[3]))) + softplus(0.7950909223530527 + -0.6380481687930124*softplus(-0.11514299287103702 + 0.5608817277242992*$(x[1]) + 0.6700218006141032*$(x[2]) + -0.9693868545458724*$(x[3])) + -0.8467501205760102*softplus(-0.5568784315984976 + 0.281764667612443*$(x[1]) + 0.42150753955170206*$(x[2]) + -0.5829321104114529*$(x[3])) + 0.5001198683847332*softplus(-0.04388549286647825 + -0.9010617152865183*$(x[1]) + -0.7540101668519821*$(x[2]) + -0.4758905699629783*$(x[3])) + -0.03855272462979853*softplus(-0.511868403698823 + 0.5804725192481515*$(x[1]) + 0.7124288843879394*$(x[2]) + -0.06834010138653612*$(x[3])) + 0.3301524112226937*softplus(-0.9281357966735606 + -0.33421670762637223*$(x[1]) + -0.5412233562926474*$(x[2]) + -0.2936515077150612*$(x[3]))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
function is_prime(i)
    for j=3:2:round(Int, sqrt(i+1))
        i % j == 0 && return 0
    end
    return 1
end

function count_primes(n)
    cnt = 1
    for i = 3:2:n
        cnt += is_prime(i)
    end
    return cnt
end

count_primes(10)
n = parse(Int, ARGS[1])
t0 = time()
v = count_primes(n)
dt = time() - t0
println("count_primes($n) = $v (run time $dt seconds)")


text:
using MinFinVol

include("DTypeTests.jl")
include("OperationTests.jl")


text:
function chooseFactorStructure(X :: Matrix, channel_layout :: MultiChannelLayout)
    n = size(X)[1]
    T = size(X)[2]
    XtX = 1/(n-1)*transpose(X) * X
    S   = 1/(n-1)*X * transpose(X)
    overall_evs = reverse(eigen(XtX.values, dims=1))
    X_by_channel = extract_horizontal_blocks(X, channel_layout.nobs_per_channel)
    channel_evs = [reverse(eigen(transpose(Xc) * Xc).values) for Xc in X_by_channel]
    D = blockdiag(extract_diagonal_blocks(S,   channel_layout.nobs_per_channel, channel_layout.nobs_per_channel))
    U = S - D

    off_diag_evs = reverse(eigen(XtX.values, dims=1)[T-n+1:T])
end

function Hstep(S :: Symmetric, params :: MCFMParams)
    U = params.factorlayout.nchannelspecific # total num of unique factors
    L = params.channellayout.nobs # total number of obs, across all channels.
    G = params.G
    Σ = params.Σ
    p = params.factorlayout.nchannelshared

    E = G * transpose(G) + Σ # eqn 24
    Eisqrt = sqrt(inv(E)) # inefficient
    Esqrt  = sqrt(E)

    Swhite = Eisqrt * S * Eisqrt # eqn 25
    Swhite_eigen = svd(Swhite)
    W = Swhite_eigen.U

    #TODO: Notation in paper is ambigious, double check defn of D.
    ds = max.(Swhite_eigen.S[1:p]  .- 1, zeros(p))
    Dsqrt = vcat(sqrt(diagm(ds)) , zeros((L-p, p)))

    B = Esqrt * W * Dsqrt # before eqn 11
    Q = LTOrthog(B, p)

    Hnew = B * Q
    Hnew = Array{Float64, 2}(Hnew)
    return Hnew
end

function _Gstep_channel(S, H, Σ, pj)
    Lj = size(H)[1]
    P = S - H * transpose(H) # eqn 44
    Σsqrt = sqrt(Σ)
    Σisqrt = sqrt(inv(Σ))
    Pwhite = Σisqrt * P * Σisqrt
    Pwhite_eigen = svd(Pwhite)
    W = Pwhite_eigen.U

    ds = max.(Pwhite_eigen.S[1:pj] .-1, zeros(pj))
    Dsqrt = vcat(sqrt(diagm(ds)) , zeros((Lj-pj, pj))) #TODO: Still questionable.

    B = Σsqrt * W * Dsqrt
    Q = LTOrthog(B, pj)
    Gnew = B * Q
    return Gnew
end

function Gstep(S :: Symmetric, params :: MCFMParams)
    Sblocks = extract_diagonal_blocks(S, params.channellayout.nobs_per_channel,params.channellayout.nobs_per_channel)
    L = params.channellayout.nobs
    Gblocks = Array{Float64, 2}[]
    for i in 1:params.channellayout.nchannels
        if params.factorlayout.nspecific_per_channel[i] == 0
            # If no channel specific factors, don't run Gstep
            Gblock = zeros(params.channellayout.nobs_per_channel[i], 0)
        else
            Gblock = _Gstep_channel(Sblocks[i],
                params.H_by_channel[i],
                params.Σ_by_channel[i],
                params.factorlayout.nspecific_per_channel[i])
            end
        push!(Gblocks, Gblock)
    end
    return blockdiag(Gblocks)
end

function Σstep(S :: Symmetric, params :: MCFMParams)
    H = params.H
    G = params.G
    Est = S - H * transpose(H) - G * transpose(G)
    Σnew = Matrix(Diagonal(Est)) # extract only the diagonal elems.
    return Σnew
end

function fit(data :: MultiChannelData, factorlayout :: MultiChannelFactorLayout; tol=1e-4, maxiter=1e6, center=true)
    # Data is assumed to be in a ragged 3d array structure: sample x channel x observation
    X = stack(data)
    if center
        X = X .- mean(X, dims=2)
    end
    # Initialize parameters
    H =  zeros(data.channellayout.nobs, factorlayout.nchannelshared)
    G =  zeros(data.channellayout.nobs, factorlayout.nchannelspecific)
    Σ =  Matrix(I(data.channellayout.nobs))
    params = MCFMParams(factorlayout, data.channellayout, H, G, Σ)
    fitpath = MCFMHistory(1, maxiter, tol, [params])

    # Compute sample cov
    S = Symmetric(1/(data.nsamples-1) * X * transpose(X))

    iter = 2
    converged = false
    while iter < maxiter && !converged
        H = Hstep(S, params)
        params = MCFMParams(factorlayout, data.channellayout, H, G, Σ)
        G = Gstep(S, params)
        params = MCFMParams(factorlayout, data.channellayout, H, G, Σ)
        Σ = Σstep(S, params)
        params = MCFMParams(factorlayout, data.channellayout, H, G, Σ)
        fitpath = MCFMHistory(iter, maxiter, tol, vcat(fitpath.parampath..., [params]))
        if _checkConverged(fitpath, tol)
            _converged = true
            fitpath
        end
        iter = iter + 1
    end
    finalfit = MCFMFit(params, fitpath)
    return finalfit
end

#
# function fitMFA(data, Ls, p, pjs; tol=1e-4, maxiter=1e6)
#
#     X, Xobs, Xchan = stack_and_view(data)
#     # Remove column means
#     #X = X .- mean(X; dims=2) # centering is wrong.
#     # Initialize parameter matrices
#     N = length(data)
#     L = sum(Ls)
#     U = sum(pjs)
#
#     H = zeros(L, p)
#     H_old = zeros(L, p)
#     G = zeros(L, U)
#     G_old = zeros(L, U)
#     Σ = Matrix(I(L))
#     Σ_old = Matrix(I(L))
#
#
#     # Compute sample cov
#     S = Symmetric(1/(N-1) * X * transpose(X))
#
#     #loop vars
#     k = 0
#     converged = false
#     while k < maxiter && !converged
#         H = Hstep(S, G, Σ, p)
#         Hblocks = extract_blocks(H, Ls, [p]; diag=false)
#         Sblocks = extract_blocks(S, Ls, Ls; diag=true)
#         Σblocks = extract_blocks(Σ, Ls, Ls; diag=true)
#         G = Gstep(Sblocks, Hblocks, Σblocks, pjs)
#         Σ = Σstep(S, H, G)
#         if _checkConverged(H_old, G_old, Σ_old, H, G, Σ; tol=tol)
#             converged = true
#             break
#         end
#         k = k + 1
#         H_old = H
#         G_old = G
#         Σ_old = Σ
#     end
#     # Construct MFA object
#     if !converged
#         @warn "Max iteration {maxiter} reached, model has not converged"
#     end
#     return H, G, Σ
# end

function _checkConverged(fitpath :: MCFMHistory, tol)
    iter = fitpath.niter
    H_old = fitpath.parampath[iter-1].H
    G_old = fitpath.parampath[iter-1].G
    Σ_old = fitpath.parampath[iter-1].Σ
    H_new = fitpath.parampath[iter].H
    G_new = fitpath.parampath[iter].G
    Σ_new = fitpath.parampath[iter].Σ
    Hdiffsize = norm(H_new - H_old)
    Gdiffsize = norm(G_new - G_old)
    Σdiffsize = norm(Σ_new - Σ_old)
    Hconv = Hdiffsize/length(H_new) <= tol
    Gconv = Gdiffsize/length(G_new) <= tol
    Σconv = Σdiffsize/size(Σ_old)[1] <= tol # diagonal
    if Hconv && Gconv && Σconv
        return true
    else
        return false
    end
end
#
# function _checkConverged(H_old, G_old, Σ_old, H_new, G_new, Σ_new; tol=1e-6)
#
#     Hdiffsize = norm(H_new - H_old)
#     Gdiffsize = norm(G_new - G_old)
#     Σdiffsize = norm(Σ_new - Σ_old)
#     Hconv = Hdiffsize/length(H_new) <= tol
#     Gconv = Gdiffsize/length(G_new) <= tol
#     Σconv = Σdiffsize/size(Σ_old)[1] <= tol # diagonal
#     if Hconv && Gconv && Σconv
#         return true
#     else
#         return false
#     end
# end

function project(data::MultiChannelData, factors::MultiChannelFactors)
    # Project each factor onto space spanned by common + specific factors
    X = stack(data)
    Xblks = extract_horizontal_blocks(X, data.channellayout.nobs_per_channel)
    Fblks = [hcat(transpose(factors.fcommon), transpose(fu)) for fu in funique_by_channel]
    Xfitteds = Matrix{Float64}(undef, data.channellayout.nchannels)
    for c in 1:data.channellayout.nchannels
        Ft = transpose(F)
        proj = Fblk[c] * pinv(Ft * F) * Ft
        Xblkfit = proj * Xblks[c]
        Xfitteds[c] = Xblkfit
    end
    Xfitted = vcat(Xfitted...)
    out = unstack(Xfitted, data.channellayout)
    return out
end

function extract_factors(data :: MultiChannelData, params ::MCFMParams; center=true)
    H = params.H
    G = params.G
    Σ = params.Σ
    H_by_channel = params.H_by_channel
    G_by_channel = params.G_by_channel
    Σ_by_channel = params.Σ_by_channel
    X = stack(data)
    if center
        X = X .- mean(X, dims=2)
    end
    R = H * transpose(H) + G * transpose(G) + Σ
    Rinv = inv(R)
    fcommon = transpose(H) * inv(R) * X
    funique = zeros(params.factorlayout.nchannelspecific, data.nsamples)
    Xblks = extract_horizontal_blocks(X, data.channellayout.nobs_per_channel)
    j = 0
    for c in 1:data.channellayout.nchannels
        if params.factorlayout.nspecific_per_channel[c] == 0
            continue
        else
            Ri = H_by_channel[c] * transpose(H_by_channel[c]) + G_by_channel[c] * transpose(G_by_channel[c]) + Σ_by_channel[c]
            fu = transpose(G_by_channel[c]) * inv(Ri) * Xblks[c]
            funique[(j+1):(j+size(fu)[1]), :] = fu
            j += size(fu)[1]
        end
    end
    factors = MultiChannelFactors(params.factorlayout, fcommon, funique)
    return factors
end

# function extract_factors(H, G, Σ, X, Ls, pjs)
#     R = H * transpose(H) + G * transpose(G) + Σ
#     #Gblks = extract_blocks(G, Ls; diag=true)
#     #Σblks = extract_blocks(G, Ls; diag=true)
#     #Xblks = extract_blocks(X, Ls; diag=false)
#     Rinv = inv(R)
#     fcommon = transpose(H) * inv(R) * X
#     Hblks = extract_blocks(H, Ls, 1; diag=false)
#     Gblks = extract_blocks(G, Ls, pjs; diag=true)
#     Σblks = extract_blocks(Σ, Ls, Ls; diag=true)
#     Xblks = extract_blocks(X, Ls, 1; diag=false)
#     funique = zeros(sum(pjs), size(X)[2])
#     j=0
#     for (Hi, Gi, Σi, Xi) in zip(Hblks, Gblks, Σblks, Xblks)
#         Ri = Hi * transpose(Hi) + Gi * transpose(Gi) + Σi
#         fu = transpose(Gi) * inv(Ri) * Xi
#         funique[(j+1):(j+size(fu)[1]), :] = fu
#         j += size(fu)[1]
#     end
#     fcommon = transpose(fcommon)
#     funique = transpose(funique)
#     return fcommon, funique
# end

function predict(factors :: MultiChannelFactors, params :: MCFMParams)
    preds = params.loading * factors.factors
    return preds
end
#
# function predict(H, G, Σ, X, Ls)
#     fcom, funi = extract_factors(H, G, Σ, X, Ls)
#     preds = H * (fcom) + G * transpose(funi)
#     return(preds)
# end


text:
struct BoardingPass
    Row::Int
    Column::Int
    SeatID::Int
end

function toBoardingPass(s::String)::BoardingPass
    s = replace(s, "F" => 0)
    s = replace(s, "B" => 1)
    s = replace(s, "L" => 0)
    s = replace(s, "R" => 1)

    row_s = s[1:length(s) - 3]
    row = parse(Int, row_s, base=2)
    col_s = s[length(s) - 2:length(s)]
    column = parse(Int, col_s, base=2)
    seat = row * 8 + column
    result = BoardingPass(row, column, seat)
    # println("parse(", s, ") = ", result)
    return result
end

# examples
@assert toBoardingPass("FBFBBFFRLR") == BoardingPass(44, 5, 357)
@assert toBoardingPass("BFFFBBFRRR") == BoardingPass(70, 7, 567)
@assert toBoardingPass("FFFBBBFRRR") == BoardingPass(14, 7, 119)
@assert toBoardingPass("BBFFBBFRLL") == BoardingPass(102, 4, 820)

function run(fname)
    # read file
    f = open(fname, "r")
    lines = readlines(f)
    close(f)

    # compute BoardingPasses
    seatIDs = map((x) -> toBoardingPass(x).SeatID, lines)
    result = maximum(seatIDs) # Warning: max() != maximum()
    println("Ran ", fname, " got ", result)
    return result
end

function run2(fname)
    # read file
    f = open(fname, "r")
    lines = readlines(f)
    close(f)

    # compute BoardingPasses
    seatIDs = map((x) -> toBoardingPass(x).SeatID, lines)
    sortedSeatIDs = sort(seatIDs)
    result = -1
    for i in 1:length(sortedSeatIDs)
        if i == 1
            continue
        end
        if sortedSeatIDs[i] - sortedSeatIDs[i - 1] == 2
            result = sortedSeatIDs[i] - 1
            break
        end
    end
    println("Ran ", fname, " got ", result)
    return result
end

run("5.txt")
run2("5.txt")


text:
using Pkg
pkg"add https://github.com/asinghvi17/NASAEarthObservations.jl"
using GeoMakie, Makie, NASAEarthObservations

using ImageMagick: load

# download data - ~140 MB
imgdir = observations("geotiff/CERES_NETFLUX_M")

date_regex(dirname, ext) = Regex("$(dirname)_(\\d{4})-(\\d{2}).$(uppercase(ext))")
imflip(img) = reverse(vec(transpose(reverse(img; dims=2))))

img = load(joinpath(imgdir, "CERES_NETFLUX_M_2006-07.TIFF"))

lons = LinRange(-89, 90, size(img)[1])
lats = LinRange(-179, 180, size(img)[2])

source = LonLat()
dest = Projection("+proj=robin")

points = GeoMakie.transform.(LonLat(), dest, GeoMakie.gridpoints(lats, lons))
faces  = GeoMakie.grid_triangle_faces(lats, lons)

scene = poly(points, faces; color = imflip(img), show_axis = false);

geoaxis!(scene, -180, 180, -90, 90; crs = (dest = dest,));

titletext = Node("07/2016")

fullsc = title(sc, titletext; fontsize = 40);

save(DataFormat{:PNG}, "d.png", fullsc); lines(rand(10))

record(fullsc, "ceres_netflux.mp4", filter!(x -> uppercase(splitext(x)[2]) == ".TIFF", sort(readdir(imgdir))); framerate = 10) do img

    year, month = match(re, img).captures

    scene.plots[1].color = imflip(ImageMagick.load(joinpath(imgdir,img)))
    titletext[] = "$month/$year"
end

lines(rand(10))