text:
"""Julia program to implement Binary Search algorithm.
Binary Search is a Divide and Conquer strategy based algorithm, where the entire array is divided into two halfs and the
required element is searched for in one of these sections as per the conditions. It is done on sorted arrays.
"""

function binary_search(arr, n, ele)
    low = 1
    high = n
    while high >= low
        mid = low + (high - low) ÷ 2
        # If the mid element is the required element, return that index
        if(arr[mid] == ele)
            return mid
        # ele is greater than mid element then ele would be present in first half of the array
        elseif (arr[mid] > ele)
            high = mid - 1
        #Else if  ele is smaller than mid element then ele would be present in last half of the array
        else
            low = mid + 1
        end
    end
    # If the element is not found return 0
    return 0
end

print("How many numbers are present in the array? ")
n = readline()
n = parse(Int, n)
if (n <= 0)
    println("Array is Empty!!!")
    exit()
end
arr = Int[]
print("Enter the numbers: ")
arr = [parse(Int, num) for num in split(readline())]
print("Which number do you want to search in the array? ")
ele = readline()
ele = parse(Int, ele)
# Sort the array in ascending order
arr = sort(arr)
res = binary_search(arr, n, ele)
if (res == 0)
    print("The number $ele is not present in the array")
else
    print("The number $ele is present in the array.")
end

"""
Time Complexity - O(log(n)), where 'n' is the size of the array
Space Complexity - O(n)

SAMPLE INPUT AND OUTPUT

SAMPLE I

How many numbers are present in the array? 5
Enter the numbers: 1 2 3 4 5
Which number do you want to search in the array? 6
The number 6 is not present in the array

SAMPLE II

How many numbers are present in the array? 3
Enter the numbers: 3 1 2
Which number do you want to search in the array? 2
The number 2 is present in the array.

"""


text:
using OrdinaryDiffEq, Test, Sundials

f(u,p,t) = 2u
u0 = 0.5
tspan = (0.0,1.0)
prob = ODEProblem(f,u0,tspan)
sol = solve(prob,Tsit5())

function f(du, u, p, t)
    du .= 2.0 * u
end
prob = ODEProblem(f,u0,tspan)
@test_throws DiffEqBase.IncompatibleInitialConditionError sol = solve(prob,Tsit5())

prob = ODEProblem{false}(f,u0,tspan)
sol = solve(prob,Tsit5())

prob = ODEProblem{false}(f,1.0+im,tspan)
@test_throws DiffEqBase.ComplexSupportError solve(prob,CVODE_Adams())

@test_throws DiffEqBase.ProblemSolverPairingError solve(prob, DFBDF())
@test_throws DiffEqBase.NoDefaultAlgorithmError solve(prob,nothing)
@test_throws DiffEqBase.NonSolverError solve(prob,5.0)

prob = ODEProblem{false}(f,u0,(nothing,nothing))
@test_throws DiffEqBase.NoTspanError solve(prob,Tsit5())

prob = ODEProblem{false}(f,Any[1.0,1f0],tspan)
@test_throws DiffEqBase.NonConcreteEltypeError solve(prob,Tsit5())

prob = ODEProblem{false}(f,(1.0,1f0),tspan)
@test_throws DiffEqBase.TupleStateError solve(prob,Tsit5())

prob = ODEProblem{false}(f,u0,(0.0+im,1.0))
@test_throws DiffEqBase.ComplexTspanError solve(prob,Tsit5())

ext: jl


text:
import LinearAlgebra

## We use ispossemidef for numbers as well as matrices.
## This follows `ishermitian`, etc.

"""
    ispossemidef(x::Number, approx_test::AbstractApprox=Equal())

Return `true` if `x` is (approximately) non-negative.
"""
function ispossemidef(x::Real, approx_test::AbstractApprox=Equal())
    return x > zero(x) || iszero(x, approx_test)
end

## Used for both pos def and pos semidef
function _isposdef(z::Complex, approx_test::AbstractApprox, posdeffunc)
    return posdeffunc(real(z), approx_test) && iszero(imag(z), approx_test)
end

function ispossemidef(z::Complex, approx_test::AbstractApprox=Equal())
    return _isposdef(z, approx_test, ispossemidef)
end

function isposdef(z::Complex, approx_test::AbstractApprox=Equal())
    return _isposdef(z, approx_test, isposdef)
end

"""
    isposdef(x::Number, approx_test::AbstractApprox=Equal())

Return `true` if `x` is approximately greater than zero.
"""
isposdef(x::Number) = isposdef(x, Equal())
isposdef(x::Real, ::Equal) = x > zero(x)
## For methods other than `Equal`, x can be zero or negative
isposdef(x::Number, approx_test::AbstractApprox) = ispossemidef(x, approx_test)

## For both positive definite and positive semidefinite
function _isposdef(m::AbstractMatrix, approx_test::AbstractApprox, posdeffunc)
    ! ishermitian(m, approx_test) && return false
    evs = LinearAlgebra.eigvals(LinearAlgebra.Hermitian(m))
    return all(x -> posdeffunc(x, approx_test), evs)
end

"""
    ispossemidef(m::AbstractMatrix, approx_test::AbstractApprox=Equal())

Return `true` if `m` is positive semidefinite.
"""
function ispossemidef(m::AbstractMatrix, approx_test::AbstractApprox=Equal())
    return _isposdef(m, approx_test, ispossemidef)
end

"""
    isposdef(m::AbstractMatrix, approx_test::AbstractApprox=Equal())

Return `true` if `m` is positive definite.
"""
function isposdef(m::AbstractMatrix, approx_test::AbstractApprox)
    return _isposdef(m, approx_test, isposdef)
end

isposdef(A::AbstractMatrix) = isposdef(A, Equal())
# copied from dense.jl
isposdef(A::AbstractMatrix, ::Equal) =
    ishermitian(A, Equal()) && LinearAlgebra.isposdef(LinearAlgebra.cholesky(LinearAlgebra.Hermitian(A); check = false))

## Compared two methods:
## a) Allocate, ie m' * m. b) iterate over columns
## Found:
## 1. iterating order, ie rows vs cols is correct
## 2. Avoiding allocation improves efficiency even for 2x2 matrices, but...
## 3. Doing the allocation, eg m' * m can be slightly faster. Eg for 100x100 dense identity matrix.
## 4. For rand(100, 100), iterating over columns is 1000 times faster. Fails on first column.
## `approx_test` is `Equal` or `EachApprox`.
function _isunitary(m::AbstractMatrix, approx_test::AbstractApprox, dotf, transposef)
    rowinds = axes(m)[2]
    for i in rowinds
        isapprox(approx_test, dotf(view(m, :, i), view(transposef(m), :, i)), 1) || return false
        for j in i+1:last(rowinds)
            isapprox(approx_test, dotf(view(m, :, i), view(transposef(m), :, j)), 0) || return false
        end
    end
    return true
end

## This is slower even for small matrices.
## Then why am I using it ? (May 2021)
function isunitary(m::AbstractMatrix, approx_test::Approx)
    return  isapprox(approx_test, m' * m, LinearAlgebra.I)
end

_identity(x) = x
"""
    isunitary(m::AbstractMatrix, approx_test::AbstractApprox=Equal())

Return `true` if `m` is unitary. If `m` is real, this tests orthogonality.
"""
isunitary(m::AbstractMatrix, approx_test::AbstractApprox=Equal()) =
    _isunitary(m, approx_test, LinearAlgebra.dot, _identity)

# abs2 is much faster, but we would need to use sqrt to adjust the tolerance, thus losing any advantage.
isunitary(x::Number, approx_test::AbstractApprox=Equal()) = isone(abs(x), approx_test)
isunitary(J::LinearAlgebra.UniformScaling, approx_test::AbstractApprox=Equal()) = isunitary(J.λ, approx_test)

"""
    _dotu(x::AbstractVector, y::AbstractVector)

Dot product with no complex conjugation of `x`. This dispatches to `LinearAlgebra.dot`
if `x` and `y` are `Real`.
"""
_dot(x::AbstractVector{<:Real}, y::AbstractVector{<:Real}) = LinearAlgebra.dot(x, y)

_dot(x::AbstractVector{Complex{T}}, y::AbstractVector{Complex{T}}) where {T <: Union{Float64, Float32}} =
    LinearAlgebra.BLAS.dotu(x, y)

## There has got to be a better way than all of these _dot methods
function _dot(x::AbstractVector{<:Complex}, y::AbstractVector{<:Complex})
    return sum(*(z...) for z in zip(x, y))
end

"""
    isinvolution(m::AbstractMatrix, approx_test::AbstractApprox=Equal())

Return `true` if `m * m == I`
"""
isinvolution(m::AbstractMatrix, approx_test::AbstractApprox=Equal()) = _isunitary(m, approx_test, _dot, transpose)

function isinvolution(m::AbstractMatrix, approx_test::Approx)
    return  isapprox(approx_test, m * m, LinearAlgebra.I)
end

function isidempotent(m::AbstractMatrix, approx_test::AbstractApprox=Equal())
    return isapprox(approx_test, m * m, m)
end

function isnormal(m::AbstractMatrix, approx_test::AbstractApprox=Equal())
    return isapprox(approx_test, m * m', m' * m)
end

"""
    commutes(X, Y, approx_test::AbstractApprox=Equal())

Return `true` if `X` and `Y` commute.
"""
commutes(X, Y, approx_test::AbstractApprox=Equal()) = isapprox(approx_test, X * Y, Y * X)

"""
    anticommutes(X, Y, approx_test::AbstractApprox=Equal())

Return `true` if `X` and `Y` anticommute.
"""
anticommutes(X, Y, approx_test::AbstractApprox=Equal()) = isapprox(approx_test, X * Y, -(Y * X))

# Some of the following from QuantumInfo.jl
# might be implemented here:
# ischannel
# iscp
# isliouvillian
# istp
# isunital

# Yao.jl
#
# iscommute
# isnormalized
# isreflexive  <-- isinvolution is a better name, See above

# From other libraries
# isidempotent
# isnormal


text:
using ProgressMeter

function simulateChain!(chain::Vector{T}, P::F, x0::T) where {F<:Function, T}
  n::Int64 = length(chain)
  x::T = x0
  for i = 1:n
    x = P(x)
    @inbounds chain[i] = x
  end
end

function simulateChain(P::F, x0::T, n::Int64) where {F<:Function, T}
  chain::Vector{T} = Vector{T}(undef, n)
  simulateChain!(chain, P, x0)
  return chain
end

function simulateChainProgress!(chain::Vector{T}, P::F, x0::T) where
  {F<:Function, T}
  n::Int64 = length(chain)
  x::T = x0
  @showprogress 10 for i = 1:n
    x = P(x)
    @inbounds chain[i] = x
  end
  return chain
end

function simulateChainProgress(P::F, x0::T, n::Int64) where {F<:Function, T}
  chain::Vector{T} = Vector{T}(undef, n)
  simulateChainProgress!(chain, P, x0)
  return chain
end

function cov(xs::Vector{SVector{d, Float64}}) where d
  xbar = mean(xs)
  Q::SMatrix{d, d, Float64} = zeros(SMatrix{d, d, Float64})
  for i = 1:length(xs)
    @inbounds Q += (xs[i] - xbar)*(xs[i] - xbar)'
  end
  return Q / (length(xs)-1)
end


text:
import FinancialDSL
using Test, Dates
import InterestRates, BusinessDays

ts_start = Dates.now()

@testset "BS" begin
    include("test_BS.jl")
end

import FinancialDSL.Currencies.BRL
import FinancialDSL.Currencies.USD
import FinancialDSL.Currencies.EUR
import FinancialDSL.Currencies.GBP
import FinancialDSL.Currencies.JPY

@testset "Currencies" begin
    include("test_currencies.jl")
end

@testset "Core" begin
    include("test_core.jl")
end

@info("Tests completed in $(Dates.now() - ts_start)")


text:
# Text-mode spinners - from https://codepen.io/eLoyyyyy/full/yPdPPW

DOTS12 = [
			"⢀⠀",
			"⡀⠀",
			"⠄⠀",
			"⢂⠀",
			"⡂⠀",
			"⠅⠀",
			"⢃⠀",
			"⡃⠀",
			"⠍⠀",
			"⢋⠀",
			"⡋⠀",
			"⠍⠁",
			"⢋⠁",
			"⡋⠁",
			"⠍⠉",
			"⠋⠉",
			"⠋⠉",
			"⠉⠙",
			"⠉⠙",
			"⠉⠩",
			"⠈⢙",
			"⠈⡙",
			"⢈⠩",
			"⡀⢙",
			"⠄⡙",
			"⢂⠩",
			"⡂⢘",
			"⠅⡘",
			"⢃⠨",
			"⡃⢐",
			"⠍⡐",
			"⢋⠠",
			"⡋⢀",
			"⠍⡁",
			"⢋⠁",
			"⡋⠁",
			"⠍⠉",
			"⠋⠉",
			"⠋⠉",
			"⠉⠙",
			"⠉⠙",
			"⠉⠩",
			"⠈⢙",
			"⠈⡙",
			"⠈⠩",
			"⠀⢙",
			"⠀⡙",
			"⠀⠩",
			"⠀⢘",
			"⠀⡘",
			"⠀⠨",
			"⠀⢐",
			"⠀⡐",
			"⠀⠠",
			"⠀⢀",
			"⠀⡀"
		];

CORNERS = ['◴', '◷', '◶', '◵'];

function crank_spinner(spinner_pos::UInt16) :: UInt16
    spinners, backspace = CORNERS, "\b";
    # spinners, backspace = DOTS12, "\b\b";
    print(backspace);
    print(spinners[spinner_pos]);
    spinner_pos = (spinner_pos) % (length(spinners)) +1;
    return spinner_pos;
end

function test_spinner()
    spinner_pos = UInt16(1);
    while (1==1)
        spinner_pos = crank_the_spinner(spinner_pos);
        sleep(0.2);
    end
end
     # test_spinner()

max_stars_repo_name: thomasn/obsidify-keep

max_issues_repo_name: thomasn/obsidify-keep

max_forks_repo_name: thomasn/obsidify-keep


text:
@testset "Search strategies" begin
    include("dfs.jl")
    include("ilds.jl")
    include("rbs.jl")
end


text:
"""
    TensorBoardBackend(logdir[, tb_overwrite];
        time=time(),
        purge_step=nothing,
        step_increment=1,
        min_level=Logging.Info)

TensorBoard backend for logging callbacks. Takes the same arguments
as [`TensorBoardLogger.TBLogger`](https://julialogging.github.io/TensorBoardLogger.jl/dev/).
"""
struct TensorBoardBackend <: LoggerBackend
    logger::TBLogger
    function TensorBoardBackend(
            logdir,
            existfn = tb_overwrite; kwargs...)
        return new(TBLogger(logdir, existfn; kwargs...))
    end
end

Base.show(io::IO, backend::TensorBoardBackend) = print(
    io, "TensorBoardBackend(", backend.logger.logdir, ")")

function log_to(backend::TensorBoardBackend, value::Loggables.Value, name, i; group = ())
    name = _combinename(name, group)
    log_value(backend.logger, name, value.data; step = i)
end

function log_to(backend::TensorBoardBackend, image::Loggables.Image, name, i; group = ())
    name = _combinename(name, group)
    im = ImageCore.clamp01nan!(collect(image.data))
    log_image(backend.logger, name, im; step = i)
end

function log_to(backend::TensorBoardBackend, text::Loggables.Text, name, i; group = ())
    name = _combinename(name, group)
    log_text(backend.logger, name, text.data; step = i)
end

function log_to(backend::TensorBoardBackend, hist::Loggables.Histogram, name, i; group=())
    name = _combinename(name, group)
    log_histogram(backend.logger, name, hist.data, step=i)
end

# Utilities

_combinename(name, group::String) = _combinename((group, name))
_combinename(name, group::Tuple) = _combinename((group..., name))
_combinename(strings::Tuple) = join(strings, '/')


text:
export NVMLError

struct NVMLError <: Exception
    code::nvmlReturn_t
end

Base.convert(::Type{nvmlReturn_t}, err::NVMLError) = err.code

Base.showerror(io::IO, err::NVMLError) =
    print(io, "NVMLError: ", description(err), " (code $(reinterpret(Int32, err.code)))")

# FIXME: errors during init cannot use `nvmlErrorString`
description(err::NVMLError) = unsafe_string(nvmlErrorString(err))

@enum_without_prefix nvmlReturn_enum NVML_

## API call wrapper

# outlined functionality to avoid GC frame allocation
@noinline function throw_api_error(res)
    throw(NVMLError(res))
end

const initialized = Ref(false)
function initialize_api()
    if !initialized[]
        nvmlInitWithFlags(0)
        atexit() do
            nvmlShutdown()
        end
        initialized[] = true
    end
end

macro check(ex)
    quote
        res = $(esc(ex))
        if res != NVML_SUCCESS
            throw_api_error(res)
        end

        nothing
    end
end


text:
module Layouts

using ..Utilities, ..LayoutUtilities

export Layout

const default_colorlist = ["#1f77b4","#ff7f0e","#2ca02c","#d62728","#9467bd","#8c564b","#e377c2","#7f7f7f","#bcbd22","#17becf"]

struct Layout <: Component
    font::Font
    title::Title
    autosize::Bool
    width::RangeValue{Int,(10,Inf),(true,true)}
    height::RangeValue{Int,(10,Inf),(true,true)}
    margin::Margin
    paper_bgcolor::Color
    plot_bgcolor::Color
    seperators::String
    hidesources::Bool
    showlegend::Bool
    colorway::Vector{Color}
    colorscale::LayoutColorScale
    # datarevision - Not implemented
    # uirevision - Not implemented
    # editrevision - Not implemented
    # template - Not implemented
    modebar::ModeBar
    meta::Vector{Any}
    transition::Transition
    clickmode::ClickMode
    dragmode::DragMode
    hovermode::HoverMode
    hoverdistance::HoverDistance
    spikedistance::SpikeDistance
    hoverlabel::HoverLabel
    selectdirection::SelectDirection
    grid::Grid
    calendar::Calendar
    xaxis::Axis
    yaxis::Axis
    function Layout(;
        font::Font = Font(),
        title::Title = Title(),
        autosize::Bool = true,
        width::Int = 700,
        height::Int = 450,
        margin::Margin = Margin(),
        paper_bgcolor::String = "#fff",
        plot_bgcolor::String = "#fff",
        separators::String = ".,",
        hidesources::Bool = false,
        showlegend::Bool = true,
        colorway::Vector{String} = default_colorlist,
        colorscale::LayoutColorScale = LayoutColorScale(),
        modebar::ModeBar = ModeBar(),
        meta::Vector{Any} = [],
        transition::Transition = Transition(),
        clickmode::String = "event",
        dragmode::String = "zoom",
        hovermode::String = "closest",
        hoverdistance::Int = 20,
        spikedistance::Int = 20,
        hoverlabel::HoverLabel = HoverLabel(),
        selectdirection::String = "any",
        grid::Grid = Grid(),
        calendar::String = "gregorian",
        xaxis::Axis = Axis(),
        yaxis::Axis = Axis())
        return new(
            font,
            title,
            autosize,
            RangeValue((10,Inf),width),
            RangeValue((10,Inf),height),
            margin,
            Color(paper_bgcolor),
            Color(plot_bgcolor),
            separators,
            hidesources,
            showlegend,
            Color.(colorway),
            colorscale,
            modebar,
            meta,
            transition,
            ClickMode(clickmode),
            DragMode(dragmode),
            HoverMode(hovermode),
            HoverDistance(hoverdistance),
            SpikeDistance(spikedistance),
            hoverlabel,
            SelectDirection(selectdirection),
            grid,
            Calendar(calendar),
            xaxis,
            yaxis)
    end
end
default["layout"] = Layout()

end


text:
module RankingDifferences

using DifferencesBase

export
    AbstractRanking,
    Ranking

include("ranking.jl")

end # module


text:
using TOML
pkgs = TOML.parsefile("../Manifest.toml")

pkg_ver(pkg) = pkgs[pkg][1]["version"]

fp = open("output/index.md", "w")
header = """# Plots Gallery

This site is an unofficial Plots.jl Gallery site. It is inspired by [MATLAB Plot Gallery](https://jp.mathworks.com/products/matlab/plot-gallery.html).

Figures are ploted by [Plots.jl](https://github.com/JuliaPlots/Plots.jl).

"""

print(fp, header)
print(fp,
"""## Tested Environment
- Julia $(string(VERSION))
- Plots $(pkg_ver("Plots"))
- PyPlot $(pkg_ver("PyPlot"))
- GR $(pkg_ver("GR"))

""")

println(fp, """
For all plots, you should import Plots.
```julia
using Plots
```

""")

function list_plots()
    str = String(read("notebook.jl"))
    arr = []
    rx = r"include\(\"(.+)\.jl\"\)"
    offset = 1
    while (m = match(rx, str, offset)) != nothing
        filename = m.captures[1]
        push!(arr, filename)
        offset = m.offset + length(filename)
    end

    titles = copy(arr)
    for (i, name) in enumerate(arr)
        titles[i] = titlecase(replace(name, "_" => " "))
    end
    return arr, titles
end

function seperate(items, n)
    arr = []
    subarr = []
    cnt = 0
    for item in items
        cnt += 1
        push!(subarr, item)
        if cnt % n == 0
            push!(arr, subarr)
            subarr = []
        end
    end
    !isempty(subarr) && push!(arr, subarr)
    return arr
end

n = 4
filenames, titles = list_plots()
filenames_jagged = seperate(filenames, n)
titles_jagged = seperate(titles, n)

for (files_item, titles_item) in zip(filenames_jagged, titles_jagged)
    len = length(files_item)
    print(fp, "| ")
    for title in titles_item
        print(fp, " $(title) |")
    end
    println(fp)

    print(fp, "| ")
    for i in 1:len
        print(fp, " ---- |")
    end
    println(fp)

    print(fp, "| ")
    for name in files_item
        print(fp, "![small_$(name).png](src/images/small_$(name).png) |")
    end
    println(fp)

    print(fp, "| ")
    for name in files_item
        print(fp, "[View source](src/$(name).md) |")
    end
    println(fp, "\n\n<br>\n")
end

close(fp)


text:
module LinRegOutliers

# using StatsModels
# using DataFrames
# using Distributions
# using Clustering
# using StatsBase
# using LinearAlgebra
# using Plots
# using Optim

# import functions from corresponding packages
import StatsModels: @formula, FormulaTerm, ModelFrame, ModelMatrix
import DataFrames: DataFrame
import Distributions: TDist, Chisq, std, cov, median
import Clustering: Hclust, hclust, cutree
import StatsBase: quantile, standardize, ZScoreTransform, mean, sample
import LinearAlgebra: inv, pinv, det, diag, eigen, norm, rank, qr, Diagonal, UpperTriangular
import Plots: scatter, title!, xlabel!, ylabel!, hline!, vline!, RGB, plot, font, text
import Optim: optimize, NelderMead
import Combinatorics: combinations

# Basis
include("basis.jl")

# Predefined datasets used in outlier detection literature
include("data.jl")

# Regression diagnostics
include("diagnostics.jl")

# Ordinary least squares type and functions
# for fast regression tasks in outlier detection algorithms
include("ols.jl")

# Hadi & Simonoff (1993) algorithm
include("hs93.jl")

# Kianifard & Swallow (1989) algorithm
include("ks89.jl")

# Sebert et. al. (1998) algorithm
include("smr98.jl")

# asm (2000) algorithm
include("asm2000.jl")

# Rousseeuw's Least Median of Squares estimator
include("lms.jl")

# Rousseeuw's Least Trimmed Squares estimator
include("lts.jl")

# Minimum Volume Ellipsoid estimator
# for robust covariance matrix
include("mve.jl")

# MVE - LTS plot for visual detection of regression outliers
include("mveltsplot.jl")

# Billor & Chatterjee & Hadi Algorithm for detecting outliers
include("bch.jl")

# Pena & Yohai (1995) algorithm
include("py95.jl")

# Satman (2013) algorithm
include("satman2013.jl")

# Satman (2015) algorithm
include("satman2015.jl")

# Least Absolute Deviations estimator
include("lad.jl")

# Least Trimmed Absolute Deviations estimator
include("lta.jl")

# Hadi (1992) detecting outliers in multivariate data
# This algorithm is not directly related with the regression
# however, detecting bad leverage points is same as detecting
# outliers in X-space
include("hadi1992.jl")

# Hadi (1994) algorithm
include("hadi1994.jl")

# Gray-scale images of distance matrices
include("dataimage.jl")

# Compact genetic algorithm
include("cga.jl")

# Genetic Algorithm
include("ga.jl")

# Modified and original Satman (2012) algorithms
include("gwlts.jl")

# RANSAC Algorithm
include("ransac.jl")

# CCF formulation and heuristic; Barratt, Angeris, and Boyd (2020)
include("ccf.jl")

# Atkinson94 Algorithm
include("atkinson94.jl")
include("atkinsonstalactiteplot.jl")

# BACON 2000 Algorithm
include("bacon.jl")

# Imon 2005 Algorithm
include("imon2005.jl")

#  Chatterjee & Mächler (1997) Algorithm
include("cm97.jl")

# All-in-one
include("summary.jl")

# Essentials from other packages
export @formula, DataFrame
export mean, quantile

# Basics
export RegressionSetting
export createRegressionSetting
export designMatrix
export responseVector
export applyColumns
export find_minimum_nonzero
export @extractRegressionSetting

# Data
export phones, hbk, stackloss
export weightloss, hs93randomdata, woodgravity
export hills, softdrinkdelivery

# Diagnostics
export dffit, dfbeta
export hatmatrix
export studentizedResiduals
export adjustedResiduals
export jacknifedS
export cooks
export mahalanobisSquaredMatrix
export covratio
export hadimeasure

# Ordinary least squares
export OLS, ols, wls, residuals, predict, coef

# Algorithms
export hs93, hs93initialset, hs93basicsubset
export ks89
export smr98, asm2000
export lms
export lts
export mve, mcd
export mveltsplot
export bch, bchplot
export py95, py95SuspectedObservations
export satman2013
export satman2015, dominates
export lad
export lta
export hadi1992
export hadi1994
export dataimage
export gwcga, galts, ga, cga, RealChromosome
export detectOutliers
export ransac
export ccf
export imon2005
export atkinson94, atkinsonstalactiteplot, generate_stalactite_plot
export bacon
export cm97

end # module

max_stars_repo_name: angeris/LinRegOutliers

max_issues_repo_name: angeris/LinRegOutliers

max_forks_repo_name: angeris/LinRegOutliers


text:
@testset "graph" begin

let A = CUDA.zeros(Int, 1)
    # ensure compilation
    A .+= 1
    @test Array(A) == [1]

    graph = capture() do
        A .+= 1
    end
    @test Array(A) == [1]

    exec = instantiate(graph)
    CUDA.launch(exec)
    @test Array(A) == [2]

    graph′ = capture() do
        A .+= 2
    end

    update(exec, graph′)
    CUDA.launch(exec)
    @test Array(A) == [4]
end

let A = CUDA.zeros(Int, 1)
    function iteration(A, val)
        # custom kernel to force compilation on the first iteration
        function kernel(a, val)
            a[] += val
            return
        end
        @cuda kernel(A, val)
        return
    end

    for i in 1:2
        @captured iteration(A, i)
    end
    @test Array(A) == [3]
end

end


text:
"""
    abstract type GridTopology{Dc,Dp}

Abstract type representing the topological information associated with a grid.

The `GridTopology` interface is defined by overloading the methods:

- [`get_faces(g::GridTopology,dimfrom::Integer,dimto::Integer)`](@ref)
- [`get_polytopes(g::GridTopology)`](@ref)
- [`get_cell_type(g::GridTopology)`](@ref)
- [`get_vertex_coordinates(g::GridTopology)`](@ref)

The `GridTopology` interface has the following traits

- [`OrientationStyle(::Type{<:GridTopology})`](@ref)
- [`RegularityStyle(::Type{<:GridTopology})`](@ref)

and tested with this function:

- [`test_grid_topology`](@ref)

"""
abstract type GridTopology{Dc,Dp} <: GridapType end

# Traits

"""
    OrientationStyle(::Type{<:GridTopology})
    OrientationStyle(::GridTopology)

`Oriented()` if has oriented faces, `NonOriented()` otherwise (default).
"""
OrientationStyle(a::GridTopology) = OrientationStyle(typeof(a))

OrientationStyle(::Type{<:GridTopology}) = NonOriented()

"""
    RegularityStyle(::Type{<:GridTopology})
    RegularityStyle(::GridTopology)

`Regular()` if no hanging-nodes default), `Irregular()` otherwise.
"""
RegularityStyle(a::GridTopology) = RegularityStyle(typeof(a))

RegularityStyle(::Type{<:GridTopology}) = Regular()

# Abstract methods

"""
    get_faces(g::GridTopology,dimfrom::Integer,dimto::Integer)
"""
function get_faces(g::GridTopology,dimfrom::Integer,dimto::Integer)
  @abstractmethod
end

"""
    get_polytopes(g::GridTopology)
"""
function get_polytopes(g::GridTopology)
  @abstractmethod
end

"""
    get_cell_type(g::GridTopology)
"""
function get_cell_type(g::GridTopology)
  @abstractmethod
end

"""
    get_vertex_coordinates(g::GridTopology)
"""
function get_vertex_coordinates(g::GridTopology)
  @abstractmethod
end

"""
    test_grid_topology(top::GridTopology)
"""
function test_grid_topology(top::GridTopology{Dc,Dp}) where {Dc,Dp}
  D = num_cell_dims(top)
  @test D == Dc
  @test num_point_dims(top) == Dp
  @test num_dims(top) == D
  cell_to_type = get_cell_type(top)
  @test isa(cell_to_type,AbstractVector{<:Integer})
  @test length(cell_to_type) == num_cells(top)
  polytopes = get_polytopes(top)
  @test isa(polytopes,Vector{<:Polytope{D}})
  get_isboundary_face(top)
  get_cell_faces(top)
  get_face_vertices(top)
  @test OrientationStyle(top) in (Oriented(), NonOriented())
  @test RegularityStyle(top) in (Regular(), Irregular())
  @test is_oriented(top) == (OrientationStyle(top) == Oriented())
  @test is_regular(top) == (RegularityStyle(top) == Regular())
  for n in 0:D
    compute_reffaces(Polytope{n},top)
    for m in 0:D
      nface_to_mfaces = get_faces(top,n,m)
      @test isa(nface_to_mfaces,AbstractArray{<:Vector{<:Integer}})
      @test length(nface_to_mfaces) == num_faces(top,n)
    end
  end
end

# Default API

"""
    num_cell_dims(::GridTopology) -> Int
    num_cell_dims(::Type{<:GridTopology}) -> Int
"""
num_cell_dims(::GridTopology{Dc,Dp}) where {Dc,Dp} = Dc
num_cell_dims(::Type{<:GridTopology{Dc,Dp}}) where {Dc,Dp} = Dc

"""
    num_point_dims(::GridTopology) -> Int
    num_point_dims(::Type{<:GridTopology}) -> Int
"""
num_point_dims(::GridTopology{Dc,Dp}) where {Dc,Dp} = Dp
num_point_dims(::Type{<:GridTopology{Dc,Dp}}) where {Dc,Dp} = Dp

"""
    num_dims(::GridTopology) -> Int
    num_dims(::Type{<:GridTopology}) -> Int

Equivalent to `num_cell_dims`.
"""
num_dims(g::GridTopology{Dc}) where Dc = Dc
num_dims(::Type{<:GridTopology{Dc}}) where Dc = Dc

"""
    num_faces(g::GridTopology,d::Integer)
    num_faces(g::GridTopology)
"""
function num_faces(g::GridTopology,d::Integer)
  D = num_cell_dims(g)
  face_to_cells = get_faces(g,d,D)
  length(face_to_cells)
end

function num_faces(g::GridTopology)
  sum((num_faces(g,d) for d in 0:num_cell_dims(g)))
end

"""
    num_cells(g::GridTopology)
"""
num_cells(g::GridTopology) = num_faces(g,num_cell_dims(g))

"""
    num_facets(g::GridTopology)
"""
num_facets(g::GridTopology) = _num_facets(g)

"""
    num_edges(g::GridTopology)
"""
num_edges(g::GridTopology) = _num_edges(g)

"""
    num_vertices(g::GridTopology)
"""
num_vertices(g::GridTopology) = _num_vertices(g)

"""
    get_dimranges(g::GridTopology)
"""
function get_dimranges(g::GridTopology)
  ranges = UnitRange{Int}[]
  k = 1
  for d in 0:num_cell_dims(g)
    nf = num_faces(g,d)
    j = k+nf-1
    push!(ranges,k:j)
    k = j+1
  end
  ranges
end

"""
    get_dimrange(g::GridTopology,d::Integer)
"""
function get_dimrange(g::GridTopology,d::Integer)
  get_dimranges(g)[d+1]
end

"""
    get_offsets(g::GridTopology)
"""
get_offsets(g::GridTopology) = _get_offsets(g)

"""
    get_offset(g::GridTopology,d::Integer)
"""
get_offset(g::GridTopology,d::Integer) = _get_offset(g,d)

"""
    get_facedims(g::GridTopology)
"""
get_facedims(g::GridTopology) =  _get_facedims(Int8,g)

"""
    get_cell_faces(g::GridTopology)

Defaults to

    compute_cell_faces(g)
"""
function get_cell_faces(g::GridTopology)
  compute_cell_faces(g::GridTopology)
end

"""
    compute_cell_faces(g::GridTopology)
"""
function compute_cell_faces(g::GridTopology)
  D = num_cell_dims(g)
  d_to_cell_to_dface = Tuple((Table(get_faces(g,D,d)) for d in 0:D))
  offsets = Tuple(get_offsets(g))
  append_tables_locally(offsets,d_to_cell_to_dface)
end

"""
    get_face_vertices(g::GridTopology,d::Integer)
"""
function get_face_vertices(g::GridTopology,d::Integer)
  get_faces(g,d,0)
end

"""
    get_face_vertices(g::GridTopology)

Defaults to

    compute_face_vertices(g)
"""
function get_face_vertices(g::GridTopology)
  compute_face_vertices(g)
end

"""
    compute_face_vertices(g::GridTopology)
"""
function compute_face_vertices(g::GridTopology)
  D = num_cell_dims(g)
  d_to_dface_to_vertices = Tuple((Table(get_faces(g,d,0)) for d in 0:D))
  append_tables_globally(d_to_dface_to_vertices...)
end

"""
    get_cell_vertices(g::GridTopology)
"""
function get_cell_vertices(g::GridTopology)
  D = num_cell_dims(g)
  get_faces(g,D,0)
end

"""
    is_simplex(p::GridTopology) -> Bool
"""
function is_simplex(p::GridTopology)
  all(map(is_simplex, get_polytopes(p)))
end

"""
    is_n_cube(p::GridTopology) -> Bool
"""
function is_n_cube(p::GridTopology)
  all(map(is_n_cube, get_polytopes(p)))
end

"""
    is_oriented(::Type{<:GridTopology}) -> Bool
    is_oriented(a::GridTopology) -> Bool
"""
is_oriented(a::GridTopology) = is_oriented(typeof(a))
is_oriented(a::Type{T}) where T <:GridTopology = OrientationStyle(T) == Oriented()

"""
    is_regular(::Type{<:GridTopology}) -> Bool
    is_regular(a::GridTopology) -> Bool
"""
is_regular(a::GridTopology) = is_regular(typeof(a))
is_regular(a::Type{T}) where T<:GridTopology = RegularityStyle(T) == Regular()

"""
    get_reffaces(::Type{Polytope{d}}, g::GridTopology) where d

By default, it calls to `compute_reffaces`.
"""
function get_reffaces(::Type{Polytope{d}}, g::GridTopology) where d
  reffaces, _ = compute_reffaces(Polytope{d},g)
  reffaces
end

"""
    get_face_type(g::GridTopology,d::Integer)

By default, it calls to `compute_reffaces`.
"""
function get_face_type(g::GridTopology,d::Integer)
  _, face_to_ftype  = compute_reffaces(Polytope{d},g)
  face_to_ftype
end

"""
    compute_reffaces(::Type{Polytope{d}}, g::GridTopology) where d
"""
function compute_reffaces(::Type{Polytope{d}}, g::GridTopology) where d
  D = num_cell_dims(g)
  ctype_to_polytope = get_polytopes(g)
  ctype_to_lftype_to_refface = [ get_reffaces(Polytope{d},polytope) for polytope in ctype_to_polytope]
  ctype_to_lface_to_lftype = [ get_face_type(polytope,d) for polytope in ctype_to_polytope]
  t = _generate_ftype_to_refface(Val{d}(),ctype_to_lftype_to_refface,ctype_to_lface_to_lftype)
  ftype_to_refface, ctype_to_lface_to_ftype = t
  cell_to_faces = get_faces(g,D,d)
  cell_to_ctype = get_cell_type(g)
  nfaces = num_faces(g,d)
  face_to_ftype = generate_face_to_face_type(
    cell_to_faces, cell_to_ctype, ctype_to_lface_to_ftype, nfaces)

  (collect1d(ftype_to_refface), face_to_ftype)
end

function compute_reffaces(::Type{Polytope{D}}, g::GridTopology{D}) where D
  (get_polytopes(g), get_cell_type(g))
end

"""
    get_reffaces(topo::GridTopology)
"""
function get_reffaces(topo::GridTopology)
  reffaces, _ , _ = compute_reffaces(topo)
  reffaces
end

"""
    get_face_type(topo::GridTopology)
"""
function get_face_type(topo::GridTopology)
  _, face_to_ftype, _ = compute_reffaces(topo)
  face_to_ftype
end

"""
    get_reffaces_offsets(topo::GridTopology)
"""
function get_reffaces_offsets(topo::GridTopology)
  _, _, offsets = compute_reffaces(topo)
  offsets
end

"""
    compute_reffaces(g::GridTopology)
"""
function compute_reffaces(g::GridTopology)
  D = num_cell_dims(g)
  d_to_refdfaces = Vector{Polytope}[]
  d_to_dface_to_ftype = Vector{Int8}[]
  for d in 0:D
    push!(d_to_refdfaces,get_reffaces(Polytope{d},g))
    push!(d_to_dface_to_ftype,get_face_type(g,d))
  end
  d_to_offset = zeros(Int,D+1)
  for d in 1:D
    d_to_offset[d+1] = d_to_offset[d] + length(d_to_refdfaces[d])
    d_to_dface_to_ftype[d+1] .+= d_to_offset[d+1]
  end
  (collect(vcat(d_to_refdfaces...)), vcat(d_to_dface_to_ftype...), d_to_offset)
end

"""
    get_isboundary_face(g::GridTopology)
"""
function get_isboundary_face(g::GridTopology)
  compute_isboundary_face(g)
end

"""
    get_isboundary_face(g::GridTopology,d::Integer)
"""
function get_isboundary_face(g::GridTopology,d::Integer)
  compute_isboundary_face(g,d)
end

"""
    compute_isboundary_face(g::GridTopology)
"""
function compute_isboundary_face(g::GridTopology)
  D = num_cell_dims(g)
  d_to_dface_to_isboundary = [compute_isboundary_face(g,d) for d in 0:D]
  vcat(d_to_dface_to_isboundary...)
end

"""
    compute_isboundary_face(g::GridTopology,d::Integer)
"""
function compute_isboundary_face(g::GridTopology,d::Integer)
  D = num_cell_dims(g)
  if d == D-1
    _compute_isboundary_facet!(g)
  else
    _compute_isboundary_face!(g,d)
  end
end

function _compute_isboundary_facet!(g)
  D = num_cell_dims(g)
  d = D-1
  face_to_cells = get_faces(g,d,D)
  generate_facet_to_isboundary(face_to_cells)
end

function _compute_isboundary_face!(g,d)
  D = num_cell_dims(g)
  cell_to_faces = get_faces(g,D,d)
  cell_to_facets = get_faces(g,D,D-1)
  cell_to_ctype = get_cell_type(g)
  polytopes = get_polytopes(g)
  ctype_to_lface_to_lfacets = map( (p)->get_faces(p,d,D-1), polytopes)
  facet_to_isboundary = get_isboundary_face(g,D-1)
  nfaces = num_faces(g,d)

  generate_face_to_isboundary_from_cells(
    facet_to_isboundary,
    cell_to_faces,
    cell_to_facets,
    cell_to_ctype,
    ctype_to_lface_to_lfacets,
    nfaces)
end

"""
    get_cell_permutations(top::GridTopology)
"""
function get_cell_permutations(top::GridTopology)
  compute_cell_permutations(top)
end

"""
    get_cell_permutations(top::GridTopology,d::Integer)
"""
function get_cell_permutations(top::GridTopology,d::Integer)
  compute_cell_permutations(top,d)
end

"""
    compute_cell_permutations(top::GridTopology)
"""
function compute_cell_permutations(top::GridTopology)
  D = num_cell_dims(top)
  tables = (compute_cell_permutations(top,d) for d in 0:D)
  append_tables_locally(tables...)
end

"""
    compute_cell_permutations(top::GridTopology,d::Integer)
"""
function compute_cell_permutations(top::GridTopology,d::Integer)

  D = num_cell_dims(top)
  cell_to_lface_to_face = Table(get_faces(top,D,d))
  data = similar(cell_to_lface_to_face.data,Int8)
  ptrs = cell_to_lface_to_face.ptrs
  cell_to_lface_to_pindex = Table(data,ptrs)

  if d == D || d == 0
    fill!(cell_to_lface_to_pindex.data,Int8(1))
    return cell_to_lface_to_pindex
  end

  face_to_fvertex_to_vertex = Table(get_faces(top,d,0))
  face_to_ftype = get_face_type(top,d)
  reffaces = get_reffaces(Polytope{d},top)
  ftype_to_pindex_to_cfvertex_to_fvertex = map(get_vertex_permutations,reffaces)
  cell_to_cvertex_to_vertex = Table(get_faces(top,D,0))
  cell_to_ctype = get_cell_type(top)
  polytopes = get_polytopes(top)
  ctype_to_lface_to_cvertices = map( (p)->get_faces(p,d,0), polytopes )

  _compute_cell_perm_indices!(
    cell_to_lface_to_pindex,
    cell_to_lface_to_face,
    cell_to_cvertex_to_vertex,
    cell_to_ctype,
    ctype_to_lface_to_cvertices,
    face_to_fvertex_to_vertex,
    face_to_ftype,
    ftype_to_pindex_to_cfvertex_to_fvertex)

  cell_to_lface_to_pindex
end

"""
    GridTopology(grid::Grid)
    GridTopology(grid::Grid, cell_to_vertices::Table, vertex_to_node::Vector)
"""
function GridTopology(grid::Grid)
  _grid = UnstructuredGrid(grid)
  UnstructuredGridTopology(_grid)
end

function GridTopology(grid::Grid, cell_to_vertices::Table, vertex_to_node::AbstractVector)
  _grid = UnstructuredGrid(grid)
  UnstructuredGridTopology(_grid,cell_to_vertices,vertex_to_node)
end

# Helpers

function  _compute_cell_perm_indices!(
  cell_to_lface_to_pindex,
  cell_to_lface_to_face,
  cell_to_cvertex_to_vertex,
  cell_to_ctype,
  ctype_to_lface_to_cvertices,
  face_to_fvertex_to_vertex,
  face_to_ftype,
  ftype_to_pindex_to_cfvertex_to_fvertex)

  ncells = length(cell_to_lface_to_face)
  for cell in 1:ncells
    ctype = cell_to_ctype[cell]
    lface_to_cvertices = ctype_to_lface_to_cvertices[ctype]
    a = cell_to_lface_to_face.ptrs[cell]-1
    c = cell_to_cvertex_to_vertex.ptrs[cell]-1
    for (lface,cfvertex_to_cvertex) in enumerate(lface_to_cvertices)
      face = cell_to_lface_to_face.data[a+lface]
      ftype = face_to_ftype[face]
      b = face_to_fvertex_to_vertex.ptrs[face]-1
      pindex_to_cfvertex_to_fvertex = ftype_to_pindex_to_cfvertex_to_fvertex[ftype]
      pindexfound = false
      for (pindex, cfvertex_to_fvertex) in enumerate(pindex_to_cfvertex_to_fvertex)
        found = true
        for (cfvertex,fvertex) in enumerate(cfvertex_to_fvertex)
          vertex1 = face_to_fvertex_to_vertex.data[b+fvertex]
          cvertex = cfvertex_to_cvertex[cfvertex]
          vertex2 = cell_to_cvertex_to_vertex.data[c+cvertex]
          if vertex1 != vertex2
            found = false
            break
          end
        end
        if found
          cell_to_lface_to_pindex.data[a+lface] = pindex
          pindexfound = true
          break
        end
      end
      @assert pindexfound "Valid pindex not found"
    end
  end

end

function generate_cells_around(
  cell_to_faces::Table,
  nfaces::Integer = maximum(cell_to_faces.data))

  data, ptrs = _face_to_cells(cell_to_faces.data,cell_to_faces.ptrs,nfaces)
  Table(data,ptrs)
end

function generate_cell_to_faces(
  cell_to_vertices::Table,
  cell_type_to_lface_to_lvertices::Vector{Vector{Vector{Int}}},
  cell_to_cell_type::AbstractVector{<:Integer},
  vertex_to_cells::Table)

  data, ptrs = _generate_cell_to_faces(
    cell_to_vertices.data,
    cell_to_vertices.ptrs,
    cell_type_to_lface_to_lvertices,
    cell_to_cell_type,
    vertex_to_cells.data,
    vertex_to_cells.ptrs)

  Table(data,ptrs)
end

function _generate_ftype_to_refface(::Val{d},ctype_to_lftype_to_refface,ctype_to_lface_to_lftype) where d

  i_to_refface = vcat( ctype_to_lftype_to_refface... )

  i = 1
  ctype_to_lftype_to_i = Vector{Int}[]
  for ctype in 1:length(ctype_to_lftype_to_refface)
    lftype_to_i = Int[]
    for lftype in length(ctype_to_lftype_to_refface[ctype])
      push!(lftype_to_i, i)
      i +=1
    end
    push!(ctype_to_lftype_to_i, lftype_to_i)
  end

  ftype_to_refface, i_to_ftype = _find_unique_with_indices(i_to_refface)

  ctype_to_lftype_to_ftype = copy(ctype_to_lftype_to_i)
  for ctype in 1:length(ctype_to_lftype_to_i)
    for lftype in 1:length(ctype_to_lftype_to_i[ctype])
      i = ctype_to_lftype_to_i[ctype][lftype]
      ftype = i_to_ftype[i]
      ctype_to_lftype_to_ftype[ctype][lftype] = ftype
    end
  end

  ctype_to_lface_to_ftype = Vector{Int}[]
  for ctype in 1:length(ctype_to_lftype_to_ftype)
    lface_to_lftype = ctype_to_lface_to_lftype[ctype]
    lftype_to_ftype = ctype_to_lftype_to_ftype[ctype]
    lface_to_ftype = lftype_to_ftype[lface_to_lftype]
    push!(ctype_to_lface_to_ftype,lface_to_ftype)
  end

  (ftype_to_refface, ctype_to_lface_to_ftype)

end

function generate_face_to_face_type(
  cell_to_faces::Table,
  cell_to_cell_type::AbstractVector{<:Integer},
  cell_type_to_lface_to_face_type::Vector{Vector{T}},
  nfaces::Integer=maximum(cell_to_faces.data)) where T<:Integer

  _generate_face_to_ftype(
    cell_to_faces.data,
    cell_to_faces.ptrs,
    cell_to_cell_type,
    cell_type_to_lface_to_face_type,
    nfaces)

end

function generate_face_to_vertices(
  cell_to_vertices::Table,
  cell_to_faces::Table,
  cell_to_ctype::AbstractVector{<:Integer},
  ctype_to_lface_to_lvertices::Vector{Vector{Vector{Int}}},
  nfaces::Integer=maximum(cell_to_faces.data))

  data, ptrs = _generate_face_to_vertices(
    cell_to_vertices.data,
    cell_to_vertices.ptrs,
    cell_to_faces.data,
    cell_to_faces.ptrs,
    cell_to_ctype,
    ctype_to_lface_to_lvertices,
    nfaces)

  Table(data,ptrs)
end

function find_cell_to_faces(
  cell_to_vertices::Table,
  cell_type_to_lface_to_lvertices::Vector{Vector{Vector{Int}}},
  cell_to_cell_type::AbstractVector{<:Integer},
  vertex_to_faces::Table)

  data, ptrs = _cell_to_faces_from_vertex_to_faces(
    cell_to_vertices.data,
    cell_to_vertices.ptrs,
    cell_type_to_lface_to_lvertices,
    cell_to_cell_type,
    vertex_to_faces.data,
    vertex_to_faces.ptrs)

  Table(data,ptrs)
end

function generate_facet_to_isboundary(face_to_cells::Table)
  _generate_facet_to_isboundary(face_to_cells.ptrs)
end

function generate_face_to_isboundary(
  facet_to_isboundary::AbstractVector{Bool},
  face_to_facets::Table)

  _generate_face_to_isboundary(
    facet_to_isboundary,
    face_to_facets.data,
    face_to_facets.ptrs)
end

function generate_face_to_isboundary_from_cells(
  facet_to_isboundary,
  cell_to_faces::Table,
  cell_to_facets::Table,
  cell_to_ctype,
  ctype_to_lface_to_lfacets,
  nfaces = maximum(cell_to_faces.data))

  face_to_isboundary = fill(false,nfaces)

  _generate_face_to_isboundary_from_cells_fill!(
    face_to_isboundary,
    facet_to_isboundary,
    cell_to_faces.data,
    cell_to_faces.ptrs,
    cell_to_facets.data,
    cell_to_facets.ptrs,
    cell_to_ctype,
    ctype_to_lface_to_lfacets)

  face_to_isboundary

end

function  _generate_face_to_isboundary_from_cells_fill!(
    face_to_isboundary,
    facet_to_isboundary,
    cell_to_faces_data,
    cell_to_faces_ptrs,
    cell_to_facets_data,
    cell_to_facets_ptrs,
    cell_to_ctype,
    ctype_to_lface_to_lfacets)

  cells = 1:length(cell_to_faces_ptrs)-1

  for cell in cells
    a = cell_to_faces_ptrs[cell]-1
    b = cell_to_facets_ptrs[cell]-1
    ctype = cell_to_ctype[cell]
    lface_to_lfacets = ctype_to_lface_to_lfacets[ctype]
    for (lface, lfacets) in enumerate(lface_to_lfacets)
      face = cell_to_faces_data[a+lface]
      for lfacet in lfacets
        facet = cell_to_facets_data[b+lfacet]
        isboundary = facet_to_isboundary[facet]
        if isboundary
          face_to_isboundary[face] = isboundary
          continue
        end
      end
    end
  end

end

function _generate_tface_to_face(
  cell_to_faces_data::AbstractVector{T},
  cell_to_faces_ptrs,
  tcell_to_tfaces_data,
  tcell_to_tfaces_ptrs,
  ltcell_to_lnodes,
  ltface_to_ltnodes,
  lface_to_lnodes,
  ntfaces) where T

  ltcell_to_lfaces =  _generate_local_face_relations(
    ltcell_to_lnodes,
    ltface_to_ltnodes,
    lface_to_lnodes)

  _generate_tface_to_face(
    cell_to_faces_data,
    cell_to_faces_ptrs,
    tcell_to_tfaces_data,
    tcell_to_tfaces_ptrs,
    ltcell_to_lfaces,
    ntfaces)

end

function find_gface_to_face(
  face_to_nodes_data,
  face_to_nodes_ptrs,
  node_to_faces_data::AbstractVector{T},
  node_to_faces_ptrs,
  gface_to_nodes_data,
  gface_to_nodes_ptrs) where T

  ngfaces = length(gface_to_nodes_ptrs) - 1
  gface_to_face = zeros(T,ngfaces)
  n = max_cells_arround_vertex(node_to_faces_ptrs)
  faces_around = fill(UNSET,n)
  faces_around_scratch = fill(UNSET,n)

  _fill_gface_to_face!(
    gface_to_face,
    face_to_nodes_data,
    face_to_nodes_ptrs,
    node_to_faces_data,
    node_to_faces_ptrs,
    gface_to_nodes_data,
    gface_to_nodes_ptrs,
    faces_around,
    faces_around_scratch)

  gface_to_face

end

# Helpers

function _generate_face_to_isboundary(
  facet_to_isboundary::AbstractVector{Bool},
  face_to_facets_data::AbstractVector{<:Integer},
  face_to_facets_ptrs::AbstractVector{<:Integer})

  nobjects = length(face_to_facets_ptrs)-1
  object_to_isboundary = fill(false,nobjects)
  _generate_object_to_isboundary_fill!(
    object_to_isboundary,
    facet_to_isboundary,
    face_to_facets_data,
    face_to_facets_ptrs)
  object_to_isboundary
end

function _generate_facet_to_isboundary(face_to_cells_ptrs::AbstractVector{<:Integer})
  nfaces = length(face_to_cells_ptrs)-1
  face_to_isboundary = fill(false,nfaces)
  _generate_face_to_isboundary_fill!(face_to_isboundary,face_to_cells_ptrs)
  face_to_isboundary
end

function _generate_face_to_isboundary_fill!(
  face_to_isboundary, face_to_cells_ptrs)
  nfaces = length(face_to_isboundary)
  for face in 1:nfaces
    ncells_around = face_to_cells_ptrs[face+1] - face_to_cells_ptrs[face]
    if ncells_around == 1
      face_to_isboundary[face] = true
    else
      face_to_isboundary[face] = false
    end
  end
end

function  _generate_object_to_isboundary_fill!(
  object_to_isboundary,
  face_to_isboundary,
  object_to_faces_data,
  object_to_faces_ptrs)

  nobjects = length(object_to_isboundary)
  for object in 1:nobjects
    a = object_to_faces_ptrs[object]-1
    b = object_to_faces_ptrs[object+1]
    nlfaces = b-(a+1)
    for lface in 1:nlfaces
      face = object_to_faces_data[a+lface]
      isboundary = face_to_isboundary[face]
      if isboundary
        object_to_isboundary[object] = true
        break
      end
    end
  end

end

function _face_to_cells(
  cell_to_faces_data::AbstractVector{L},
  cell_to_faces_ptrs::AbstractVector{P},
  nfaces) where {L,P}

  face_to_cells_ptrs = zeros(P,nfaces+1)

  _face_to_cells_count!(
    face_to_cells_ptrs, cell_to_faces_data)

  length_to_ptrs!(face_to_cells_ptrs)

  ndata = face_to_cells_ptrs[end]-1

  face_to_cells_data = Vector{L}(undef,ndata)

  _face_to_cells_fill!(
    face_to_cells_data, face_to_cells_ptrs,
    cell_to_faces_data, cell_to_faces_ptrs)

  rewind_ptrs!(face_to_cells_ptrs)

  (face_to_cells_data, face_to_cells_ptrs)

end

function _face_to_cells_count!(
    face_to_cells_ptrs, cell_to_faces_data)
  for i=1:length(cell_to_faces_data)
    face_to_cells_ptrs[cell_to_faces_data[i]+1]+=1
  end
end

function _face_to_cells_fill!(
    face_to_cells_data, face_to_cells_ptrs,
    cell_to_faces_data, cell_to_faces_ptrs)

  ncells = length(cell_to_faces_ptrs) - 1
  for cell in 1:ncells
    a = cell_to_faces_ptrs[cell]
    b = cell_to_faces_ptrs[cell+1]-1
    for p in a:b
      face = cell_to_faces_data[p]
      q = face_to_cells_ptrs[face]
      face_to_cells_data[q] = cell
      face_to_cells_ptrs[face] += 1
    end
  end

end

function _generate_cell_to_faces(
    cell_to_vertices_data::AbstractVector{L},
    cell_to_vertices_ptrs::AbstractVector{P},
    ctype_to_lface_to_lvertices,
    cell_to_ctype,
    vertex_to_cells_data,
    vertex_to_cells_ptrs) where {L,P}

  cell_to_faces_ptrs = _cell_to_faces_count(
    cell_to_vertices_ptrs,
    cell_to_ctype,
    ctype_to_lface_to_lvertices)

  length_to_ptrs!(cell_to_faces_ptrs)

  ndata = cell_to_faces_ptrs[end]-1

  cell_to_faces_data = fill(L(UNSET),ndata)

  nvertices = max_nvertices_in_lface(ctype_to_lface_to_lvertices)
  vertices = fill(UNSET,nvertices)
  vertices_scratch = fill(UNSET,nvertices)

  nvertices = max_cells_arround_vertex(vertex_to_cells_ptrs)
  cells_around = fill(UNSET,nvertices)
  cells_around_scratch = fill(UNSET,nvertices)
  lface_of_cells_around = fill(UNSET,nvertices)

  _cell_to_faces_fill!(
    cell_to_faces_data,
    cell_to_faces_ptrs,
    cell_to_vertices_data,
    cell_to_vertices_ptrs,
    ctype_to_lface_to_lvertices,
    cell_to_ctype,
    vertex_to_cells_data,
    vertex_to_cells_ptrs,
    vertices,
    vertices_scratch,
    cells_around,
    cells_around_scratch,
    lface_of_cells_around)

  (cell_to_faces_data, cell_to_faces_ptrs)

end

function  _cell_to_faces_from_vertex_to_faces(
  cell_to_vertices_data::AbstractVector{L},
  cell_to_vertices_ptrs::AbstractVector{P},
  ctype_to_lface_to_lvertices,
  cell_to_ctype,
  vertex_to_faces_data,
  vertex_to_faces_ptrs) where {L,P}

  cell_to_faces_ptrs = _cell_to_faces_count(
    cell_to_vertices_ptrs,
    cell_to_ctype,
    ctype_to_lface_to_lvertices)

  length_to_ptrs!(cell_to_faces_ptrs)

  ndata = cell_to_faces_ptrs[end]-1

  cell_to_faces_data = fill(L(UNSET),ndata)

  nvertices = max_nvertices_in_lface(ctype_to_lface_to_lvertices)
  vertices = fill(UNSET,nvertices)

  nvertices = max_cells_arround_vertex(vertex_to_faces_ptrs)
  faces_around = fill(UNSET,nvertices)
  faces_around_scratch = fill(UNSET,nvertices)

  _cell_to_faces_fill!(
    cell_to_faces_ptrs,
    cell_to_faces_data,
    cell_to_vertices_data,
    cell_to_vertices_ptrs,
    ctype_to_lface_to_lvertices,
    cell_to_ctype,
    vertex_to_faces_data,
    vertex_to_faces_ptrs,
    vertices,
    faces_around,
    faces_around_scratch)

  (cell_to_faces_data, cell_to_faces_ptrs)
end

function _cell_to_faces_count(
  cell_to_vertices_ptrs::AbstractVector{P},
  cell_to_ctype,
  ctype_to_lface_to_lvertices) where P

  ncells = length(cell_to_vertices_ptrs) - 1

  cell_to_faces_ptrs = zeros(P,ncells+1)

  type_to_nlfaces = _type_to_nlfaces(ctype_to_lface_to_lvertices)

  _cell_to_faces_count!(
    cell_to_faces_ptrs,
    type_to_nlfaces,
    cell_to_ctype)

  cell_to_faces_ptrs

end

function _cell_to_faces_fill!(
  cell_to_faces_ptrs,
  cell_to_faces_data,
  cell_to_vertices_data,
  cell_to_vertices_ptrs,
  ctype_to_lface_to_lvertices,
  cell_to_ctype,
  vertex_to_faces_data,
  vertex_to_faces_ptrs,
  vertices,
  faces_around,
  faces_around_scratch)

  ncells = length(cell_to_ctype)

  for cell in 1:ncells

    ctype = cell_to_ctype[cell]
    lface_to_lvertices = ctype_to_lface_to_lvertices[ctype]
    nlfaces = length(lface_to_lvertices)
    a = cell_to_faces_ptrs[cell]-1

    for lface in 1:nlfaces

      _fill_vertices_in_lface!(
        vertices,
        lface,
        lface_to_lvertices,
        cell,
        cell_to_vertices_data,
        cell_to_vertices_ptrs)

      _find_cells_around_vertices!(
        faces_around,
        faces_around_scratch,
        vertices,
        vertex_to_faces_data,
        vertex_to_faces_ptrs)

      for face in faces_around
        if face != UNSET
          cell_to_faces_data[a+lface] = face
          break
        end
      end

    end

  end

end

function max_nvertices_in_lface(ctype_to_lface_to_lvertices)
  n = 0
  for lface_to_lvertices in ctype_to_lface_to_lvertices
    for lvertices in lface_to_lvertices
      m = length(lvertices)
      n = max(n,m)
    end
  end
  n
end

function max_cells_arround_vertex(vertex_to_cells_ptrs)
  n = 0
  nvertices = length(vertex_to_cells_ptrs)-1
  for vertex in 1:nvertices
    m = vertex_to_cells_ptrs[vertex+1]-vertex_to_cells_ptrs[vertex]
    n = max(n,m)
  end
  n
end

function _type_to_nlfaces(ctype_to_lface_to_lvertices)
  [ length(lface_to_lvertices)
   for lface_to_lvertices in ctype_to_lface_to_lvertices]
end

function _cell_to_faces_count!(
    cell_to_faces_ptrs,
    type_to_nlfaces,
    cell_to_ctype)

  ncells = length(cell_to_ctype)
  for cell in 1:ncells
    ctype = cell_to_ctype[cell]
    nlfaces = type_to_nlfaces[ctype]
    cell_to_faces_ptrs[cell+1] = nlfaces
  end

end

function  _cell_to_faces_fill!(
    cell_to_faces_data,
    cell_to_faces_ptrs,
    cell_to_vertices_data,
    cell_to_vertices_ptrs,
    ctype_to_lface_to_lvertices,
    cell_to_ctype,
    vertex_to_cells_data,
    vertex_to_cells_ptrs,
    vertices,
    vertices_scratch,
    cells_around,
    cells_around_scratch,
    lface_of_cells_around)

  face = 1

  ncells = length(cell_to_ctype)

  for cell in 1:ncells

    ctype = cell_to_ctype[cell]
    lface_to_lvertices = ctype_to_lface_to_lvertices[ctype]
    nlfaces = length(lface_to_lvertices)
    a = cell_to_faces_ptrs[cell]-1

    for lface in 1:nlfaces

      if cell_to_faces_data[a+lface] != UNSET
        continue
      end

      _fill_vertices_in_lface!(
        vertices,
        lface,
        lface_to_lvertices,
        cell,
        cell_to_vertices_data,
        cell_to_vertices_ptrs)

      _find_cells_around_vertices!(
        cells_around,
        cells_around_scratch,
        vertices,
        vertex_to_cells_data,
        vertex_to_cells_ptrs)

      _find_lface_of_cells_around_lface!(
        lface_of_cells_around,
        cells_around,
        lface,
        ctype_to_lface_to_lvertices,
        cell,
        cell_to_vertices_data,
        cell_to_vertices_ptrs,
        cell_to_ctype,
        vertices,
        vertices_scratch)

      _fill_face_in_cells_arround!(
        cell_to_faces_data,
        cell_to_faces_ptrs,
        face,
        cells_around,
        lface_of_cells_around)

      face += 1
    end

  end

end

function _fill_face_in_cells_arround!(
  cell_to_faces_data,
  cell_to_faces_ptrs,
  face,
  cells_around,
  lface_of_cells_around)

  for icell_around in 1:length(cells_around)
    cell_around = cells_around[icell_around]
    if cell_around == UNSET
      continue
    end
    lface = lface_of_cells_around[icell_around]
    f = cell_to_faces_ptrs[cell_around]-1
    cell_to_faces_data[f+lface] = face
  end

end

function _find_cells_around_vertices!(
  cells_around,
  cells_around_scratch,
  vertices,
  vertex_to_cells_data,
  vertex_to_cells_ptrs)

  ncells_around = UNSET
  ncells_around_scratch = UNSET

  for ivertex in 1:length(vertices)
    vertex = vertices[ivertex]
    if vertex == UNSET
      continue
    end
    if ivertex == 1
      ncells_around = _fill_cells_around_scratch!(
        cells_around,
        vertex,
        vertex_to_cells_data,
        vertex_to_cells_ptrs)
    else
      ncells_around_scratch = _fill_cells_around_scratch!(
        cells_around_scratch,
        vertex,
        vertex_to_cells_data,
        vertex_to_cells_ptrs)
      _set_intersection!(
        cells_around,cells_around_scratch,
        ncells_around,ncells_around_scratch)
    end
  end

end

function _fill_cells_around_scratch!(
  cells_around_scratch,
  vertex,
  vertex_to_cells_data,
  vertex_to_cells_ptrs)

  cells_around_scratch .= UNSET
  d = vertex_to_cells_ptrs[vertex] - 1
  ncells_around = vertex_to_cells_ptrs[vertex+1] - (d + 1)
  for icell_around in 1:ncells_around
    cell_around = vertex_to_cells_data[d+icell_around]
    cells_around_scratch[icell_around] = cell_around
  end
  ncells_around

end

function _set_intersection!(
  cells_around,cells_around_scratch,
  ncells_around,ncells_around_scratch)
  for i in 1:ncells_around
    if cells_around[i] == UNSET
      continue
    end
    _find_eq!(i,cells_around,cells_around_scratch,ncells_around_scratch)
  end
end

function _find_eq!(i,cells_around,cells_around_scratch,ncells_around_scratch)
  for j in 1:ncells_around_scratch
    if cells_around[i] == cells_around_scratch[j]
      return
    end
  end
  cells_around[i] = UNSET
  return
end

function _find_lface_of_cells_around_lface!(
  lface_of_cells_around,
  cells_around,
  lface,
  ctype_to_lface_to_lvertices,
  cell,
  cell_to_vertices_data,
  cell_to_vertices_ptrs,
  cell_to_ctype,
  vertices,
  vertices_scratch)

  lface_of_cells_around .= UNSET

  for icell_around in 1:length(cells_around)
    cell_around = cells_around[icell_around]
    if cell_around == UNSET
      continue
    end

    lface_of_cell_around = _find_lface_with_same_vertices(
      vertices,
      vertices_scratch,
      cell_around,
      cell_to_vertices_data,
      cell_to_vertices_ptrs,
      cell_to_ctype,
      ctype_to_lface_to_lvertices)

    lface_of_cells_around[icell_around] = lface_of_cell_around

  end

end

function _find_lface_with_same_vertices(
  vertices,
  vertices_scratch,
  cell,
  cell_to_vertices_data,
  cell_to_vertices_ptrs,
  cell_to_ctype,
  ctype_to_lface_to_lvertices)

  ctype = cell_to_ctype[cell]
  lface_to_lvertices = ctype_to_lface_to_lvertices[ctype]

  nlfaces = length(lface_to_lvertices)
  c = cell_to_vertices_ptrs[cell]-1

  for lface in 1:nlfaces

    _fill_vertices_in_lface!(
      vertices_scratch,
      lface,
      lface_to_lvertices,
      cell,
      cell_to_vertices_data,
      cell_to_vertices_ptrs)

    if _set_equal(vertices,vertices_scratch)
      return lface
    end

  end

  return UNSET

end

function _fill_vertices_in_lface!(
  vertices,
  lface,
  lface_to_lvertices,
  cell,
  cell_to_vertices_data,
  cell_to_vertices_ptrs)

  vertices .= UNSET
  c = cell_to_vertices_ptrs[cell] - 1
  lvertices = lface_to_lvertices[lface]
  nlfvertex = length(lvertices)
  for lfvertex in 1:nlfvertex
    lvertex = lvertices[lfvertex]
    vertex = cell_to_vertices_data[c+lvertex]
    vertices[lfvertex] = vertex
  end

end

function _set_equal(vertices,vertices_scratch)

  b = _is_subset(vertices,vertices_scratch)
  if b == false; return false; end
  b = _is_subset(vertices_scratch,vertices)
  if b == false; return false; end
  return true

end

function _is_subset(vertices,vertices_scratch)
  for i in 1:length(vertices)
    v = vertices[i]
    if v == UNSET
      continue
    end
    b = _find_eq(v,vertices_scratch)
    if b == false; return false; end
  end
  return true
end

function _find_eq(v,vertices_scratch)
  for vs in vertices_scratch
    if v == vs
      return true
    end
  end
  return false
end

function _generate_face_to_ftype(
  cell_to_faces_data,
  cell_to_faces_ptrs,
  cell_to_ctype,
  ctype_to_lface_to_ftype::Vector{Vector{T}},
  nfaces) where T

  face_to_ftype = fill(Int8(UNSET),nfaces)
  _face_to_ftype_fill!(
    face_to_ftype,
    cell_to_faces_data,
    cell_to_faces_ptrs,
    cell_to_ctype,
    ctype_to_lface_to_ftype)
  face_to_ftype
end

function _face_to_ftype_fill!(
  face_to_ftype,
  cell_to_faces_data,
  cell_to_faces_ptrs,
  cell_to_ctype,
  ctype_to_lface_to_ftype)

  ncells = length(cell_to_ctype)
  for cell in 1:ncells
    ctype = cell_to_ctype[cell]
    lface_to_ftype = ctype_to_lface_to_ftype[ctype]
    a = cell_to_faces_ptrs[cell]-1
    nlfaces = cell_to_faces_ptrs[cell+1] - (a + 1)
    for lface in 1:nlfaces
      face = cell_to_faces_data[a+lface]
      if face_to_ftype[face] != UNSET
        continue
      end
      ftype = lface_to_ftype[lface]
      face_to_ftype[face] = ftype
    end
  end

end

function _generate_face_to_vertices(
  cell_to_vertices_data::AbstractVector{L},
  cell_to_vertices_ptrs::AbstractVector{P},
  cell_to_faces_data,
  cell_to_faces_ptrs,
  cell_to_ctype,
  ctype_to_lface_to_lvertices,
  nfaces) where {L,P}

  face_to_vertices_ptrs = fill(P(UNSET),nfaces+1)

  _face_to_vertices_count!(
    face_to_vertices_ptrs,
    cell_to_faces_data,
    cell_to_faces_ptrs,
    cell_to_ctype,
    ctype_to_lface_to_lvertices)

  length_to_ptrs!(face_to_vertices_ptrs)
  ndata = face_to_vertices_ptrs[end]-1
  face_to_vertices_data = fill(L(UNSET),ndata)

  _face_to_vertices_fill!(
    face_to_vertices_data,
    face_to_vertices_ptrs,
    cell_to_vertices_data,
    cell_to_vertices_ptrs,
    cell_to_faces_data,
    cell_to_faces_ptrs,
    cell_to_ctype,
    ctype_to_lface_to_lvertices)

  (face_to_vertices_data, face_to_vertices_ptrs)

end

function _face_to_vertices_count!(
  face_to_vertices_ptrs,
  cell_to_faces_data,
  cell_to_faces_ptrs,
  cell_to_ctype,
  ctype_to_lface_to_lvertices)

  ncells = length(cell_to_ctype)
  for cell in 1:ncells
    ctype = cell_to_ctype[cell]
    lface_to_lvertices = ctype_to_lface_to_lvertices[ctype]
    a = cell_to_faces_ptrs[cell]-1
    nlfaces = cell_to_faces_ptrs[cell+1] - (a + 1)
    for lface in 1:nlfaces
      face = cell_to_faces_data[a+lface]
      if face_to_vertices_ptrs[face+1] != UNSET
        continue
      end
      nfvertices = length(lface_to_lvertices[lface])
      face_to_vertices_ptrs[face+1] = nfvertices
    end
  end

end

function _face_to_vertices_fill!(
  face_to_vertices_data,
  face_to_vertices_ptrs,
  cell_to_vertices_data,
  cell_to_vertices_ptrs,
  cell_to_faces_data,
  cell_to_faces_ptrs,
  cell_to_ctype,
  ctype_to_lface_to_lvertices)

  ncells = length(cell_to_ctype)
  for cell in 1:ncells
    ctype = cell_to_ctype[cell]
    lface_to_lvertices = ctype_to_lface_to_lvertices[ctype]
    a = cell_to_faces_ptrs[cell]-1
    c = cell_to_vertices_ptrs[cell]-1
    nlfaces = cell_to_faces_ptrs[cell+1] - (a + 1)
    for lface in 1:nlfaces
      face = cell_to_faces_data[a+lface]
      v = face_to_vertices_ptrs[face]-1
      lvertices = lface_to_lvertices[lface]
      nfvertices = length(lvertices)
      if nfvertices==0 || face_to_vertices_data[v+1] != UNSET
        continue
      end
      for lfvertex in 1:nfvertices
        lvertex = lvertices[lfvertex]
        vertex = cell_to_vertices_data[c+lvertex]
        face_to_vertices_data[v+lfvertex] = vertex
      end
    end
  end

end

function _generate_local_face_relations(
  ltcell_to_lnodes, ltface_to_ltnodes, lface_to_lnodes)

  nltcells = length(ltcell_to_lnodes)
  nltfaces = length(ltface_to_ltnodes)
  nlfaces = length(lface_to_lnodes)

  ltcell_ltface_to_lface = [ zeros(Int,nltfaces) for ltcell in 1:nltcells  ]

  for ltcell in 1:nltcells

    ltnode_to_lnode = ltcell_to_lnodes[ltcell]

    for ltface in 1:nltfaces
      ltnodes = ltface_to_ltnodes[ltface]
      for lface in 1:nlfaces
        lnodes = lface_to_lnodes[lface]
        allin = true
        for ltnode in ltnodes
          lnode = ltnode_to_lnode[ltnode]
          if !(lnode in lnodes)
            allin = false
            break
          end
        end
        if allin
          ltcell_ltface_to_lface[ltcell][ltface] = lface
          break
        end
      end
    end

  end

  ltcell_ltface_to_lface

end

function _generate_tface_to_face(
  cell_to_faces_data::AbstractVector{T},
  cell_to_faces_ptrs,
  tcell_to_tfaces_data,
  tcell_to_tfaces_ptrs,
  ltcell_to_lfaces,
  ntfaces) where T

  tface_to_face = zeros(T,ntfaces)

  _generate_tface_to_face!(
    tface_to_face,
    cell_to_faces_data,
    cell_to_faces_ptrs,
    tcell_to_tfaces_data,
    tcell_to_tfaces_ptrs,
    ltcell_to_lfaces)

  tface_to_face

end

function  _generate_tface_to_face!(
  tface_to_face,
  cell_to_faces_data,
  cell_to_faces_ptrs,
  tcell_to_tfaces_data,
  tcell_to_tfaces_ptrs,
  ltcell_to_lfaces)

  ncells = length(cell_to_faces_ptrs)-1
  nltcells = length(ltcell_to_lfaces)

  tcell = 1
  for cell in 1:ncells
    c = cell_to_faces_ptrs[cell]-1
    for ltcell in 1:nltcells
      a = tcell_to_tfaces_ptrs[tcell]-1
      b = tcell_to_tfaces_ptrs[tcell+1]
      nltfaces = b - (a+1)
      ltface_to_lface = ltcell_to_lfaces[ltcell]
      for ltface in 1:nltfaces
        tface = tcell_to_tfaces_data[a+ltface]
        lface = ltface_to_lface[ltface]
        if lface != 0
          face = cell_to_faces_data[c+lface]
          tface_to_face[tface] = face
        end
      end
      tcell += 1
    end
  end

end

function  _fill_gface_to_face!(
  gface_to_face,
  face_to_nodes_data,
  face_to_nodes_ptrs,
  node_to_faces_data,
  node_to_faces_ptrs,
  gface_to_nodes_data,
  gface_to_nodes_ptrs,
  faces_around,
  faces_around_scratch)

  ngfaces = length(gface_to_nodes_ptrs) - 1

  nfaces_around = UNSET
  nfaces_around_scratch = UNSET

  for gface in 1:ngfaces

    a = gface_to_nodes_ptrs[gface]-1
    b = gface_to_nodes_ptrs[gface+1]
    nlnodes = b-(a+1)

    for lnode in 1:nlnodes
      node = gface_to_nodes_data[lnode+a]
      if lnode == 1
        nfaces_around = _fill_cells_around_scratch!(
          faces_around,
          node,
          node_to_faces_data,
          node_to_faces_ptrs)
      else
        nfaces_around_scratch = _fill_cells_around_scratch!(
          faces_around_scratch,
          node,
          node_to_faces_data,
          node_to_faces_ptrs)
        _set_intersection!(
          faces_around,faces_around_scratch,
          nfaces_around,nfaces_around_scratch)
      end
    end

    for face in faces_around
      if face != UNSET
        gface_to_face[gface] = face
        break
      end
    end

  end

end


text:
module MyTemplate

# TODO: [Develop](https://invenia.github.io/PkgTemplates.jl/dev/user/#PkgTemplates.Develop)
# TODO: [Consider giving test its own project](https://invenia.github.io/PkgTemplates.jl/dev/user/#PkgTemplates.Tests)
# TODO: [Benchmark](https://github.com/invenia/PkgTemplates.jl/issues/86)
# TODO: [Hooks](https://github.com/invenia/PkgTemplates.jl/issues/11)

using PkgTemplates: Codecov, Documenter, Git, GitHubActions, License, Template, Tests, TravisCI

import ..FlexibleLicense
import ..TestsInSrc

"""
    template(kwargs...)

PkgTemplate templates the way I configure 'em.

Assumes use of GitHub.
Turns onn GitHub services if `services`.
Turns on free-for-open-source services if `services` and `!is_proprietary`.
Marks all rights reserved by owner if `is_proprietary`.
Sets up testing to look for *_test.jl files within src directory.
"""
function template(;
    parent_dir::String,
    is_pkg::Bool,
    is_proprietary::Bool,
    owner::Union{Nothing,String}=nothing,
    min_julia_version::VersionNumber=v"1.6",
    services::Bool=true,
    platforms=(; linux=true, osx=false, windows=false, x64=true, x86=false, arm64=false),
)
    test_arm = platforms.arm64
    @assert !is_proprietary || !test_arm

    # TODO: "UNLICENSED" or "Proprietary"?
    license_name = is_proprietary ? "UNLICENSED" : "MIT"

    plugins = Any[
        Git(; manifest=!is_pkg, ssh=true),
        !License,
        FlexibleLicense(; owner, name=license_name),
        Tests(; file=TestsInSrc.TEMPLATE_PATH, project=false),
    ]

    if services
        gh_platforms = (; filter(kv -> first(kv) != :arm64, pairs(platforms))...)

        @info """The GitHub Actions 'Documentation' job will not trigger gh-pages deploys until you
        (having admin rights) have manually pushed to the created 'gh-pages' branch at least once."""

        push!(plugins, Documenter{GitHubActions}(), GitHubActions(; gh_platforms...))

        if test_arm
            # n.b. In addition to TravisCI, DroneCI also supports arm64 (and more) and is also free for open source projects
            push!(plugins, TravisCI(;
                osx=false,
                windows=false,
                x64=false,
                arm64=true,
                coverage=false,
            ))
        end

        if !is_proprietary
            push!(plugins, Codecov())

            @info """Using coverage plugins, don't forget to manually add your API tokens as secrets,
            as described [here](https://help.github.com/en/actions/automating-your-workflow-with-github-actions/creating-and-using-encrypted-secrets#creating-encrypted-secrets)."""
        end
    end

    Template(;
        dir=parent_dir,
        plugins=plugins,
        julia=min_julia_version,
    )
end

end # module


text:
# global_variables.jl
@deprecate disable_combine_contract!() ITensors.disable_combine_contract()
@deprecate disable_tblis!() ITensors.disable_tblis()
@deprecate disable_warn_order!() ITensors.disable_warn_order()
@deprecate enable_combine_contract!() ITensors.enable_combine_contract()
@deprecate enable_tblis!() ITensors.enable_tblis()
@deprecate reset_warn_order!() ITensors.reset_warn_order()
@deprecate set_warn_order!(N) ITensors.set_warn_order(N)
@deprecate use_combine_contract() ITensors.using_combine_contract()
@deprecate use_debug_checks() ITensors.using_debug_checks()

# index.jl
@deprecate getindex(i::Index, n::Int) (i => n)

# indexset.jl
@deprecate store(is::IndexSet) data(is)
@deprecate firstintersect(is...; kwargs...) getfirst(intersect(is...); kwargs...)
@deprecate firstsetdiff(is...; kwargs...) getfirst(setdiff(is...); kwargs...)

# itensor.jl
@deprecate commonindex(args...; kwargs...) commonind(args...; kwargs...)
@deprecate emptyITensor(::Type{Any}) emptyITensor()
@deprecate findindex(args...; kwargs...) firstind(args...; kwargs...)
@deprecate findinds(args...; kwargs...) inds(args...; kwargs...)
@deprecate linkindex(args...; kwargs...) linkind(args...; kwargs...)
@deprecate matmul(A::ITensor, B::ITensor) product(A, B)
@deprecate replaceindex!(args...; kwargs...) replaceind!(args...; kwargs...)
@deprecate siteindex(args...; kwargs...) siteind(args...; kwargs...)
@deprecate store(A::ITensor) storage(A)
@deprecate setstore!(T::ITensor, st) setstorage!(T, st) false
@deprecate setstore(T::ITensor, st) setstorage(T, st) false
@deprecate uniqueindex(args...; kwargs...) uniqueind(args...; kwargs...)

# mps/abstractmps.jl
@deprecate orthoCenter(args...; kwargs...) orthocenter(args...; kwargs...)
@deprecate store(m::AbstractMPS) data(m)
@deprecate replacesites!(args...; kwargs...) ITensors.replace_siteinds!(args...; kwargs...)
@deprecate applyMPO(args...; kwargs...) contract(args...; kwargs...)
@deprecate applympo(args...; kwargs...) contract(args...; kwargs...)
@deprecate errorMPOprod(args...; kwargs...) error_contract(args...; kwargs...)
@deprecate error_mpoprod(args...; kwargs...) error_contract(args...; kwargs...)
@deprecate error_mul(args...; kwargs...) error_contract(args...; kwargs...)
@deprecate multMPO(args...; kwargs...) contract(args...; kwargs...)
@deprecate sum(A::AbstractMPS, B::AbstractMPS; kwargs...) add(A, B; kwargs...)
@deprecate multmpo(args...; kwargs...) contract(args...; kwargs...)
@deprecate set_leftlim!(args...; kwargs...) ITensors.setleftlim!(args...; kwargs...)
@deprecate set_rightlim!(args...; kwargs...) ITensors.setrightlim!(args...; kwargs...)
@deprecate tensors(args...; kwargs...) ITensors.data(args...; kwargs...)
@deprecate primelinks!(args...; kwargs...) ITensors.prime_linkinds!(args...; kwargs...)
@deprecate simlinks!(args...; kwargs...) ITensors.sim_linkinds!(args...; kwargs...)
@deprecate mul(A::AbstractMPS, B::AbstractMPS; kwargs...) contract(A, B; kwargs...)

# mps/mpo.jl
@deprecate MPO(A::MPS; kwargs...) outer(A, A; kwargs...)

# mps/mps.jl
@deprecate randomMPS(sites::Vector{<:Index}, linkdims::Integer) randomMPS(
  sites; linkdims=linkdims
)
@deprecate randomMPS(ElType::Type, sites::Vector{<:Index}, linkdims::Integer) randomMPS(
  ElType, sites; linkdims=linkdims
)
@deprecate randomMPS(sites::Vector{<:Index}, state, linkdims::Integer) randomMPS(
  sites, state; linkdims=linkdims
)

# physics/autompo.jl
@deprecate toMPO(args...; kwargs...) MPO(args...; kwargs...)

# qn/qn.jl
@deprecate store(qn::QN) data(qn)


text:
using StatsBase

#NETWORK PROPERTIES
LAYER_SIZE = 3
SEQUENCE_LEN = 25
LEARN_RATE = 1e-1

#NETWORK WEIGHTS, BIASES
function weights(vsize)
	W_XH = randn(LAYER_SIZE,vsize)*0.01
	W_HH = randn(LAYER_SIZE,LAYER_SIZE)*0.01
	W_HY = randn(vsize,LAYER_SIZE)*0.01
	return W_XH, W_HH, W_HY
end

function biases(vsize)
	B_H = zeros(LAYER_SIZE,1)
	B_Y = zeros(vsize,1)
	return B_H, B_Y
end

#LOSS FUNCTION

function loss(inputs,targets,vsize,W_XH,W_HH,W_HY,B_H,B_Y,LAST_H)
	xs, hs, ys, ps = Dict(), Dict(), Dict(), Dict()
	hs[0] = copy(LAST_H)
	error = 0

	for t in 1:length(inputs)
		xs[t] = zeros(vsize,1)
		xs[t][inputs[t]] = 1
		hs[t] = tanh((W_XH*xs[t]) + (W_HH*hs[t-1]) + B_H)
		ys[t] = (W_HY*hs[t]) + B_Y
		ps[t] = exp(ys[t])/sum(exp(ys[t]))
		error += -log(ps[t][targets[t],1])
	end

	D_W_XH, D_W_HH, D_W_HY = zeros(W_XH), zeros(W_HH), zeros(W_HY)
	D_B_H, D_B_Y = zeros(B_H), zeros(B_Y)
	H_NEXT = zeros(hs[0])

	for t in reverse(1:length(inputs))
		dy = copy(ps[t])
		dy[targets[t]] -= 1
		D_W_HY += dy*hs[t].'
		D_B_Y += dy
		dh = (W_HY.'*dy) + H_NEXT
		dh_raw = (1 - hs[t].*hs[t]).*dh
		D_B_H += dh_raw
		D_W_XH += dh_raw*xs[t].'
		D_W_HH += dh_raw*hs[t-1].'
		H_NEXT = W_HH'*dh_raw
	end

	D_W_XH = clamp!(D_W_XH,-5,5)
	D_W_HH = clamp!(D_W_HH,-5,5)
	D_W_HY = clamp!(D_W_HY,-5,5)
	D_B_H = clamp!(D_B_H,-5,5)
	D_B_Y = clamp!(D_B_Y,-5,5)

	return error, D_W_XH, D_W_HH, D_W_HY, D_B_H, D_B_Y, hs[length(inputs)]

end

#SAMPLER

function sample(h, seed, vsize, n, W_XH, W_HH, W_HY, B_H, B_Y)
	x = zeros(vsize,1)
	x[seed] = 1
	ixes = []

	for t in 1:n
		h = tanh(W_XH*x + W_HH*h + B_H)
		y = W_HY*h + B_Y
		p = exp(y)/sum(exp(y))
		w = StatsBase.weights(p)
		ix = StatsBase.sample(1:vsize,w)
		x = zeros(vsize,1)
		x[ix] = 1
		append!(ixes,ix)
	end

	return ixes
end

#WRITE TO FILE

function tofile(n, text)
	f = open("output.txt","a")
		write(f,"----------\n$n\n----------\n\n$text\n\n")
	close(f)
end

#INITIALIZE THE NETWORK
function run(contents)
	n,p = 1,1
	chars = Set([ch for ch in contents])
	vsize = length(chars)
	ix_char = Dict(ch for ch in enumerate(chars))
	char_ix = map(reverse,ix_char)
	errs = []

	smooth_loss = -log(1.0/vsize) * SEQUENCE_LEN

	W_XH, W_HH, W_HY = weights(vsize)
	B_H, B_Y = biases(vsize)

	mXh, mHh, mHy, mBh, mBy = zeros(W_XH), zeros(W_HH), zeros(W_HY), zeros(B_H), zeros(B_Y)

	while true
		if p + SEQUENCE_LEN + 1 >= length(contents) || n == 1
			LAST_H = zeros(LAYER_SIZE,1)
			p = 1
		end

		inputs = [char_ix[ch] for ch in contents[p:p+SEQUENCE_LEN]]
		targets = [char_ix[ch] for ch in contents[p+1:p+1+SEQUENCE_LEN]]

		error, D_W_XH, D_W_HH, D_W_HY, D_B_H, D_B_Y, LAST_H = loss(inputs, targets, vsize, W_XH, W_HH, W_HY, B_H, B_Y, LAST_H)

		append!(errs,error)

		if n%100 == 0
			println("$n: ", sum(errs)/length(errs))
		end

		if n%1000 == 0
			text = sample(LAST_H,inputs[1],vsize,200,W_XH,W_HH,W_HY,B_H,B_Y)
			tofile(n, join([ix_char[ch] for ch in text]))
		end

		#MANUAL WEIGHTING
		#XH
		mXh += D_W_XH.*D_W_XH
		W_XH += -LEARN_RATE * D_W_XH ./ sqrt(mXh + 1e-8)
		#HH
		mHh += D_W_HH.*D_W_HH
		W_HH += -LEARN_RATE * D_W_HH ./ sqrt(mHh + 1e-8)
		#HY
		mHy += D_W_HY.*D_W_HY
		W_HY += -LEARN_RATE * D_W_HY ./ sqrt(mHy + 1e-8)
		#BH
		mBh += D_B_H.*D_B_H
		B_H += -LEARN_RATE * D_B_H ./ sqrt(mBh + 1e-8)
		#BY
		mBy += D_B_Y.*D_B_Y
		B_Y += -LEARN_RATE * D_B_Y ./ sqrt(mBy + 1e-8)

		p += SEQUENCE_LEN
		n += 1
	end

end

#LOAD FILE

f = open("texts/gatsby_ord.txt")
	contents = readstring(f)
close(f)

y = run(contents)

max_stars_repo_name: dluman/char-rnn

max_issues_repo_name: dluman/char-rnn

max_forks_repo_name: dluman/char-rnn


text:
abstract type AbstractFEASolver end

abstract type AbstractDisplacementSolver <: AbstractFEASolver end

@params mutable struct DirectDisplacementSolver{T,dim,TP<:AbstractPenalty{T}} <:
                       AbstractDisplacementSolver
    problem::StiffnessTopOptProblem{dim,T}
    globalinfo::GlobalFEAInfo{T}
    elementinfo::ElementFEAInfo{dim,T}
    u::AbstractVector{T}
    lhs::AbstractVector{T}
    rhs::AbstractVector{T}
    vars::AbstractVector{T}
    penalty::TP
    prev_penalty::TP
    xmin::T
    qr::Bool
end
function Base.show(::IO, ::MIME{Symbol("text/plain")}, x::DirectDisplacementSolver)
    return println("TopOpt direct solver")
end
function DirectDisplacementSolver(
    sp::StiffnessTopOptProblem{dim,T};
    xmin=T(1) / 1000,
    penalty=PowerPenalty{T}(1),
    prev_penalty=deepcopy(penalty),
    quad_order=default_quad_order(sp),
    qr=false,
) where {dim,T}
    elementinfo = ElementFEAInfo(sp, quad_order, Val{:Static})
    globalinfo = GlobalFEAInfo(sp)
    u = zeros(T, ndofs(sp.ch.dh))
    lhs = similar(u)
    rhs = similar(u)
    vars = fill(one(T), getncells(sp.ch.dh.grid) - sum(sp.black) - sum(sp.white))
    varind = sp.varind
    return DirectDisplacementSolver(
        sp, globalinfo, elementinfo, u, lhs, rhs, vars, penalty, prev_penalty, xmin, qr
    )
end
function (s::DirectDisplacementSolver{T})(
    ::Type{Val{safe}}=Val{false},
    ::Type{newT}=T;
    assemble_f=true,
    reuse_chol=false,
    rhs=assemble_f ? s.globalinfo.f : s.rhs,
    lhs=assemble_f ? s.u : s.lhs,
    kwargs...,
) where {T,safe,newT}
    globalinfo = s.globalinfo
    N = size(globalinfo.K, 1)
    assemble!(
        globalinfo,
        s.problem,
        s.elementinfo,
        s.vars,
        getpenalty(s),
        s.xmin;
        assemble_f=assemble_f,
    )
    K = globalinfo.K
    if safe
        m = meandiag(K)
        for i in 1:size(K, 1)
            if K[i, i] ≈ zero(T)
                K[i, i] = m
            end
        end
    end
    nans = false
    if !reuse_chol
        try
            if T === newT
                if s.qr
                    globalinfo.qrK = qr(K.data)
                else
                    globalinfo.cholK = cholesky(Symmetric(K))
                end
            else
                if s.qr
                    globalinfo.qrK = qr((newT.(K)).data)
                else
                    globalinfo.cholK = cholesky(Symmetric(newT.(K)))
                end
            end
        catch err
            lhs .= T(NaN)
            nans = true
        end
    end
    if !nans
        if T === newT
            if s.qr
                lhs .= globalinfo.qrK \ rhs
            else
                lhs .= globalinfo.cholK \ rhs
            end
        else
            if s.qr
                lhs .= globalinfo.qrK \ newT.(rhs)
            else
                lhs .= globalinfo.cholK \ newT.(rhs)
            end
        end
    end
    return nothing
end


text:
@inline ∇_dot_qᶜ(i, j, k, grid::AbstractGrid{FT}, ::Nothing, args...) where FT = zero(FT)
@inline ∂ⱼ_τ₁ⱼ(i, j, k, grid::AbstractGrid{FT}, ::Nothing, args...) where FT = zero(FT)
@inline ∂ⱼ_τ₂ⱼ(i, j, k, grid::AbstractGrid{FT}, ::Nothing, args...) where FT = zero(FT)
@inline ∂ⱼ_τ₃ⱼ(i, j, k, grid::AbstractGrid{FT}, ::Nothing, args...) where FT = zero(FT)

calculate_diffusivities!(diffusivities, ::Nothing, args...) = nothing
calculate_diffusivities!(::Nothing, ::Nothing, args...) = nothing


text:
mutable struct ExternalWorld
    time::Int64
    function ExternalWorld()
        new(0)
    end
end

function update!(world)

end


text:
struct Peek{Rows}
    rows::Rows
    maximum_length::Int
end

export Peek
"""
    Peek(rows, maximum_length = 4)

Peek an iterator which returns named tuples. Will show no more than `maximum_length` rows.

```jldoctest Peek
julia> using LightQuery

julia> Peek(Rows(a = 1:5, b = 5:-1:1))
Showing 4 of 5 rows
|   a |   b |
| ---:| ---:|
|   1 |   5 |
|   2 |   4 |
|   3 |   3 |
|   4 |   2 |
```
"""
function Peek(rows)
    Peek(rows, 4)
end

function make_any(values)
    Any[values...]
end

function show(output::IO, peek::Peek)
    rows = peek.rows
    maximum_length = peek.maximum_length
    if isa(IteratorSize(rows), Union{HasLength,HasShape})
        if length(rows) > maximum_length
            println(output, "Showing $(maximum_length) of $(length(rows)) rows")
        end
    else
        println(output, "Showing at most $(maximum_length) rows")
    end
    columns = make_columns(take(rows, maximum_length))
    rows = map(make_any, zip(columns...))
    pushfirst!(rows, make_any(String.(propertynames(columns))))
    show(output, MD(Table(rows, make_any(map(function (column)
        :r
    end, columns)))))
end

"""
    function reduce_rows(rows, a_function, columns...)

Reduce a function over each of `columns` in `rows`.

```jldoctest
julia> using LightQuery

julia> using Test: @inferred

julia> @inferred reduce_rows(Rows(a = [1, 1], b = [1.0, 1.0]), +, name"a", name"b")
(a = 2, b = 2.0)
```
"""
function reduce_rows(rows, a_function, columns...)
    reduce(
        let a_function = a_function, columns = columns
            function (row1, row2)
                map(
                    a_function,
                    columns(row1),
                    columns(row2),
                )
            end
        end,
        rows,
    )
end

export reduce_rows


text:
# Automatically generated using Clang.jl wrap_c, version 0.0.0

using Compat

# Skipping MacroDefinition: InvalidOid ( ( Oid ) 0 )

const OID_MAX = UINT_MAX
const PG_DIAG_SEVERITY = 'S'
const PG_DIAG_SQLSTATE = 'C'
const PG_DIAG_MESSAGE_PRIMARY = 'M'
const PG_DIAG_MESSAGE_DETAIL = 'D'
const PG_DIAG_MESSAGE_HINT = 'H'
const PG_DIAG_STATEMENT_POSITION = 'P'
const PG_DIAG_INTERNAL_POSITION = 'p'
const PG_DIAG_INTERNAL_QUERY = 'q'
const PG_DIAG_CONTEXT = 'W'
const PG_DIAG_SCHEMA_NAME = 's'
const PG_DIAG_TABLE_NAME = 't'
const PG_DIAG_COLUMN_NAME = 'c'
const PG_DIAG_DATATYPE_NAME = 'd'
const PG_DIAG_CONSTRAINT_NAME = 'n'
const PG_DIAG_SOURCE_FILE = 'F'
const PG_DIAG_SOURCE_LINE = 'L'
const PG_DIAG_SOURCE_FUNCTION = 'R'
const PG_COPYRES_ATTRS = 0x01
const PG_COPYRES_TUPLES = 0x02
const PG_COPYRES_EVENTS = 0x04
const PG_COPYRES_NOTICEHOOKS = 0x08

# Skipping MacroDefinition: PQsetdb ( M_PGHOST , M_PGPORT , M_PGOPT , M_PGTTY , M_DBNAME ) PQsetdbLogin ( M_PGHOST , M_PGPORT , M_PGOPT , M_PGTTY , M_DBNAME , NULL , NULL )
# Skipping MacroDefinition: PQfreeNotify ( ptr ) PQfreemem ( ptr )

const PQnoPasswordSupplied = "fe_sendauth: no password supplied\n"

typealias Oid UInt32
typealias pg_int64 Clong

# begin enum ANONYMOUS_1
typealias ANONYMOUS_1 UInt32
const CONNECTION_OK = (UInt32)(0)
const CONNECTION_BAD = (UInt32)(1)
const CONNECTION_STARTED = (UInt32)(2)
const CONNECTION_MADE = (UInt32)(3)
const CONNECTION_AWAITING_RESPONSE = (UInt32)(4)
const CONNECTION_AUTH_OK = (UInt32)(5)
const CONNECTION_SETENV = (UInt32)(6)
const CONNECTION_SSL_STARTUP = (UInt32)(7)
const CONNECTION_NEEDED = (UInt32)(8)
# end enum ANONYMOUS_1

# begin enum ConnStatusType
typealias ConnStatusType UInt32
const CONNECTION_OK = (UInt32)(0)
const CONNECTION_BAD = (UInt32)(1)
const CONNECTION_STARTED = (UInt32)(2)
const CONNECTION_MADE = (UInt32)(3)
const CONNECTION_AWAITING_RESPONSE = (UInt32)(4)
const CONNECTION_AUTH_OK = (UInt32)(5)
const CONNECTION_SETENV = (UInt32)(6)
const CONNECTION_SSL_STARTUP = (UInt32)(7)
const CONNECTION_NEEDED = (UInt32)(8)
# end enum ConnStatusType

# begin enum ANONYMOUS_2
typealias ANONYMOUS_2 UInt32
const PGRES_POLLING_FAILED = (UInt32)(0)
const PGRES_POLLING_READING = (UInt32)(1)
const PGRES_POLLING_WRITING = (UInt32)(2)
const PGRES_POLLING_OK = (UInt32)(3)
const PGRES_POLLING_ACTIVE = (UInt32)(4)
# end enum ANONYMOUS_2

# begin enum PostgresPollingStatusType
typealias PostgresPollingStatusType UInt32
const PGRES_POLLING_FAILED = (UInt32)(0)
const PGRES_POLLING_READING = (UInt32)(1)
const PGRES_POLLING_WRITING = (UInt32)(2)
const PGRES_POLLING_OK = (UInt32)(3)
const PGRES_POLLING_ACTIVE = (UInt32)(4)
# end enum PostgresPollingStatusType

# begin enum ANONYMOUS_3
typealias ANONYMOUS_3 UInt32
const PGRES_EMPTY_QUERY = (UInt32)(0)
const PGRES_COMMAND_OK = (UInt32)(1)
const PGRES_TUPLES_OK = (UInt32)(2)
const PGRES_COPY_OUT = (UInt32)(3)
const PGRES_COPY_IN = (UInt32)(4)
const PGRES_BAD_RESPONSE = (UInt32)(5)
const PGRES_NONFATAL_ERROR = (UInt32)(6)
const PGRES_FATAL_ERROR = (UInt32)(7)
const PGRES_COPY_BOTH = (UInt32)(8)
const PGRES_SINGLE_TUPLE = (UInt32)(9)
# end enum ANONYMOUS_3

# begin enum ExecStatusType
typealias ExecStatusType UInt32
const PGRES_EMPTY_QUERY = (UInt32)(0)
const PGRES_COMMAND_OK = (UInt32)(1)
const PGRES_TUPLES_OK = (UInt32)(2)
const PGRES_COPY_OUT = (UInt32)(3)
const PGRES_COPY_IN = (UInt32)(4)
const PGRES_BAD_RESPONSE = (UInt32)(5)
const PGRES_NONFATAL_ERROR = (UInt32)(6)
const PGRES_FATAL_ERROR = (UInt32)(7)
const PGRES_COPY_BOTH = (UInt32)(8)
const PGRES_SINGLE_TUPLE = (UInt32)(9)
# end enum ExecStatusType

# begin enum ANONYMOUS_4
typealias ANONYMOUS_4 UInt32
const PQTRANS_IDLE = (UInt32)(0)
const PQTRANS_ACTIVE = (UInt32)(1)
const PQTRANS_INTRANS = (UInt32)(2)
const PQTRANS_INERROR = (UInt32)(3)
const PQTRANS_UNKNOWN = (UInt32)(4)
# end enum ANONYMOUS_4

# begin enum PGTransactionStatusType
typealias PGTransactionStatusType UInt32
const PQTRANS_IDLE = (UInt32)(0)
const PQTRANS_ACTIVE = (UInt32)(1)
const PQTRANS_INTRANS = (UInt32)(2)
const PQTRANS_INERROR = (UInt32)(3)
const PQTRANS_UNKNOWN = (UInt32)(4)
# end enum PGTransactionStatusType

# begin enum ANONYMOUS_5
typealias ANONYMOUS_5 UInt32
const PQERRORS_TERSE = (UInt32)(0)
const PQERRORS_DEFAULT = (UInt32)(1)
const PQERRORS_VERBOSE = (UInt32)(2)
# end enum ANONYMOUS_5

# begin enum PGVerbosity
typealias PGVerbosity UInt32
const PQERRORS_TERSE = (UInt32)(0)
const PQERRORS_DEFAULT = (UInt32)(1)
const PQERRORS_VERBOSE = (UInt32)(2)
# end enum PGVerbosity

# begin enum ANONYMOUS_6
typealias ANONYMOUS_6 UInt32
const PQPING_OK = (UInt32)(0)
const PQPING_REJECT = (UInt32)(1)
const PQPING_NO_RESPONSE = (UInt32)(2)
const PQPING_NO_ATTEMPT = (UInt32)(3)
# end enum ANONYMOUS_6

# begin enum PGPing
typealias PGPing UInt32
const PQPING_OK = (UInt32)(0)
const PQPING_REJECT = (UInt32)(1)
const PQPING_NO_RESPONSE = (UInt32)(2)
const PQPING_NO_ATTEMPT = (UInt32)(3)
# end enum PGPing

type pg_conn
end

typealias PGconn Void

type pg_result
end

typealias PGresult Void

type pg_cancel
end

typealias PGcancel Void

type pgNotify
    relname::Cstring
    be_pid::Cint
    extra::Cstring
    next::Ptr{pgNotify}
end

type PGnotify
    relname::Cstring
    be_pid::Cint
    extra::Cstring
    next::Ptr{pgNotify}
end

typealias PQnoticeReceiver Ptr{Void}
typealias PQnoticeProcessor Ptr{Void}
typealias pqbool UInt8

type _PQprintOpt
    header::pqbool
    align::pqbool
    standard::pqbool
    html3::pqbool
    expanded::pqbool
    pager::pqbool
    fieldSep::Cstring
    tableOpt::Cstring
    caption::Cstring
    fieldName::Ptr{Cstring}
end

type PQprintOpt
    header::pqbool
    align::pqbool
    standard::pqbool
    html3::pqbool
    expanded::pqbool
    pager::pqbool
    fieldSep::Cstring
    tableOpt::Cstring
    caption::Cstring
    fieldName::Ptr{Cstring}
end

type _PQconninfoOption
    keyword::Cstring
    envvar::Cstring
    compiled::Cstring
    val::Cstring
    label::Cstring
    dispchar::Cstring
    dispsize::Cint
end

type PQconninfoOption
    keyword::Cstring
    envvar::Cstring
    compiled::Cstring
    val::Cstring
    label::Cstring
    dispchar::Cstring
    dispsize::Cint
end

type PQArgBlock
    len::Cint
    isint::Cint
    u::Void
end

type pgresAttDesc
    name::Cstring
    tableid::Oid
    columnid::Cint
    format::Cint
    typid::Oid
    typlen::Cint
    atttypmod::Cint
end

type PGresAttDesc
    name::Cstring
    tableid::Oid
    columnid::Cint
    format::Cint
    typid::Oid
    typlen::Cint
    atttypmod::Cint
end

typealias pgthreadlock_t Ptr{Void}

max_stars_repo_name: rbt-lang/rbt-proto

max_issues_repo_name: rbt-lang/rbt-proto

max_forks_repo_name: rbt-lang/rbt-proto


text:
function g2_semiinf(τ, B, ρ; nmed = 1.33, ndet = 1.5, μsp = 10.0, μa = 0.1, λ = 700)

    αDB::Float64 = B[1]
    β::Float64 = B[2]
    n::Float64 = nmed/ndet

    D::Float64 = 1/3μsp

    if n > 1.0
        A = 504.332889 - 2641.00214n + 5923.699064n^2 - 7376.355814n^3 +
        5507.53041n^4 - 2463.357945n^5 + 610.956547n^6 - 64.8047n^7
    elseif n < 1.0
        A = 3.084635 - 6.531194n + 8.357854n^2 - 5.082751n^3
    else
        A = 1.0
    end

    z0::Float64 = 1/(μsp + μa)
    zb::Float64 = 2A*D
    k0::Float64 = (n*2*pi)/(λ*10^-7)

    r1 = sqrt(ρ^2 +z0^2)
    r2 = sqrt(ρ^2 + (z0 +2zb)^2)

    C1 = 3*μsp*μa
    C2 = 6*μsp^2*k0^2*αDB

    A = zeros(length(τ))

    Threads.@threads for n in eachindex(τ)
        Kτ = sqrt(C1 + C2*τ[n])

        A[n] = exp(-r1*Kτ)/r1 - exp(-r2.*Kτ)/r2
    end

    g1 = (A)./(A[1])
    g2 = 1 .+ β*g1.^2
end

function test(τ, B::Array{Float64,1}, ρ::Float64; nmed::Float64 = 1.33, ndet::Float64 = 1.5, μsp::Float64 = 10.0, μa::Float64 = 0.1, λ::Float64 = 700.0)

    αDB::Float64 = B[1]
    β::Float64 = B[2]
    n::Float64 = nmed/ndet

    D::Float64 = 1/3μsp

    if n > 1.0
        A = 504.332889 - 2641.00214n + 5923.699064n^2 - 7376.355814n^3 +
        5507.53041n^4 - 2463.357945n^5 + 610.956547n^6 - 64.8047n^7
    elseif n < 1.0
        A = 3.084635 - 6.531194n + 8.357854n^2 - 5.082751n^3
    else
        A = 1.0
    end

    z0::Float64 = 1/(μsp + μa)
    zb::Float64 = 2A*D
    k0::Float64 = (n*2*pi)/(λ*10^-7)

    r1 = sqrt(ρ^2 +z0^2)
    r2 = sqrt(ρ^2 + (z0 +2zb)^2)

    C1 = 3*μsp*μa
    C2 = 6*μsp^2*k0^2*αDB

    A = zeros(length(τ))

    Threads.@threads for n in eachindex(τ)
        Kτ = sqrt(C1 + C2*τ[n])

        A[n] = exp(-r1*Kτ)/r1 - exp(-r2.*Kτ)/r2
    end

    g1 = (A)./(A[1])
    g2 = 1 .+ β*g1.^2
end

beta = 0.5
BFi = 2e-8
τ = 1e-8:1e-6:1e-2
test(τ, [BFi, beta], 1.)

t = 0:0.01:10

ydata =  test(τ, [BFi, beta], 1.) + 0.05*randn(length(τ))

function curvetestff(ydata, τ)

    fit = curve_fit((τ, β) -> test(τ, β, 1.0), τ, ydata, [5e-9, 0.3]; autodiff=:finiteforward)
end

#=
@benchmark test(τ, [BFi, beta], 1.)
BenchmarkTools.Trial:
  memory estimate:  704.81 KiB
  allocs estimate:  19576
  --------------
  minimum time:     96.266 μs (0.00% GC)
  median time:      101.922 μs (0.00% GC)
  mean time:        146.992 μs (23.95% GC)
  maximum time:     8.491 ms (97.80% GC)
  --------------
  samples:          10000
  evals/sample:     1

max_stars_repo_name: sammzer/LightPropagation

max_issues_repo_name: sammzer/LightPropagation

max_forks_repo_name: sammzer/LightPropagation


text:
module TreatmentPanels

include("TreatmentPanel.jl")
include("show_and_plot.jl")

# Export main types
export TreatmentPanel
export BalancedPanel, UnbalancedPanel

# Export treatment description types
export UnitTreatmentType
export SingleUnitTreatment, MultiUnitSimultaneousTreatment, MultiUnitStaggeredTreatment
export TreatmentDurationType
export ContinuousTreatment, StartEndTreatment

# Export utility functions
export treated_ids, treated_labels, first_treated_period_ids, first_treated_period_labels, length_T₀, length_T₁

end


text:
# Experimental AD codes

import ChainRulesCore, ChainRules
import ChainRulesCore: rrule, NoTangent

struct EVAL{TM, TX}
   m::TM
   X::TX
end

struct EVAL_D{TM, TX}
   m::TM
   X::TX
end

function (y::EVAL)(params)
   set_params!(y.m, params)
   return evaluate(y.m, y.X).val
end

function rrule(y::EVAL, params)
   set_params!(y.m, params)
   val = evaluate(y.m, y.X).val
   adj = dp -> ( NoTangent(), dp * getproperty.(grad_params(y.m, y.X), :val))
   return val, adj
end

function (y::EVAL_D)(params)
   set_params!(y.m, params)
   return grad_config(y.m, y.X)
end

function rrule(y::EVAL_D, params)
   set_params!(y.m, params)
   val = grad_config(y.m, y.X)
   adj = dp -> ( NoTangent(), adjoint_EVAL_D(y.m, y.X, dp) )
   return val, adj
end


text:
const REG_MODELS =  (
    SigmaRidgeRegressor,
    SingleGroupRidgeRegressor,
    MultiGroupRidgeRegressor,
    LooRidgeRegressor,
)

_path(M) =  "SigmaRidgeRegression.$(MMI.name(M))"
_descr(::Type{<:SigmaRidgeRegressor}) = "σ-Ridge regression"
_descr(::Type{<:SingleGroupRidgeRegressor}) = "Vanilla ridge regression"
_descr(::Type{<:MultiGroupRidgeRegressor}) = "Group-regularized ridge regression"
_descr(::Type{<:LooRidgeRegressor}) = "Ridge regressor tuned by leave-one-out cross-validation"

MMI.metadata_pkg.(
    REG_MODELS,
    package_name = "SigmaRidgeRegression",
    package_uuid = "e566bdc2-f852-4ed0-870a-7389015b131d",
    package_url = "https://github.com/nignatiadis/SigmaRidgeRegression.jl",
    is_pure_julia = true,
    package_license = "MIT",
    is_wrapper = false,
)

for M in REG_MODELS
    MMI.metadata_model(M,
        input_scitype=MMI.Table(MMI.Continuous),
        target_scitype=AbstractVector{MMI.Continuous},
        supports_weights=false,
        docstring=_descr(M),
        load_path=_path(M))
end


text:
# Worked with Kasra and Waleed #

cd("D:\\Semester 3_Fall 2020\\ECON 6453 - Econometrics III\\ProblemSets\\PS4-mixture")
using Random, Distributions, JLD2, DataFrames, CSV, LinearAlgebra, Statistics, FreqTables, TexTables, Optim, HTTP, GLM, ForwardDiff

#:::::::::::::::::::::::::::::::::::::::::::::::::::
# Question 1: Multinomial logit with Alternate Specific Covariates
#:::::::::::::::::::::::::::::::::::::::::::::::::::

url = "https://raw.githubusercontent.com/OU-PhD-Econometrics/fall-2020/master/ProblemSets/PS4-mixture/nlsw88t.csv"
df = CSV.read(HTTP.get(url).body)

X = [df.age df.white df.collgrad]
Z = hcat(df.elnwage1, df.elnwage2, df.elnwage3, df.elnwage4, df.elnwage5, df.elnwage6, df.elnwage7, df.elnwage8)
y = df.occ_code
# Using PS3's code for Multinomial Logit #
function mlogit(alpha, X, Z, y)
    beta = alpha[1:end-1] # Indexing the parameter so coefficient on X is the first element
    gamma = alpha[end] # Indexing the parameter so coefficient on Z is the last element
    K = size(X,2) # K is the no. of covariates in X, which in this case are 3
    J = length(unique(y)) # J is the no. of coefficients, which are 8 in this case (So J-1 = 7)
    N = length(y) # N is a vector with the same length as y (not sure why we are defining N)
    bigY = zeros(N,J) # Simply defining the matrix of dependent variables Y, which are going to be dummies in case of Multinomial Logit

    for j=1:J # Simply defining j to be something, next step, j represents column
        bigY[:,j] = y.==j # Making a dummies for each of the corresponding occupation in y
    end
    bigBeta = [reshape(beta,K,J-1) zeros(K)] # Reshaping into a matrix with K rows and J-1 columns
    # Adding a column of zeros in the end to make it conform with bigY and also to normalize the parameters of the last occupation to be zero

    T = promote_type(eltype(X),eltype(alpha)) # To get the standard errors we promote the type
    num   = zeros(T,N,J) # Define numerator T-type , N-rows, J-columns
    dem   = zeros(T,N) # Define denominator T-type , N-rows
    for j=1:J
        num[:,j] = exp.(X*bigBeta[:,j] .+ (Z[:,j] .- Z[:,J])*gamma) # Numerator Likelihood function
        dem .+= num[:,j] # Denominator Likelihood function: Simply the Sum of all the numerators
    end

    P = num./repeat(dem,1,J) # Putting together numerator and denominator

    loglike = -sum( bigY.*log.(P) ) # Since Optim minimizes the function, we take its negative to maximize

    return loglike # This is used for storing the function parameters
end

startvals = [2*rand(7*size(X,2)).-1; .1]
td = TwiceDifferentiable(alpha -> mlogit(alpha, X, Z, y), startvals; autodiff = :forward)

# Use Optim #
α_hat_multinomial = optimize(td, startvals, LBFGS(), Optim.Options(g_tol = 1e-5, iterations=100_000, show_trace=true, show_every=50))
α_hat_multinomial_ad = α_hat_multinomial.minimizer

# Find Hessian #
H  = Optim.hessian!(td, α_hat_multinomial_ad) # Finding the Hessian Matrix
α_hat_multinomial_ad_se = sqrt.(diag(inv(H))) # Standard Errors of alpha_hat

println([α_hat_multinomial_ad α_hat_multinomial_ad_se])

#:::::::::::::::::::::::::::::::::::::::::::::::::::
# Question 2: Interpretation of gamma_hat
#:::::::::::::::::::::::::::::::::::::::::::::::::::
# gamma_hat make more sense as compared to PS3 because it's sign is positive #

#:::::::::::::::::::::::::::::::::::::::::::::::::::
# Question 3:
#:::::::::::::::::::::::::::::::::::::::::::::::::::

### Part (A) ###
using Distributions
include("lgwt.jl") # make sure the function gets read in

# define distribution
d = Normal(0,1) # mean=0, standard deviation=1

# get quadrature nodes and weights for 7 grid points, with +_ 4SD
nodes, weights = lgwt(7,-4,4)

# now compute the integral over the density and verify it's 1
sum(weights.*pdf.(d,nodes)) # d here is normal dist. as defined earlier #

# now compute the expectation and verify it's 0
sum(weights.*nodes.*pdf.(d,nodes)) # it seems as if the nodes here is x in: intergral x phi(x) dx #

### Part (B) ###

#Practice 1#
d = Normal(0,2)
nodes, weights = lgwt(7,-5,5)
sum(weights.*(nodes.^2).*pdf.(d,nodes))

#Practice 2#
nodes, weights = lgwt(10,-5,5)
sum(weights.*(nodes.^2).*pdf.(d,nodes))

# Comment: Increasing the quadrature points gives us better estimates (of variace of f(x)) #

### Part (C) ###

#Practice 1#
Random.seed!(1234)
d = Normal(0,2)
N = 1000000
nodes, weights = lgwt(10,-5,5)
x = rand(Uniform(-10,10), N)
A = (10-(-10))*(1/N)*sum((x.^2).*pdf(d,x)) # A is approximately equal to 4 #

#Practice 2#
A1 = (10-(-10))*(1/N)*sum((x).*pdf(d,x)) # A1 is approx. equal to 0 #
#Practice 3#
A2 = (10-(-10))*(1/N)*sum(pdf(d,x)) # A2 is approx. equal to 1 #

#Practice 4#
D = 1000
x = rand(Uniform(-10,10), D)

A = (10-(-10))*(1/D)*sum((x.^2).*pdf(d,x))

A1 = (10-(-10))*(1/D)*sum((x).*pdf(d,x))

A2 = (10-(-10))*(1/D)*sum(pdf(d,x))
# Reducing D to 1000 more likely less precise estimates #

### Part (D) ###

#:::::::::::::::::::::::::::::::::::::::::::::::::::
# Question 4: With Gauss-Lagendre Quadrature
#:::::::::::::::::::::::::::::::::::::::::::::::::::
function mixedlogit_quadrature(theta, X, Z, y, rn)

        alpha = theta[1:end-1]
        gamma = theta[end]
        K = size(X,2)
        J = length(unique(y))
        N = length(y)
        bigY = zeros(N,J)
        for j=1:J
            bigY[:,j] = y.==j
        end
        bigAlpha = [reshape(alpha,K,J-1) zeros(K)]

        nodes, weights = lgwt(rn,-4,4);
        d=Normal(0,1)

        T = promote_type(eltype(X),eltype(theta))
        num   = zeros(T,N,J)
        dem   = zeros(T,N)
        for j=1:J
            num[:,j] = sum(exp.(X*bigAlpha[:,j] .+ (Z[:,j] .- Z[:,J])*gamma).*weights.*pdf.(d,nodes))
            dem .+= num[:,j]
        end

        P = num./repeat(dem,1,J)

        loglike = -sum( log.(P.^bigY) )

        return loglike
    end

#:::::::::::::::::::::::::::::::::::::::::::::::::::
# Question 5: With Monte Carlo Simulations
#:::::::::::::::::::::::::::::::::::::::::::::::::::
function mixedlogit_simulation(theta, X, Z, y,rn)

        alpha = theta[1:end-1]
        gamma = theta[end]
        K = size(X,2)
        J = length(unique(y))
        N = length(y)
        bigY = zeros(N,J)
        for j=1:J
            bigY[:,j] = y.==j
        end
        bigAlpha = [reshape(alpha,K,J-1) zeros(K)]

        d=Normal(0,1)
        rx=rand(Uniform(0,1 ),rn)

        T = promote_type(eltype(X),eltype(theta))
        num   = zeros(T,N,J)
        dem   = zeros(T,N)
        for j=1:J
            num[:,j] = sum(exp.(X*bigAlpha[:,j] .+ (Z[:,j] .- Z[:,J])*gamma).*(pdf.(d,x)))
            dem .+= num[:,j]
        end

        P = num./repeat(dem,1,J)

        loglike = -sum( log.(P.^bigY) )

        return loglike
    end

                                ### End ###


text:
"An abstract type for the problems from [NSDE.jl](https://github.com/antonuccig/NSDE.jl)."
abstract type AbstractNSDEProblem end

"An abstract type for the solvers from [NSDE.jl](https://github.com/antonuccig/NSDE.jl)."
abstract type AbstractNSDESolver end

"An abstract type for the solutions from [NSDE.jl](https://github.com/antonuccig/NSDE.jl)."
abstract type AbstractNSDESolution end

"An abstract type for initial value problems."
abstract type AbstractInitialValueProblem <: AbstractNSDEProblem end

"An abstract type for the solvers of [`AbstractInitialValueProblem`](@ref)s."
abstract type AbstractInitialValueSolver <: AbstractNSDESolver end

"An abstract type for the solution of an [`AbstractInitialValueProblem`](@ref) obtained with an [`AbstractInitialValueSolver`](@ref)."
abstract type AbstractInitialValueSolution <: AbstractNSDESolution end

"An abstract type for the right-hand side of an [`AbstractInitialValueProblem`](@ref)."
abstract type AbstractRightHandSideFunction end