text:
module PovertyAndInequalityMeasures

using Base: Real
using IterableTables
using IteratorInterfaceExtensions
using TableTraits
using DataFrames

export DEFAULT_ATKINSON_ES, DEFAULT_ENTROPIES, DEFAULT_FGT_ALPHAS
export make_gini, make_poverty, make_inequality, binify, add_decomposed_theil
export PovertyMeasures, InequalityMeasures, ineqs_equal, povs_equal

const WEIGHT          = 1
const INCOME          = 2
const WEIGHTED_INCOME = 3
const POPN_ACCUM      = 4
const INCOME_ACCUM    = 5
const DEFAULT_FGT_ALPHAS = [ 0.0, 0.50, 1.0, 1.50, 2.0, 2.5 ];
const DEFAULT_ATKINSON_ES = [ 0.25, 0.50, 0.75, 1.0, 1.25, 1.50, 1.75, 2.0, 2.25 ];
const DEFAULT_ENTROPIES = [ 1.25, 1.50, 1.75, 2.0, 2.25, 2.5 ];

mutable struct PovertyMeasures{T<:Real}
    headcount :: T
    gap :: T
    watts :: T
    sen :: T
    shorrocks :: T
    fgt_alphas :: Vector{T}
    foster_greer_thorndyke :: Vector{T}
    time_to_exit :: T
    gini_amongst_poor :: T
    poverty_gap_gini :: T
end

function povs_equal( p1 :: PovertyMeasures, p2 :: PovertyMeasures ) :: Bool
    return (p1.headcount ≈ p2.headcount) &&
        (p1.gap ≈ p2.gap) &&
        (p1.watts ≈ p2.watts) &&
        (p1.sen ≈ p2.sen) &&
        (p1.shorrocks ≈ p2.shorrocks) &&
        (p1.fgt_alphas ≈ p2.fgt_alphas) &&
        (p1.foster_greer_thorndyke ≈ p2.foster_greer_thorndyke) &&
        (p1.time_to_exit ≈ p2.time_to_exit) &&
        (p1.gini_amongst_poor ≈ p2.gini_amongst_poor) &&
        (p1.poverty_gap_gini ≈ p2.poverty_gap_gini)
end

mutable struct InequalityMeasures{T<:Real}
    atkinson_es :: Vector{T}
    atkinson :: Vector{T}
    generalised_entropy_alphas :: Vector{T}
    generalised_entropy :: Vector{T}
    hoover :: T
    theil_l :: T
    theil_t :: T

    gini :: T
    palma :: T
    median :: T
    total_income :: T
    average_income :: T
    total_population :: T
    deciles :: Matrix{T}
    negative_or_zero_income_count :: Int
end

function ineqs_equal( i1 :: InequalityMeasures, i2 :: InequalityMeasures; include_populations :: Bool = true ) :: Bool
    eq = ( i1.atkinson_es ≈ i2.atkinson_es ) &&
        ( i1.atkinson ≈ i2.atkinson ) &&
        ( i1.generalised_entropy_alphas ≈ i2.generalised_entropy_alphas ) &&
        ( i1.generalised_entropy ≈ i2.generalised_entropy ) &&
        ( i1.hoover ≈ i2.hoover ) &&
        ( i1.theil_l ≈ i2.theil_l ) &&
        ( i1.theil_t ≈ i2.theil_t ) &&
        ( i1.gini ≈ i2.gini ) &&
        ( i1.palma ≈ i2.palma ) &&
        ( i1.median ≈ i2.median )

    if include_populations
        eq2 = ( i1.total_income ≈ i2.total_income ) &&
        ( i1.average_income ≈ i2.average_income ) &&
        ( i1.total_population ≈ i2.total_population ) &&
        ( i1.deciles ≈ i2.deciles )
        eq = eq && eq2
    end
    return eq
end

#
# FIXME this shouldn't need to return a copy!
#
function sortAndAccumulate(
    aug :: Matrix, sort_data :: Bool, nrows :: Integer )
    if sort_data
        aug = sortslices( aug, alg=QuickSort, dims=1,lt=((x,y)->isless(x[INCOME],y[INCOME])))
    end
    cumulative_weight = 0.0
    cumulative_income = 0.0
    for row in 1:nrows
            cumulative_weight += aug[row,WEIGHT]
            cumulative_income += aug[row,WEIGHTED_INCOME]
            aug[row,POPN_ACCUM] = cumulative_weight
            aug[row,INCOME_ACCUM] = cumulative_income
    end
    return aug
end # sortAndAccumulate

"
internal function that makes a sorted array
with cumulative income and population added
"
function make_augmented(
    data,
    weightcol :: Symbol,
    incomecol :: Symbol
    ;
    sort_data  :: Bool = true,
    delete_negatives :: Bool = false ) :: Matrix
    @assert TableTraits.isiterabletable( data ) "data needs to implement IterableTables"
    data = DataFrame(data) # this just makes iteration&handling missings easier
    T = eltype( data[!,incomecol] )
    if ! ( T <: AbstractFloat )
	    T = Float64 # jam on some real type regardless of the datatype of the frame if not some type of real
    end
	# better if int then float otherwise leave
    #
    # iter = IteratorInterfaceExtensions.getiterator(data)
    nrows = size(data)[1]
    aug = zeros( T, nrows, 5 )
    r = 0
    for row in eachrow(data)
        if ( ismissing( row[weightcol])) || ( ismissing( row[incomecol])) ||
            ( delete_negatives && row[incomecol] < 0.0 )
            ;
        else
            r += 1
            aug[r,WEIGHT] = T(row[weightcol])
            aug[r,INCOME] = T(row[incomecol])
            aug[r,WEIGHTED_INCOME] = aug[r,WEIGHT]*aug[r,INCOME]
        end # not missing or negative
    end
    aug = aug[1:r,:]
    return sortAndAccumulate( aug, sort_data, r )
end

"
internal function that makes a sorted array
with cumulative income and population added
"
function make_augmented(
    data            :: Matrix,
    weightpos       :: Integer = 1,
    incomepos       :: Integer = 2;

    sort_data        :: Bool = true,
    delete_negatives :: Bool = false ) :: Matrix
    T = eltype( data )
    if ! ( T <: AbstractFloat )
	    T = Float64 # jam on some real type regardless of the datatype of the data matrix if not some type of real
    end

    nrows = size( data )[1]
    aug = zeros( T, nrows, 5 )
    r = 0
    for row in 1:nrows
        if( delete_negatives && data[row,incomepos] < 0.0 )
            ;
        else
            r += 1
            aug[r,WEIGHT] = data[row,weightpos]
            aug[r,INCOME] = data[row,incomepos]
            aug[r,WEIGHTED_INCOME] = data[row,incomepos]*data[row,weightpos]
        end # not neg or negs accepted
    end # row loop
    return sortAndAccumulate( aug[1:r,:], sort_data, r )
    # return aug[1:r,:]
end

"""
calculate a Gini coefficient on one of our sorted arrays
"""
function make_gini( data :: Matrix ) :: Real
    lorenz = 0.0
    nrows = size( data )[1]
    if nrows == 0
        return 0.0
    end
    lastr = data[nrows,:]
    for row in 1:nrows
        lorenz += data[row, WEIGHT]*((2.0*data[row,INCOME_ACCUM]) - data[row,WEIGHTED_INCOME])
    end
    return 1.0-(lorenz/lastr[INCOME_ACCUM])/lastr[POPN_ACCUM]
end

"""
generate a subset of one of our datasets with just the elements whose incomes
are below the line. Probably possible in 1 line, once I get the hang of this
a bit more.
"""
function make_all_below_line( data :: Matrix, line :: Real ) :: Matrix
    nrows = size( data )[1]
    ncols = size( data )[2]
    T = eltype( data )
    outa = zeros( T, nrows, ncols )
    @assert ncols == 5 "data should have 5 cols"
    nout = 0
    for row in 1:nrows
        if data[row,INCOME] < line
            nout += 1
            outa[ nout, : ] .= data[row,:]
        end
    end
    return outa[1:nout,:]
end

"""
Create a dictionary of poverty measures.

This is based on the [World Bank's Poverty Handbook](http://documents.worldbank.org/curated/en/488081468157174849/Handbook-on-poverty-and-inequality)
by  Haughton and Khandker.

Arguments:
* `rawdata` - an nxm array of real nunbers; each row is an observation; one col should be a weight, another is income;
   positions assumed to be 1 and 2 unless weight and incomepos are supplied
* `line` - a poverty line, assumed same for all obs (so income is assumed to be already equivalised)
* `foster_greer_thorndyke_alphas` - coefficients for FGT poverty measures; note that FGT(0)
   corresponds to headcount and FGT(1) to gap; count and gap are computed directly anyway
   but it's worth checking one against the other.
* `growth` is (e.g.) 0.01 for 1% per period, and is used for 'time to exit' measure.

Output is  dictionary with an entry for each measure.

The measures are:

* `headcount`;
* `gap`;
* `Foster Greer Thorndyke`, for each of the values in `foster_greer_thorndyke_alphas`;
* `Watts`;
* `time to exit`, for the supplied growth rate;
* `Shorrocks`;
* `Sen`.

See World Bank, ch. 4.

"""
function make_poverty(
    rawdata                       :: Matrix,
    line                          :: Real,
    growth                        :: Real,
    weightpos                     :: Integer = 1,
    incomepos                     :: Integer = 2,
    foster_greer_thorndyke_alphas :: Vector = DEFAULT_FGT_ALPHAS ) :: PovertyMeasures

    data = make_augmented( rawdata, weightpos, incomepos )
    make_povertyinternal(
        data = data,
        line = line,
        growth = growth,
        foster_greer_thorndyke_alphas = foster_greer_thorndyke_alphas )
end

"""
As above, but using the QueryVerse IterableTables interface
rawdata - basically anything resembling a dataframe; see: [https://github.com/queryverse/IterableTables.jl]
throws an exception if rawdata doesn't support iterabletable interface
"""
function make_poverty(
    rawdata,
    line                          :: Real,
    growth                        :: Real,
    weightcol                     :: Symbol,
    incomecol                     :: Symbol,
    foster_greer_thorndyke_alphas :: Vector = DEFAULT_FGT_ALPHAS,
     ) :: PovertyMeasures
    @assert TableTraits.isiterabletable( rawdata ) "data needs to implement IterableTables"
    data = make_augmented( rawdata, weightcol, incomecol )
    make_povertyinternal(
        data = data,
        line = line,
        growth = growth,
        foster_greer_thorndyke_alphas = foster_greer_thorndyke_alphas )
end

function calc_positive_inc_popn( data :: Matrix ) :: Real
    pinc = 0.0
    nrows = size( data )[1]
    for r in 1:nrows
        if data[INCOME] > 0
            pinc += data[r,WEIGHT]
        end
    end
    return pinc;
end # postive_inc_population

"""
Internal version, once we have our datatset
"""
function make_povertyinternal(
    ;
    data                          :: Matrix,
    line                          :: Real,
    growth                        :: Real = 0.0,
    foster_greer_thorndyke_alphas :: Vector = DEFAULT_FGT_ALPHAS ) :: PovertyMeasures

    nfgs = size( foster_greer_thorndyke_alphas )[1]
    T = eltype( foster_greer_thorndyke_alphas )
    z = zero(T)
    pv = PovertyMeasures(
        z,
        z,
        z,
        z,
        z,
        foster_greer_thorndyke_alphas,
        zeros( nfgs ),
        z,
        z,
        z )

    nrows = size( data )[1]
    ncols = size( data )[2]
    population = data[ nrows, POPN_ACCUM ]
    total_income = data[ nrows, INCOME_ACCUM ]
    positive_inc_popn = calc_positive_inc_popn( data )

    @assert ncols == 5 "data should have 5 cols"

    belowline = make_all_below_line( data, line )
    nbrrows = size( belowline )[1]

    pv.gini_amongst_poor = make_gini( belowline )
    for row in 1:nbrrows
        inc :: T= belowline[row,INCOME]
        weight :: T = belowline[row,WEIGHT]
        gap :: T = line - inc
        @assert gap >= 0 "poverty gap must be postive"
        pv.headcount  += weight
        pv.gap  += weight*gap/line
        if belowline[row,INCOME ] > 0
            pv.watts  += weight*log(line/inc)
        end
        for p in 1:nfgs
            fg = foster_greer_thorndyke_alphas[p]
            pv.foster_greer_thorndyke[p] += weight*((gap/line)^fg)
        end
    end # main loop
    pv.watts /= population
    if growth > 0.0
        pv.time_to_exit = pv.watts/growth
    end
    pv.gap /= population
    pv.headcount /= population
    pv.foster_greer_thorndyke ./= population
    #
    # Gini of poverty gaps; see: WB pp 74-5
    #
    # create a 'Gini of the Gaps'
    # the sort routine in make_augmented does a really
    # bad job here either because the data
    # is mostly zeros or because it's reverse sorted
    # (smallest income -> biggest gap)
    # we we can just create the dataset in reverse
    # and use that
    gdata = zeros( T, nrows, 5 )
    for row in 1:nrows
        gap = max( 0.0, line - data[row,INCOME] )
        gpos = nrows - row + 1
        gdata[gpos,INCOME] = gap;
        gdata[gpos,WEIGHT] = data[row,WEIGHT]
    end
    gdata = make_augmented( gdata, 1, 2, sort_data=false )
    pv.poverty_gap_gini = make_gini( gdata )

    pv.sen = pv.headcount*pv.gini_amongst_poor+pv.gap*(1.0-pv.gini_amongst_poor)
    pv.shorrocks  = pv.headcount*pv.gap*(1.0+pv.poverty_gap_gini)
    return pv
end # make_poverty

"""
Make a wee dict with :theil_between and :theil_within
See WB eqns 6.7/6.8.
TODO
1. there are some papers on decomposing Atkinson, but I
don't understand them ..
2. the over time 3-part version

 popindc : Inequal for the population as a whole
 subindices : an array of dics, one for each subgroup of interest
"""
function add_decomposed_theil( popindc :: InequalityMeasures, subindices :: Vector{InequalityMeasures} ) :: NamedTuple
    popn = popindc.total_population
    income = popindc.total_income
    avinc = popindc.average_income
    within_l = 0.0
    within_t = 0.0

    between_l = 0.0
    between_t = 0.0
    totalpop = 0.0
    totalinc = 0.0
    for ind in subindices
        popshare = ind.total_population/popn
        incshare = ind.total_income/income
        totalpop += popshare
        totalinc += incshare

        within_l += ind.theil_l*popshare
        between_l += popshare*log(avinc/ind.average_income)

        within_t += ind.theil_t*incshare
        between_t += incshare*log(incshare/popshare)
    end
    @assert totalpop ≈ 1.0
    @assert totalinc ≈ 1.0
    @assert within_l+between_l ≈ popindc.theil_l  "within_l=$(within_l) between_l=$(between_l) ≈ popindc.theil_l=$(popindc.theil_l) "
    @assert within_t+between_t ≈ popindc.theil_t "within_t=$(within_t) between_t=$(between_t) ≈ popindc.theil_t=$(popindc.theil_t) "
    ( between_l = between_l, within_l = within_l,  between_t = between_t, within_t = within_t)
end

"""
Make a struct of inequality measures.
This is mainly taken from chs 5 and 6 of the World Bank book.

1. `rawdata` a matrix with cols with weights and incomes
2. `atkinson_es` inequality aversion values for the Atkinson indexes
3. `generalised_entropy_alphas`
4. `weightpos` - column with weights
5. `incomepos` - column with incomes

Returned is a Dict of inequality measures with:

* `Gini`;
* `Atkinson`, for each value in `atkinson_es`;
* `Theil`;
* `generalised_entropy`;
* `Hoover`;
* `Theil`;
* `Palma`.

See WB chs 5 an 6, and Cobham and Sumner on the Palma.

Also in the struct are:

* `total_income`
* `total_population`
* `average_income`
* `deciles`.

"""
function make_inequality(
    rawdata                    :: Matrix,
    weightpos                  :: Integer = 1,
    incomepos                  :: Integer = 2,
    atkinson_es                :: Vector{<:AbstractFloat} = DEFAULT_ATKINSON_ES,
    generalised_entropy_alphas :: Vector{<:AbstractFloat} = DEFAULT_ENTROPIES ) :: InequalityMeasures
    data = make_augmented( rawdata, weightpos, incomepos )
    return make_inequalityinternal(
        data = data,
        atkinson_es = atkinson_es,
        generalised_entropy_alphas = generalised_entropy_alphas
    )
end

"""
As above but using the iterable table interface; see: https://github.com/queryverse/IterableTables.jl
weightcol and incomecol are names of columns
"""
function make_inequality(
    rawdata,
    weightcol                  :: Symbol,
    incomecol                  :: Symbol,
    atkinson_es                :: Vector{<:AbstractFloat} = DEFAULT_ATKINSON_ES,
    generalised_entropy_alphas :: Vector{<:AbstractFloat} = DEFAULT_ENTROPIES ) :: InequalityMeasures
    @assert TableTraits.isiterabletable( rawdata ) "data needs to implement IterableTables"
    data = make_augmented( rawdata, weightcol, incomecol )
    return make_inequalityinternal(
        data = data,
        atkinson_es = atkinson_es,
        generalised_entropy_alphas = generalised_entropy_alphas
    )

end

"""
Chop a dataset with populations and incomes
into numbins groups in a form suitable for
e.g. a Gini curve.
* col1 is cumulative population,
*    2 cumulative income/whatever,
*    3 threshold income level.
"""
function binify(
    rawdata   :: Matrix,
    numbins   :: Integer,
    weightpos :: Integer = 1,
    incomepos :: Integer = 2 ) :: Matrix
    data = make_augmented( rawdata, weightpos, incomepos )
    return binifyinternal( data, numbins )
end

"""
As above, but using any DataFrame like thing that supports the TableTraits.isiterabletable interface
"""
function binify(
    rawdata,
    numbins   :: Integer,
    weightcol :: Symbol,
    incomecol :: Symbol ) :: Matrix
    @assert TableTraits.isiterabletable( rawdata ) "data needs to implement IterableTables"
    data = make_augmented( rawdata, weightcol, incomecol )
    return binifyinternal( data, numbins )
end

## FIXME this can go horribly wrong with small unbalanced numbers of rows -
# try pop=[22,1,1,1,1,1,1,1,1,1 inc==[1,2,2,2,2,2,2,2,2,2]
function binifyinternal(
    data      :: Matrix,
    numbins   :: Integer ) :: Matrix
    T = eltype( data )

    nrows = size( data )[1]
    ncols = size( data )[2]
    out = zeros( T, numbins, 3 )
    total_population = data[ nrows, POPN_ACCUM ]
    total_income = data[ nrows, INCOME_ACCUM ]
    bin_size :: T = 1.0/numbins
    bno = 0
    thresh = bin_size
    popsharelast = 0.0
    incomelast = 0.0
    incomeaccumlast = 0.0
    incomeaccum = 0.0
    # print( "total_population $total_population total_income $total_income \n")
    for row in 1:nrows
        income = data[row,INCOME]
        incomeaccum = data[row,INCOME_ACCUM]/total_income
        popshare = data[row,POPN_ACCUM]/total_population
        if popshare ≈ thresh
            bno += 1
            out[bno,1] = popshare
            out[bno,2] = incomeaccum
            out[bno,3] = income
            # print( "row $row popshare $popshare ≈ thresh $thresh incomeaccum $incomeaccum \n")
            thresh += bin_size
        elseif( popsharelast < thresh ) && ( popshare > thresh)
            bno += 1
            # print( "row $row popsharelast $popsharelast thresh $thresh popshare $popshare\n" )
            pgap = popshare - popsharelast
            p1 = (thresh - popsharelast)
            p2 = (popshare - thresh)
            out[bno,1] = ((popshare*p2)+(popsharelast*p1))/pgap
            out[bno,2] = ((incomeaccum*p2)+(incomeaccumlast*p1))/pgap
            out[bno,3] = ((income*p2)+(incomelast*p1))/pgap
            thresh += bin_size
        end
        popsharelast = popshare
        incomelast = income
        incomeaccumlast = incomeaccum
    end
    return out
end

"""
Make a dictionary of inequality measures.
This is mainly taken from chs 5 and 6 of the World Bank book.

1. `rawdata` a matrix with cols with weights and incomes
2. `atkinson_es` inequality aversion values for the Atkinson indexes
3. `generalised_entropy_alphas`
4. `weightpos` - column with weights
5. `incomepos` - column with incomes

"""
function make_inequalityinternal(
    ;
    data                       :: Matrix,
    atkinson_es                :: Vector = DEFAULT_ATKINSON_ES,
    generalised_entropy_alphas :: Vector = DEFAULT_ENTROPIES ) :: InequalityMeasures
    nrows = size( data )[1]
    ncols = size( data )[2]
    @assert ncols == 5 "data should have 5 cols but has $ncols"
    T = eltype( generalised_entropy_alphas )
    z = zero(T)
    nats = size( atkinson_es )[1]
    neps = size( generalised_entropy_alphas )[1]
    iq =InequalityMeasures(
        atkinson_es,
        [i == 1.0 ? 1.0 : 0.0 for i in atkinson_es],
        generalised_entropy_alphas,
        zeros(T,neps),
        z,
        z,
        z,
        z,
        z,
        z,
        z,
        z,
        z,
        zeros( T, 10, 3 ),
        0 )

    # initialise atkinsons; 1 for e = 1 0 otherwise

    iq.gini = make_gini( data )
    positive_inc_popn = calc_positive_inc_popn( data )
    total_income = data[nrows,INCOME_ACCUM]
    total_population = data[nrows,POPN_ACCUM]
    y_bar = total_income/total_population
    iq.total_income = total_income
    iq.total_population = total_population
    iq.average_income = y_bar
    popsharelast = 0.0
    incomelast = 0.0
    for row in 1:nrows
        income = data[row,INCOME]
        weight = data[row,WEIGHT]
        y_yb  :: T = income/y_bar
        yb_y  :: T = y_bar/income
        iq.hoover += weight*abs( income - y_bar )
        for i in 1:neps
            alpha :: T = iq.generalised_entropy_alphas[i]
            iq.generalised_entropy[i] += weight*(y_yb^alpha)
        end # entropies
        # atkinson kinda sorta needs
        # to be over +ives only since otherwise atk(1) = 1 always
        # since the inner sum goes to 0
        if income > 0.0
            ln_y_yb :: T = log( y_yb )
            ln_yb_y :: T = log( yb_y )
            iq.theil_l += weight*ln_yb_y
            iq.theil_t += weight*y_yb*ln_y_yb
            for i in 1:nats
                es :: T = iq.atkinson_es[i]
                if es != 1.0
                    iq.atkinson[i] += weight*(y_yb^(1.0-es))
                else
                    iq.atkinson[i] *= (income)^(weight/positive_inc_popn)
                end # e = 1 case
            end # atkinsons
        else
            iq.negative_or_zero_income_count += weight
        end # positive income
    end # main loop

    @assert (iq.negative_or_zero_income_count + positive_inc_popn) ≈ iq.total_population "nzc $(iq.negative_or_zero_income_count) pip $positive_inc_popn tp $(iq.total_population)"
    deciles = binify( data, 10 )
    iq.median = deciles[5,3]
    # top 10/bottom 40
    if deciles[4,2] > 0
        iq.palma = (1.0-deciles[9,2])/deciles[4,2]
    end
    iq.deciles = deciles
    iq.hoover /= 2.0*total_income
    for i in 1:neps
        alpha :: T = iq.generalised_entropy_alphas[i]
        aq :: T = (1.0/(alpha*(alpha-1.0)))
        iq.generalised_entropy[i] =
            aq*((iq.generalised_entropy[i]/iq.total_population)-1.0)
    end # entropies
    for i in 1:nats
        es :: T = iq.atkinson_es[i]
        if es != 1.0
            iq.atkinson[i] = 1.0 - ( iq.atkinson[i]/positive_inc_popn)^(1.0/(1.0-es))
        else
            iq.atkinson[i] = 1.0 - ( iq.atkinson[i]/y_bar )
        end # e = 1
    end
    iq.theil_l /= positive_inc_popn
    iq.theil_t /= positive_inc_popn
    return iq
end # make_inequalityinternal

end


text:
export null, isnull, all_ones, UndefValue, PointerNull

null(typ::LLVMType) = Value(API.LLVMConstNull(typ))

all_ones(typ::LLVMType) = Value(API.LLVMConstAllOnes(typ))

isnull(val::Value) = convert(Core.Bool, API.LLVMIsNull(val))

abstract type Constant <: User end

# forward declarations
@checked mutable struct Module
    ref::API.LLVMModuleRef
end
abstract type Instruction <: User end

@checked struct PointerNull <: Constant
    ref::API.LLVMValueRef
end
identify(::Type{Value}, ::Val{API.LLVMConstantPointerNullValueKind}) = PointerNull

PointerNull(typ::PointerType) = PointerNull(API.LLVMConstPointerNull(typ))

@checked struct UndefValue <: Constant
    ref::API.LLVMValueRef
end
identify(::Type{Value}, ::Val{API.LLVMUndefValueValueKind}) = UndefValue

UndefValue(typ::LLVMType) = UndefValue(API.LLVMGetUndef(typ))

## scalar

export ConstantInt, ConstantFP

@checked struct ConstantInt <: Constant
    ref::API.LLVMValueRef
end
identify(::Type{Value}, ::Val{API.LLVMConstantIntValueKind}) = ConstantInt

# NOTE: fixed set for dispatch, also because we can't rely on sizeof(T)==width(T)
const WideInteger = Union{Int64, UInt64}
ConstantInt(typ::IntegerType, val::WideInteger, signed=false) =
    ConstantInt(API.LLVMConstInt(typ, reinterpret(Culonglong, val),
                convert(Bool, signed)))
const SmallInteger = Union{Int8, Int16, Int32, UInt8, UInt16, UInt32}
ConstantInt(typ::IntegerType, val::SmallInteger, signed=false) =
    ConstantInt(typ, convert(Int64, val), signed)

function ConstantInt(typ::IntegerType, val::Integer, signed=false)
    valbits = ceil(Int, log2(abs(val))) + 1
    numwords = ceil(Int, valbits / 64)
    words = Vector{Culonglong}(undef, numwords)
    for i in 1:numwords
        words[i] = (val >> 64(i-1)) % Culonglong
    end
    return ConstantInt(API.LLVMConstIntOfArbitraryPrecision(typ, numwords, words))
end

# NOTE: fixed set where sizeof(T) does match the numerical width
const SizeableInteger = Union{Int8, Int16, Int32, Int64, Int128, UInt8, UInt16, UInt32, UInt64, UInt128}
function ConstantInt(val::T, ctx::Context) where T<:SizeableInteger
    typ = IntType(sizeof(T)*8, ctx)
    return ConstantInt(typ, val, T<:Signed)
end

Base.convert(::Type{T}, val::ConstantInt) where {T<:Unsigned} =
    convert(T, API.LLVMConstIntGetZExtValue(val))

Base.convert(::Type{T}, val::ConstantInt) where {T<:Signed} =
    convert(T, API.LLVMConstIntGetSExtValue(val))

@checked struct ConstantFP <: Constant
    ref::API.LLVMValueRef
end
identify(::Type{Value}, ::Val{API.LLVMConstantFPValueKind}) = ConstantFP

ConstantFP(typ::FloatingPointType, val::Real) =
    ConstantFP(API.LLVMConstReal(typ, Cdouble(val)))

ConstantFP(val::Float16, ctx::Context) =
    ConstantFP(HalfType(ctx), val)
ConstantFP(val::Float32, ctx::Context) =
    ConstantFP(FloatType(ctx), val)
ConstantFP(val::Float64, ctx::Context) =
    ConstantFP(DoubleType(ctx), val)

Base.convert(::Type{T}, val::ConstantFP) where {T<:AbstractFloat} =
    convert(T, API.LLVMConstRealGetDouble(val, Ref{API.LLVMBool}()))

## aggregate

export ConstantAggregateZero

@checked struct ConstantAggregateZero <: Constant
    ref::API.LLVMValueRef
end
identify(::Type{Value}, ::Val{API.LLVMConstantAggregateZeroValueKind}) = ConstantAggregateZero

# there currently seems to be no function in the LLVM-C interface which returns a
# ConstantAggregateZero value directly, but values can occur through calls to LLVMConstNull

## constant expressions

export ConstantExpr, ConstantAggregate, ConstantArray, ConstantStruct, ConstantVector, InlineAsm

@checked struct ConstantExpr <: Constant
    ref::API.LLVMValueRef
end
identify(::Type{Value}, ::Val{API.LLVMConstantExprValueKind}) = ConstantExpr

abstract type ConstantAggregate <: Constant end

@checked struct ConstantArray <: ConstantAggregate
    ref::API.LLVMValueRef
end
identify(::Type{Value}, ::Val{API.LLVMConstantArrayValueKind}) = ConstantArray
identify(::Type{Value}, ::Val{API.LLVMConstantDataArrayValueKind}) = ConstantArray

ConstantArray(typ::LLVMType, data::Vector{T}) where {T<:Constant} =
    ConstantArray(API.LLVMConstArray(typ, data, length(data)))
ConstantArray(typ::IntegerType, data::Vector{T}) where {T<:Integer} =
    ConstantArray(typ, map(x->ConstantInt(convert(T,x),context(typ)), data))
ConstantArray(typ::FloatingPointType, data::Vector{T}) where {T<:AbstractFloat} =
    ConstantArray(typ, map(x->ConstantFP(convert(T,x),context(typ)), data))

Base.getindex(ca::ConstantArray, idx::Integer) =
    API.LLVMGetElementAsConstant(ca, idx-1)
Base.length(ca::ConstantArray) = length(llvmtype(ca))
Base.eltype(ca::ConstantArray) = eltype(llvmtype(ca))
Base.convert(::Type{Array{T,1}}, ca::ConstantArray) where {T<:Integer} =
    [convert(T,ConstantInt(ca[i])) for i in 1:length(ca)]
Base.convert(::Type{Array{T,1}}, ca::ConstantArray) where {T<:AbstractFloat} =
    [convert(T,ConstantFP(ca[i])) for i in 1:length(ca)]

@checked struct ConstantStruct <: ConstantAggregate
    ref::API.LLVMValueRef
end
identify(::Type{Value}, ::Val{API.LLVMConstantStructValueKind}) = ConstantStruct

@checked struct ConstantVector <: ConstantAggregate
    ref::API.LLVMValueRef
end
identify(::Type{Value}, ::Val{API.LLVMConstantVectorValueKind}) = ConstantVector

@checked struct InlineAsm <: Constant
    ref::API.LLVMValueRef
end
identify(::Type{Value}, ::Val{API.LLVMInlineAsmValueKind}) = InlineAsm

InlineAsm(typ::FunctionType, asm::String, constraints::String,
          side_effects::Core.Bool, align_stack::Core.Bool=false) =
    InlineAsm(API.LLVMConstInlineAsm(typ, asm, constraints,
                                     convert(Bool, side_effects),
                                     convert(Bool, align_stack)))

## global values

abstract type GlobalValue <: Constant end

export GlobalValue,
       isdeclaration,
       linkage, linkage!,
       section, section!,
       visibility, visibility!,
       dllstorage, dllstorage!,
       unnamed_addr, unnamed_addr!,
       alignment, alignment!

parent(val::GlobalValue) = Module(API.LLVMGetGlobalParent(val))

isdeclaration(val::GlobalValue) = convert(Core.Bool, API.LLVMIsDeclaration(val))

linkage(val::GlobalValue) = API.LLVMGetLinkage(val)
linkage!(val::GlobalValue, linkage::API.LLVMLinkage) =
    API.LLVMSetLinkage(val, linkage)

function section(val::GlobalValue)
  #=
  The following started to fail on LLVM 4.0:
    Context() do ctx
      LLVM.Module("SomeModule", ctx) do mod
        st = LLVM.StructType("SomeType", ctx)
        ft = LLVM.FunctionType(st, [st])
        fn = LLVM.Function(mod, "SomeFunction", ft)
        section(fn) == ""
      end
      end
  =#
  section_ptr = API.LLVMGetSection(val)
  return section_ptr != C_NULL ? unsafe_string(section_ptr) : ""
end
section!(val::GlobalValue, sec::String) = API.LLVMSetSection(val, sec)

visibility(val::GlobalValue) = API.LLVMGetVisibility(val)
visibility!(val::GlobalValue, viz::API.LLVMVisibility) =
    API.LLVMSetVisibility(val, viz)

dllstorage(val::GlobalValue) = API.LLVMGetDLLStorageClass(val)
dllstorage!(val::GlobalValue, storage::API.LLVMDLLStorageClass) =
    API.LLVMSetDLLStorageClass(val, storage)

unnamed_addr(val::GlobalValue) = convert(Core.Bool, API.LLVMHasUnnamedAddr(val))
unnamed_addr!(val::GlobalValue, flag::Core.Bool) =
    API.LLVMSetUnnamedAddr(val, convert(Bool, flag))

const AlignedValue = Union{GlobalValue,Instruction}   # load, store, alloca
alignment(val::AlignedValue) = API.LLVMGetAlignment(val)
alignment!(val::AlignedValue, bytes::Integer) = API.LLVMSetAlignment(val, bytes)

## global variables

abstract type GlobalObject <: GlobalValue end

export GlobalVariable, unsafe_delete!,
       initializer, initializer!,
       isthreadlocal, threadlocal!,
       threadlocalmode, threadlocalmode!,
       isconstant, constant!,
       isextinit, extinit!

@checked struct GlobalVariable <: GlobalObject
    ref::API.LLVMValueRef
end
identify(::Type{Value}, ::Val{API.LLVMGlobalVariableValueKind}) = GlobalVariable

GlobalVariable(mod::Module, typ::LLVMType, name::String) =
    GlobalVariable(API.LLVMAddGlobal(mod, typ, name))

GlobalVariable(mod::Module, typ::LLVMType, name::String, addrspace::Integer) =
    GlobalVariable(API.LLVMAddGlobalInAddressSpace(mod, typ,
                                                   name, addrspace))

unsafe_delete!(::Module, gv::GlobalVariable) = API.LLVMDeleteGlobal(gv)

function initializer(gv::GlobalVariable)
    init = API.LLVMGetInitializer(gv)
    init == C_NULL ? nothing : Value(init)
end
initializer!(gv::GlobalVariable, val::Constant) =
  API.LLVMSetInitializer(gv, val)
initializer!(gv::GlobalVariable, ::Nothing) =
  API.LLVMSetInitializer(gv, C_NULL)

isthreadlocal(gv::GlobalVariable) = convert(Core.Bool, API.LLVMIsThreadLocal(gv))
threadlocal!(gv::GlobalVariable, bool) =
  API.LLVMSetThreadLocal(gv, convert(Bool, bool))

isconstant(gv::GlobalVariable) = convert(Core.Bool, API.LLVMIsGlobalConstant(gv))
constant!(gv::GlobalVariable, bool) =
  API.LLVMSetGlobalConstant(gv, convert(Bool, bool))

threadlocalmode(gv::GlobalVariable) = API.LLVMGetThreadLocalMode(gv)
threadlocalmode!(gv::GlobalVariable, mode) =
  API.LLVMSetThreadLocalMode(gv, mode)

isextinit(gv::GlobalVariable) =
  convert(Core.Bool, API.LLVMIsExternallyInitialized(gv))
extinit!(gv::GlobalVariable, bool) =
  API.LLVMSetExternallyInitialized(gv, convert(Bool, bool))

ext: jl


text:
shared_parameters = """    Xtrain::M where M <: AbstractMatrix{<:AbstractFloat}
Matrix of training points.

    Xtest::M where M <: AbstractMatrix{<:AbstractFloat}
Matrix of test points.

    k::Int
Nuumber of nearest neighbors to find in the training points.

    batch_size::Int
Batch size to split the training and/or test data.

    metric::PreMetric
(Pre-) Metric to use for distance computation.

    convert::Function
Conversion function to convert the input data (`Xtrain` and `Xtest`) to an appropriate format, e.g. `CUDA.cu`."""

"""
    knn_hybrid_batch_test(Xtrain, Xtest, k, batch_size; metric, convert)

A parallel KNN implementation to determine the nearest neighbors of a number of points, where each column in  `Xtrain`
and `Xtest` corresponds to a point. This implementation batches points of `Xtest`, and calculates the nearest neighbors
to all points in `Xtrain`.

Parameters
----------
$shared_parameters
"""
function knn_hybrid_batch_test(Xtrain::M, Xtest::M, k::Int, batch_size::Int = guess_batchsize(Xtest, k);
    metric::PreMetric = Euclidean(), convert::Function = identity) where {M<:input_type}
    n_train = size(Xtrain, 2)
    n_test = size(Xtest, 2)
    @assert batch_size <= n_test "Batch size must be smaller or equal the number of test points."

    batches = Iterators.partition(1:n_test, batch_size)
    idxs = [Vector{Int}(undef, k) for _ in 1:n_test]
    dists = [Vector{eltype(Xtrain)}(undef, k) for _ in 1:n_test]
    workspace_idxs = [collect(1:n_train) for _ in 1:batch_size]

    # Convert train points according to conversion function, e.g. `cu`
    Xtrain = convert(Xtrain)

    # TODO: Can we pin / page-lock this matrix generically for AMD/CUDA?
    mcpu = Matrix{eltype(Xtrain)}(undef, n_train, batch_size)
    mgpu = similar(Xtrain, n_train, batch_size)
    sort_finished = Task(() -> nothing) |> schedule

    # Convert the metric to a kernel
    distance_kernel = metric_to_kernel(metric)

    for (i, batch) in enumerate(batches)
        n_batch = length(batch)
        mgpu[:, 1:n_batch] .= distance_kernel(convert(Xtest[:, batch]), Xtrain)
        wait(sort_finished)
        copyto!(mcpu, mgpu)
        sort_finished = async_cpu_colsort!(idxs, dists, workspace_idxs, mcpu, batch_size, n_batch, i, k)
    end
    wait(sort_finished)
    return idxs, dists
end

"""
    knn_hybrid_batch_all(Xtrain, Xtest, k, batch_size; metric, convert)

A mixed GPU/CPU implementation to determine the nearest neighbors of a number of points, where each column in  `Xtrain`
and `Xtest` corresponds to a point. This implementation batches points of `Xtrain` and `Xtest` simultaneously, saves the
temporary nearest neighbors and reduces the temporary neighbors after all other points in `Xtrain` have been observed.

Parameters
----------
$shared_parameters
"""
function knn_hybrid_batch_all(Xtrain::M, Xtest::M, k::Int, batch_size::Int = guess_batchsize(Xtest, k);
    metric::PreMetric = Euclidean(), convert::Function = identity) where {M<:input_type}
    n_train = size(Xtrain, 2)
    n_test = size(Xtest, 2)
    @assert batch_size >= k "Batch size must be larger or equal to k."
    @assert batch_size <= n_train "Batch size must be smaller or equal the number of training points."

    # generate the partitions
    train_batches = Iterators.partition(1:n_train, batch_size)
    test_batches = Iterators.partition(1:n_test, batch_size)

    # get sizes to preallocate the temporary store of all possible nearest neighbors
    temp_length = map(train_batches) do batch
        len = length(batch)
        len > k ? k : len
    end |> sum

    # pre-allocate the final outputs
    idxs = [Vector{Int}(undef, k) for _ in 1:n_test]
    dists = [Vector{eltype(Xtrain)}(undef, k) for _ in 1:n_test]

    # pre-allocate the temporary store of the distance matrix
    mcpu = Matrix{eltype(Xtrain)}(undef, batch_size, batch_size)
    mgpu = similar(convert(Xtrain[1:1]), batch_size, batch_size)

    # pre-allocate the temporary store of the indices and distances
    temp_idxs = Matrix{Int}(undef, temp_length, n_test)
    temp_dists = Matrix{eltype(Xtrain)}(undef, temp_length, n_test)

    # pre-allocate the temporary sorting indices (there are never more than batch_size indices needed)
    outer_workspace = [collect(1:temp_length) for _ in 1:n_test]
    inner_workspace = [collect(1:batch_size) for _ in 1:batch_size]

    # this will hold a task that determinies if sorting is already done
    sort_finished = Task(() -> nothing) |> schedule
    reduction_finished = Task(() -> nothing) |> schedule

    # Convert the metric to a kernel
    distance_kernel = metric_to_kernel(metric)

    for batch_test in test_batches
        n_batch_test = length(batch_test)
        for (idx_train, batch_train) in enumerate(train_batches)
            n_train_batch = length(batch_train)

            # calculate the distance matrix, it's important that we restrict the pre-allocated matrix to the
            # possible points in the batch (which might be less then the full `batchsize x batchsize` matrix)
            mgpu[1:n_train_batch, 1:n_batch_test] .=
                distance_kernel(convert(Xtest[:, batch_test]), convert(Xtrain[:, batch_train]))

            # we calculate the distances in parallel, but we need to wait for the sorting before we can copy
            wait(sort_finished)
            copyto!(mcpu, mgpu)

            # prepare ranges of indices for the currently calculated neighbors, batches might be smaller than `k`
            nbatch_or_k = n_train_batch > k ? k : n_train_batch
            range_lower = (idx_train - 1) * k + 1
            range_upper = range_lower + nbatch_or_k - 1

            # we are working with a fixed-size square matrix of `batch_size x batch_size`, but only
            # some values might be used for the last batches (that are smaller than batch_size)
            sort_finished = async_cpu_batchsort!(temp_idxs, temp_dists, inner_workspace, mcpu,
                range_lower:range_upper, batch_train, batch_test, n_train_batch, n_batch_test, nbatch_or_k)
        end
        # Note: It's important to wait for the tasks in the corresponding loops the are started from, otherwise (only
        # waiting after the loop) it is not guaranteed that all tasks finished, but only the last one.
        wait(sort_finished)
        wait(reduction_finished)
        reduction_finished = async_cpu_reduction!(idxs, dists, temp_idxs, temp_dists, outer_workspace, batch_test, k)
    end
    wait(reduction_finished)
    return idxs, dists
end

"""
    knn_full(Xtrain, Xtest, k, batch_size, distance_kernel)

A trivial implementation to determine the nearest neighbors of a number of points, where each column in  `Xtrain` and
`Xtest` corresponds to a point. This implementation calculates all distances at once, sorting afterwards. Note: because
there is no `partialsortperm` GPU kernel available, this implementation inefficiently uses `sortperm` to sort the
resulting distance matrix.

Parameters
----------
    Xtrain::M where M <: AbstractGPUArray{<:AbstractFloat}
Matrix of training points.

    Xtest::M where M <: AbstractGPUArray{<:AbstractFloat}
Matrix of test points.

    k::Int
Nuumber of nearest neighbors to find in the training points.

    metric::PreMetric
Metric to use for distance computation.

    convert::Function
Conversion function to convert the input data (`Xtrain` and `Xtest`) to an appropriate format, e.g. `CUDA.cu`. In this
implementation, both `Xtrain` and `Xtest` are converted beforehand.
"""
function knn_full(Xtrain::M, Xtest::M, k::Int;
    metric::PreMetric = Euclidean(), convert::Function = identity) where {M<:input_type}
    n_test = size(Xtest, 2)
    idxs = [Vector{Int}(undef, k) for _ in 1:n_test]
    dists = [Vector{eltype(convert(Xtrain[1:1]))}(undef, k) for _ in 1:n_test]
    distance_kernel = metric_to_kernel(metric)
    distance_matrix = distance_kernel(convert(Xtest), convert(Xtrain))
    colsort!(idxs, dists, distance_matrix, k)
    return idxs, dists
end

function colsort!(idxs, dists, distance_matrix, k)
    for (i, col) in enumerate(eachcol(distance_matrix))
        # TODO: currently, there is no partialsortperm kernel available
        sortres = view(sortperm(col), 1:k)
        copyto!(idxs[i], sortres)
        copyto!(dists[i], col[sortres])
    end
end

function async_cpu_colsort!(idxs, dists, workspace, distance_matrix, batch_size, n_batch, current_batch, k)
    return @async Threads.@threads for col in 1:n_batch
        idx = ((current_batch - 1) * batch_size) + col
        sorted_idxs = partialsortperm!(workspace[col],
            view(distance_matrix, :, col), 1:k, rev = false, initialized = true)
        copyto!(idxs[idx], sorted_idxs)
        copyto!(dists[idx], distance_matrix[sorted_idxs, col])
        nothing # JuliaLang/julia#40626
    end
end

function async_cpu_batchsort!(idxs, dists, workspace, distance_matrix, idx_range, batch_train, batch_test,
    n_batch_train, n_batch_test, n_or_k_j)
    return @async Threads.@threads for col in 1:n_batch_test
        # TODO: this sort could probably be sped up, we currently have to rely on resizing and initializing
        # because not all batches must have the same size, if all have the same size, however we can assume
        # `initalized = true` and don't need to resize
        idx = partialsortperm!(resize!(workspace[col], n_batch_train),
            view(distance_matrix, 1:n_batch_train, col), 1:n_or_k_j, rev = false, initialized = false)
        @inbounds idxs[idx_range, batch_test[col]] = view(batch_train, idx)
        @inbounds dists[idx_range, batch_test[col]] = view(distance_matrix, idx, col)
        nothing # JuliaLang/julia#40626
    end
end

function async_cpu_reduction!(idxs, dists, temp_idxs, temp_dists, workspace, batch_test, k)
    return @async Threads.@threads for i in batch_test
        idx = partialsortperm!(workspace[i], view(temp_dists, :, i), 1:k, rev = false, initialized = true)
        copyto!(idxs[i], temp_idxs[idx, i])
        copyto!(dists[i], temp_dists[idx, i])
        nothing # JuliaLang/julia#40626
    end
end


text:
module Utils

using CUDA

"""
    unsafe_free_all!(namespace::Module)

Call `CUDA.unsafe_free!` on all `CuArray`s in `namespace`.
"""
function unsafe_free_all!(namespace::Module)
    for n in names(namespace; all = true)
        x = try
            getfield(namespace, n)
        catch err
            err isa UndefVarError || rethrow()
            continue
        end
        if x isa CuArray
            @debug "Free: `$namespace.$n`"
            CUDA.unsafe_free!(x)
        end
    end
end

end  # module


text:
# Streaming XML Reader
# ====================

struct _TextReader
    # type tag
end

"""
A streaming XML reader type.
"""
mutable struct StreamReader
    ptr::Ptr{_TextReader}
    input::Union{IO,Cvoid}

    function StreamReader(ptr::Ptr{_TextReader}, input=nothing)
        @assert ptr != C_NULL
        return new(ptr, input)
    end
end

@inline function Base.getproperty(reader::StreamReader, name::Symbol)
    name == :type      ? nodetype(reader)                                         :
    name == :depth     ? nodedepth(reader)                                        :
    name == :name      ? (hasnodename(reader)    ? nodename(reader)    : nothing) :
    name == :content   ? (hasnodecontent(reader) ? nodecontent(reader) : nothing) :
    name == :namespace ? namespace(reader)                                        :
    Core.getfield(reader, name)
end

function Base.show(io::IO, reader::StreamReader)
    prefix = isdefined(Main, :StreamReader) ? "StreamReader" : "EzXML.StreamReader"
    @printf(io, "%s(<%s@%p>)", prefix, repr(nodetype(reader)), reader.ptr)
end

# Reader type (enum xmlReaderTypes).
if sizeof(Cint) == 2
    primitive type ReaderType <: Integer 16 end
elseif sizeof(Cint) == 4
    primitive type ReaderType <: Integer 32 end
elseif sizeof(Cint) == 8
    primitive type ReaderType <: Integer 64 end
else
    @assert false "invalid Cint size"
end

function ReaderType(x::Integer)
    return convert(ReaderType, x)
end

function Base.convert(::Type{ReaderType}, x::Integer)
    return reinterpret(ReaderType, convert(Cint, x))
end

function Base.convert(::Type{T}, x::ReaderType) where {T<:Integer}
    return convert(T, reinterpret(Cint, x))
end

Base.hash(x::ReaderType, h::UInt) = hash(convert(Cint,x), h)

function Base.convert(::Type{ReaderType}, x::ReaderType)
    return x
end

function Base.promote_rule(::Type{ReaderType}, ::Type{T}) where {T<:Union{Cint,Int}}
    return T
end

const READER_NONE                   = ReaderType( 0)
const READER_ELEMENT                = ReaderType( 1)
const READER_ATTRIBUTE              = ReaderType( 2)
const READER_TEXT                   = ReaderType( 3)
const READER_CDATA                  = ReaderType( 4)
const READER_ENTITY_REFERENCE       = ReaderType( 5)
const READER_ENTITY                 = ReaderType( 6)
const READER_PROCESSING_INSTRUCTION = ReaderType( 7)
const READER_COMMENT                = ReaderType( 8)
const READER_DOCUMENT               = ReaderType( 9)
const READER_DOCUMENT_TYPE          = ReaderType(10)
const READER_DOCUMENT_FRAGMENT      = ReaderType(11)
const READER_NOTATION               = ReaderType(12)
const READER_WHITESPACE             = ReaderType(13)
const READER_SIGNIFICANT_WHITESPACE = ReaderType(14)
const READER_END_ELEMENT            = ReaderType(15)
const READER_END_ENTITY             = ReaderType(16)
const READER_XML_DECLARATION        = ReaderType(17)

function Base.show(io::IO, x::ReaderType)
    if x == READER_NONE
        print(io, "READER_NONE")
    elseif x == READER_ELEMENT
        print(io, "READER_ELEMENT")
    elseif x == READER_ATTRIBUTE
        print(io, "READER_ATTRIBUTE")
    elseif x == READER_TEXT
        print(io, "READER_TEXT")
    elseif x == READER_CDATA
        print(io, "READER_CDATA")
    elseif x == READER_ENTITY_REFERENCE
        print(io, "READER_ENTITY_REFERENCE")
    elseif x == READER_ENTITY
        print(io, "READER_ENTITY")
    elseif x == READER_PROCESSING_INSTRUCTION
        print(io, "READER_PROCESSING_INSTRUCTION")
    elseif x == READER_COMMENT
        print(io, "READER_COMMENT")
    elseif x == READER_DOCUMENT
        print(io, "READER_DOCUMENT")
    elseif x == READER_DOCUMENT_TYPE
        print(io, "READER_DOCUMENT_TYPE")
    elseif x == READER_DOCUMENT_FRAGMENT
        print(io, "READER_DOCUMENT_FRAGMENT")
    elseif x == READER_NOTATION
        print(io, "READER_NOTATION")
    elseif x == READER_WHITESPACE
        print(io, "READER_WHITESPACE")
    elseif x == READER_SIGNIFICANT_WHITESPACE
        print(io, "READER_SIGNIFICANT_WHITESPACE")
    elseif x == READER_END_ELEMENT
        print(io, "READER_END_ELEMENT")
    elseif x == READER_END_ENTITY
        print(io, "READER_END_ENTITY")
    elseif x == READER_XML_DECLARATION
        print(io, "READER_XML_DECLARATION")
    else
        @assert false "unknown reader type"
    end
end

function Base.print(io::IO, x::ReaderType)
    print(io, convert(Cint, x))
end

function Base.string(x::ReaderType)
    return sprint(print, x)
end

function StreamReader(input::IO)
    readcb = make_read_callback()
    closecb = C_NULL
    context = input
    uri = C_NULL
    encoding = C_NULL
    options = 0
    reader_ptr = @check ccall(
        (:xmlReaderForIO, libxml2),
        Ptr{_TextReader},
        (Ptr{Cvoid}, Ptr{Cvoid}, Ref{IO}, Cstring, Cstring, Cint),
        readcb, closecb, context, uri, encoding, options) != C_NULL
    return StreamReader(reader_ptr, input)
end

function Base.open(::Type{StreamReader}, filename::AbstractString)
    encoding = C_NULL
    options = 0
    reader_ptr = @check ccall(
        (:xmlReaderForFile, libxml2),
        Ptr{_TextReader},
        (Cstring, Cstring, Cint),
        filename, encoding, options) != C_NULL
    return StreamReader(reader_ptr)
end

function Base.close(reader::StreamReader)
    ccall(
        (:xmlFreeTextReader, libxml2),
        Cvoid,
        (Ptr{Cvoid},),
        reader.ptr)
    reader.ptr = C_NULL
    if reader.input isa IO
        close(reader.input)
    end
    return nothing
end

function Base.eltype(::Type{StreamReader})
    return ReaderType
end

function Base.IteratorSize(::Type{StreamReader})
    return Base.SizeUnknown()
end

function Base.iterate(reader::StreamReader, state=nothing)
    if read_node(reader)
        nothing
    else
        nodetype(reader), nothing
    end
end

# Read a next node and return `true` iff finished.
function read_node(reader)
    ret = @check ccall(
        (:xmlTextReaderRead, libxml2),
        Cint,
        (Ptr{Cvoid},),
        getfield(reader, :ptr)) ≥ 0
    return ret == 0
end

"""
    nodedepth(reader::StreamReader)

Return the depth of the current node of `reader`.
"""
function nodedepth(reader::StreamReader)
    ret = ccall(
        (:xmlTextReaderDepth, libxml2),
        Cint,
        (Ptr{Cvoid},),
        getfield(reader, :ptr))
    return Int(ret)
end

"""
    nodetype(reader::StreamReader)

Return the type of the current node of `reader`.
"""
function nodetype(reader::StreamReader)
    typ = ccall(
        (:xmlTextReaderNodeType, libxml2),
        Cint,
        (Ptr{Cvoid},),
        getfield(reader, :ptr))
    return convert(ReaderType, typ)
end

"""
    hasnodename(reader::StreamReader)

Return if the current node of `reader` has a node name.
"""
function hasnodename(reader::StreamReader)
    return ccall(
        (:xmlTextReaderConstName, libxml2),
        Cstring,
        (Ptr{Cvoid},),
        getfield(reader, :ptr)) != C_NULL
end

"""
    nodename(reader::StreamReader)

Return the name of the current node of `reader`.
"""
function nodename(reader::StreamReader)
    name_ptr = ccall(
        (:xmlTextReaderConstName, libxml2),
        Cstring,
        (Ptr{Cvoid},),
        getfield(reader, :ptr))
    if name_ptr == C_NULL
        throw(ArgumentError("no node name"))
    end
    return unsafe_string(name_ptr)
end

"""
    hasnodecontent(reader::StreamReader)

Return if the current node of `reader` has content.
"""
function hasnodecontent(reader::StreamReader)
    if nodetype(reader) == READER_ATTRIBUTE
        return false
    end
    # TODO: this allocates memory; any way to avoid it?
    ptr = ccall(
        (:xmlTextReaderReadString, libxml2),
        Cstring,
        (Ptr{Cvoid},),
        getfield(reader, :ptr))
    if ptr == C_NULL
        return false
    else
        Libc.free(ptr)
        return true
    end
end

"""
    nodecontent(reader::StreamReader)

Return the content of the current node of `reader`.
"""
function nodecontent(reader::StreamReader)
    if nodetype(reader) == READER_ATTRIBUTE
        throw(ArgumentError("no content"))
    end
    content_ptr = ccall(
        (:xmlTextReaderReadString, libxml2),
        Cstring,
        (Ptr{Cvoid},),
        getfield(reader, :ptr))
    if content_ptr == C_NULL
        throw(ArgumentError("no content"))
    end
    content = unsafe_string(content_ptr)
    Libc.free(content_ptr)
    return content
end

"""
    hasnodevalue(reader::StreamReader)

Return if the current node of `reader` has a value.
"""
function hasnodevalue(reader::StreamReader)
    r = ccall(
       (:xmlTextReaderHasValue, libxml2),
       Cint,
       (Ptr{Cvoid},),
       getfield(reader, :ptr))
    return r == 1
end

"""
    nodevalue(reader::StreamReader)

Return the value of the current node of `reader`.

This can be different from `nodecontent`.
"""
function nodevalue(reader::StreamReader)
    value_ptr = ccall(
        (:xmlTextReaderConstValue, libxml2),
        Cstring,
        (Ptr{Cvoid},),
        getfield(reader, :ptr))
    if value_ptr == C_NULL
        throw(ArgumentError("no node value"))
    end
    return unsafe_string(value_ptr)
end

"""
    hasnodeattributes(reader::StreamReader)

Return if the current node of 'reader' has attributes
"""
function hasnodeattributes(reader::StreamReader)
    r = ccall(
       (:xmlTextReaderHasAttributes, libxml2),
       Cint,
       (Ptr{Cvoid},),
       getfield(reader, :ptr))
    @assert r ≥ 0 "XML Error Detected"
    return r == 1
end

struct AttributeReader
    reader::StreamReader

    function AttributeReader(reader::StreamReader)
        if nodetype(reader) != READER_ELEMENT
            throw(ArgumentError("Reader not an Element Node"))
        end
        return new(reader)
    end
end

function Base.eltype(::Type{AttributeReader})
    return StreamReader
end

function Base.IteratorSize(::Type{AttributeReader})
    return Base.SizeUnknown()
end

function Base.iterate(attrs::AttributeReader, state=nothing)
    r = ccall(
        (:xmlTextReaderMoveToNextAttribute, libxml2),
        Cint,
        (Ptr{Cvoid},),
        getfield(attrs.reader, :ptr))
    if r == 1
        return attrs.reader, nothing
    end
    if nodetype(attrs.reader) == READER_ATTRIBUTE
        r = ccall(
            (:xmlTextReaderMoveToElement, libxml2),
            Cint,
            (Ptr{Cvoid},),
            getfield(attrs.reader, :ptr))
        @assert r == 1
    end
    return nothing
end

"""
    eachattribute(reader::StreamReader)

Return an `AttributeReader` object for the current node of `reader`
"""
eachattribute(reader::StreamReader) = AttributeReader(reader)

"""
    countattributes(reader::StreamReader)

Count the number of attributes in the current node of `reader`.
"""
function countattributes(reader::StreamReader)
    r = ccall(
        (:xmlTextReaderAttributeCount, libxml2),
        Cint,
        (Ptr{Cvoid},),
        getfield(reader, :ptr))
    return Int(r)
end

"""
    nodeattributes(reader::StreamReader)

Return a dictionary of the attributes in the current node of `reader`.
"""
function nodeattributes(reader::StreamReader)
    attrs = Dict{String,String}()
    for attr in eachattribute(reader)
        attrs[nodename(attr)] = nodevalue(attr)
    end
    return attrs
end

function attribute_ptr(reader::StreamReader, name::AbstractString)
    value_ptr = ccall(
        (:xmlTextReaderGetAttribute, libxml2),
        Cstring,
        (Ptr{Cvoid}, Cstring),
        getfield(reader, :ptr), name)
end

function attribute_ptr(reader::StreamReader, no::Integer)
    value_ptr = ccall(
        (:xmlTextReaderGetAttributeNo, libxml2),
        Cstring,
        (Ptr{Cvoid}, Cint),
        getfield(reader, :ptr), Cint(no - 1))
end

"""
    haskey(reader::StreamReader, key::Union{Integer,AbstractString})

Check if current node of `reader` has attribute `key`.
"""
function Base.haskey(reader::StreamReader, key::Union{Integer,AbstractString})
    value_ptr = attribute_ptr(reader, key)
    ret = value_ptr != C_NULL
    Libc.free(value_ptr)
    return ret
end

"""
    getindex(reader::StreamReader, key::Union{Integer,AbstractString})

Get attribute `key` at current node of `reader`.
"""
function Base.getindex(reader::StreamReader, key::Union{Integer,AbstractString})
    value_ptr = attribute_ptr(reader, key)
    if value_ptr == C_NULL
        throw(KeyError(key))
    end
    value = unsafe_string(value_ptr)
    Libc.free(value_ptr)
    return value
end

"""
    namespace(reader::StreamReader)

Return the namespace of the current node of `reader`.
"""
function namespace(reader::StreamReader)
    ns_ptr = ccall(
        (:xmlTextReaderConstNamespaceUri, libxml2),
        Cstring,
        (Ptr{Cvoid},),
        getfield(reader, :ptr))
    if ns_ptr == C_NULL
        throw(ArgumentError("no namespace"))
    end
    return unsafe_string(ns_ptr)
end

"""
    expandtree(reader::StreamReader)

Expand the current node of `reader` into a full subtree that will be available
until the next read of node.

Note that the expanded subtree is a read-only and temporary object. You cannot
modify it or keep references to any nodes of it after reading the next node.

Currently, namespace functions and XPath query will not work on the expanded
subtree.
"""
function expandtree(reader::StreamReader)
    node_ptr = @check ccall(
        (:xmlTextReaderExpand, libxml2),
        Ptr{_Node},
        (Ptr{Cvoid},),
        getfield(reader, :ptr)) != C_NULL
    # do not automatically free memories
    return Node(node_ptr, false)
end


text:
function create_model(connectivity)
  N = size(connectivity, 1) # Number of nodes
  E = count(!iszero, connectivity)# Number of edges
  @parameters α[1:N], β[1:N], γ[1:E], κ[1:E], η[1:E]
  @variables t, a[1:N](t)
  rxs = Reaction[]
  for i=1:size(connectivity, 1)
    push!(rxs, Reaction(β[i]*α[i], nothing, [a[i]]))
    push!(rxs, Reaction(β[i], [a[i]], nothing))
  end
  e = 1 # Edge count
  for (i, row) in enumerate(eachrow(connectivity))
    for (j, val) in enumerate(row)
      if val == -1
        push!(rxs, Reaction(hill(abs(a[j]),γ[e],κ[e],η[e]), [a[i]], nothing))
        e += 1
      elseif val == 1
        push!(rxs, Reaction((1.0 - a[i]) * hill(abs(a[j]),γ[e],κ[e],η[e]), nothing, [a[i]]))
        e += 1
      end
    end
  end
  @named model = ReactionSystem(rxs, t)
  return model
end


text:
export bisect, cross_edges

"""
`bisect(G::SimpleGraph)` partitions the vertex set of `G` using the
eigenvector associated with the second smallest eigenvalue of the
graph's Laplacian matrix (called `x` below).

This can be invoked as follows:

+ `bisect(G,"user",pivot)` splits the vertices `v` depending on
  `x[v] >= pivot` vs. `x[v] < pivot`.

+ `bisect(G,"zero")` is the same as `bisect(G,"user", 0.0)`.

+ `bisect(G,"median")` is equivalent to `bisect(G,"user",m)` where `m`
  is the median value of `x`.

+ `bisect(G,"equal")` creates a partition in which the two parts have sizes
  the differ by at most 1.

A plain call to `bisect(G)` is equivalent to `bisect(G,"zero")` (which
is the same as `bisect(G,"user", 0.0)`).
"""
function bisect(G::SimpleGraph,
                where::AbstractString="zero",
                pivot::Real=0.0
               )

    verbose = false

    if verbose
        println("G = $G")
        println("where = $where")
        println("pivot = $pivot")
    end

    T = vertex_type(G)
    VV = vlist(G)
    n  = NV(G)
    L  = laplace(G)
    # x  = collect( eig(L)[2][:,2] )
    x = collect( eigen(L).vectors[:,2])

    if verbose
        println(sort(x))
    end

    piv = 0.0

    if where=="equal"
        pairs = sort(collect(zip(x,VV)))
        vtcs  = [ p[2] for p in pairs ]

        if verbose
            println(pairs)
            println(vtcs)
        end

        mid = floor(Int, n/2)
        if verbose
            println("Equal partition of the vertex set: $mid and $(n-mid)")
        end
        A = Set{T}(vtcs[1:mid])
        B = Set{T}(vtcs[mid+1:end])
        return A,B
    end

    if where=="median"
        piv = median(x)
    elseif where=="user"
        piv = pivot
    elseif where=="zero"
        piv = 0.0
    else
        error("Unknown \"where\" specifier: $where")
    end

    if verbose
        println("pivot set to $piv")
    end

    A = Set{T}()
    B = Set{T}()

    for k=1:n
        v = VV[k]
        if x[k] >= piv
            push!(A,v)
        else
            push!(B,v)
        end
    end

    return A,B
end

# import IterTools.product

"""
`cross_edges(G::SimpleGraph,A,B)` returns the set of edges of `G` with
one end in `A` and one end in `B`. Here `A` and `B` are collections
of vertices of `G`.
"""
function cross_edges(G::SimpleGraph, A, B)

    AB = Base.Iterators.product(A,B)

    result = Set(filter(e -> has(G,e[1],e[2]), collect(AB)))
    return result
end


text:
# Copyright 2020 Google LLC
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

using Plots, LaTeXStrings

function run_algs(test_function::SmoothPlusProximalFunction, v_0::Vector{Float64},
    parameters::Basic_agd_parameters)
  println("Running basic_AGD")
  final_point, basic_history = basic_AGD(test_function,
    v_0,
    parameters)

  println("Running function_value_restart_AGD")
  final_point, fval_restart_history =
    function_value_restart_AGD(test_function, v_0, parameters
     )

  println("Running our_adaptive_scheme_AGD")
  final_point, our_adaptive_history =
    our_adaptive_restart_AGD(test_function, v_0, parameters)

  return basic_history, fval_restart_history, our_adaptive_history
end

function plot_algs(basic_history, fval_restart_history, our_adaptive_history;
  y_min::Float64=1e-8)
 plot(basic_history.fval_history, yaxis=:log, line=(2,:solid),
    label="No restarts", color=1)
 plot_history!(fval_restart_history,
    label="O'D. and C.", color=4)
 plot_history!(our_adaptive_history, label="Our adaptive", color=2)

 f_max = max(
  maximum(basic_history.fval_history),
  maximum(fval_restart_history.fval_history),
  maximum(our_adaptive_history.fval_history),
 )
 ylims!((y_min, 2.0 * f_max))
 plot!(xlabel="Iteration", ylabel=L"$f(x) - f^{*}$")
end

function plot_history!(history::AGD_history;
    label::String="",
    color="black",
    line=(2,:solid),
  )
  plot!(history.fval_history,
    yaxis=:log,
    label=label,
    color=color,
    line=line,
  )
  scatter!(history.restart_intervals,
    history.fval_history[history.restart_intervals],
    yaxis=:log,
    label="",
    color=color,
  )
end

function power_of_two_list(list::Array{Int64,1})
  power_list = Array{Int64,1}()
  for i in list
    push!(power_list, 2^i)
  end
  return power_list
end

function pick_best_fixed_restart_frequency(test_function, v_0, parameters;
    set_of_frequencies::Array{Int64,1},
    approximate_optimality_threshold::Float64)
  best_number_of_iter = Inf
  best_fval = Inf
  best_history = nothing
  best_freq = NaN
  best_point = nothing
  for freq in set_of_frequencies
    final_point, fixed_frequency_history =
      fixed_frequency_AGD(test_function, v_0, parameters, freq)
    this_frequency_is_best = false
    # find the first iteration that is below approximate_optimality_threshold
    lower_iterations = findall(x -> x < approximate_optimality_threshold, fixed_frequency_history.fval_history)
    if length(lower_iterations) > 1
      num_iter_to_approx_opt = lower_iterations[1]
    else
      num_iter_to_approx_opt = length(fixed_frequency_history.fval_history)
    end
    if num_iter_to_approx_opt <= best_number_of_iter
      final_fval = fixed_frequency_history.fval_history[num_iter_to_approx_opt]
      if num_iter_to_approx_opt < best_number_of_iter ||
        final_fval < best_fval
         best_number_of_iter = num_iter_to_approx_opt
         best_fval = final_fval
         best_history = fixed_frequency_history
         best_freq = freq
         best_point = final_point
      end
    end
  end
  if best_freq == maximum(set_of_frequencies)
    @warn "best frequency is equal to maximum"
  end
  if best_freq == minimum(set_of_frequencies)
    @warn "best frequency is equal to minimum"
  end
  return (frequency=best_freq, history=best_history, point=best_point)
end


text:
using SpikingNetworks
@unit

parseParameters"""
τₘ = 20*ms
τₑ = 5*ms
τᵢ = 10*ms
Vₜ = -50*mV
Vᵣ = -60*mV
Eₗ = -49*mV
N=4000;M=4*N÷5
wₑ=zeros(N,N)
wᵢ=zeros(N,N)
wₑ[1:N,1:M]=(rand(N,M).<0.02)*(60*0.27/10)*mV
wᵢ[1:N,M+1:N]=(rand(N,N-M).<0.02)*(-20*4.5/10)*mV"""

parseEquations"""
dv/dt=(gₑ+gᵢ-(v-Eₗ))/τₘ
dgₑ/dt=-gₑ/τₑ
dgᵢ/dt=-gᵢ/τᵢ"""

net=init"""
N=N
v=Vᵣ+rand(N)*(Vₜ-Vᵣ)
gₑ=zeros(N)
gᵢ=zeros(N)
wₑ=wₑ
wᵢ=wᵢ"""

spike"v.>Vₜ"

reset"""
v[i]=Vᵣ"""

synapse"""
gₑ+=wₑ
gᵢ+=wᵢ"""

net.dt=0.1
@time run!(net,100/net.dt)
t,neuron=rasterPlot(net)
using MatlabPlot
mplot(t,neuron,".k")


text:
# ============================================================
# I / O Methods
# ============================================================

export read_cfg, read_dict
export write_sol, write_jld, write_bson, write_vtk, write_tec

include("io_read.jl")
include("io_write.jl")
include("io_plot.jl")


text:
using OctTrees
import OctTrees: modify, stop_cond
using GeometricalPredicates
import GeometricalPredicates:getx, gety, getz
using Test

##################################################################
#
#   testing QUADTREES
#
##################################################################

q=QuadTree(100)

OctTrees.insert!(q, Point(0.1, 0.1))
OctTrees.insert!(q, Point(0.9, 0.9))

@test q.number_of_nodes_used == 5

q=QuadTree(100)

OctTrees.insert!(q, Point(0.1, 0.1))
OctTrees.insert!(q, Point(0.9, 0.9))

tot=0
for i in 1:q.number_of_nodes_used
    !isfullleaf(q.nodes[i]) && continue
    global tot+=1
end
@test tot == 2

@test !q.head.lxhy.is_divided
@test q.head.lxhy.is_empty
@test !q.head.hxly.is_divided
@test q.head.hxly.is_empty
@test !q.head.lxly.is_divided
@test !q.head.lxly.is_empty
@test !q.head.hxhy.is_divided
@test !q.head.hxhy.is_empty
@test q.head.is_divided
@test q.head.is_empty

OctTrees.insert!(q, Point(0.55, 0.9))

@test !q.head.hxhy.hxhy.is_divided
@test !q.head.hxhy.hxhy.is_empty
@test !q.head.hxhy.lxhy.is_divided
@test !q.head.hxhy.lxhy.is_empty
@test !q.head.hxhy.lxly.is_divided
@test q.head.hxhy.lxly.is_empty
@test !q.head.hxhy.hxly.is_divided
@test q.head.hxhy.hxly.is_empty

OctTrees.insert!(q, Point(0.9, 0.55))

@test !q.head.hxhy.hxhy.is_divided
@test !q.head.hxhy.hxhy.is_empty
@test !q.head.hxhy.lxhy.is_divided
@test !q.head.hxhy.lxhy.is_empty
@test !q.head.hxhy.lxly.is_divided
@test q.head.hxhy.lxly.is_empty
@test !q.head.hxhy.hxly.is_divided
@test !q.head.hxhy.hxly.is_empty

##################################################################

struct Part <: AbstractPoint2D
	_x::Float64
	_y::Float64
	Part(x,y) = new(x,y)
end
Part() = Part(0., 0.)
getx(p::Part) = p._x
gety(p::Part) = p._y

q=QuadTree(Part; n=4000100)

pa = [Part(rand(), rand()) for i in 1:1000000]
function insert_unsorted_array(pa::Array{Part,1}, q::QuadTree)
	for p in pa
		OctTrees.insert!(q, p)
	end
end

@time insert_unsorted_array(pa,q)

##################################################################

pa = [Point(rand(), rand()) for i in 1:1000000]
function insert_unsorted_array(pa::Array{Point2D,1}, q::QuadTree)
	for p in pa
		OctTrees.insert!(q, p)
	end
end
q=QuadTree(4000100)
@time insert_unsorted_array(pa,q)
clear!(q)
@time insert_unsorted_array(pa,q)

# a massive particle
struct Particle <: AbstractPoint2D
	_x::Float64
	_y::Float64
	_m::Float64
	Particle(x,y,m) = new(x,y,m)
end
Particle(x::Float64, y::Float64) = Particle(x, y, 1.)
Particle() = Particle(0., 0., 0.)
getx(p::Particle) = p._x
gety(p::Particle) = p._y

q=QuadTree(Particle; n=100)

function modify(q::QuadTreeNode{Particle}, p::Particle)
	total_mass = q.point._m + p._m
	newx = (q.point._x*q.point._m + p._x)/total_mass
	newy = (q.point._y*q.point._m + p._y)/total_mass
	q.point = Particle(newx, newy, total_mass)
end

@test q.head.is_empty == true

OctTrees.insert!(q, Particle(0.1, 0.1), Modify)

@test q.head.is_empty == false
@test q.head.point._m == 1.0
@test q.head.point._x == 0.1
@test q.head.point._y == 0.1

OctTrees.insert!(q, Particle(0.9, 0.9), Modify)

@test q.head.is_empty == true
@test q.head.point._m == 2.0
@test q.head.point._x == (0.1+0.9)/2
@test q.head.point._y == (0.1+0.9)/2
@test q.head.lxly.point._m == 1.0
@test q.head.lxly.point._x == 0.1
@test q.head.lxly.point._y == 0.1
@test q.head.hxhy.point._m == 1.0
@test q.head.hxhy.point._x == 0.9
@test q.head.hxhy.point._y == 0.9

cond_satisfied = false
function stop_cond(q::QuadTreeNode{Particle}, cond_data::Int64)
	q.point._m <= 1.1 && return false
	global cond_satisfied = true
	@test q.point._m == 2.0
	@test cond_data==1
	true
end

OctTrees.map(q, 1)

@test cond_satisfied == true

float_cond_satisfied = false
function stop_cond(q::QuadTreeNode{Particle}, cond_data::Float64)
	q.point._m <= 1.1 && return false
	global float_cond_satisfied = true
	@test q.point._m == 2.0
	@test cond_data==1.0
	true
end

OctTrees.map(q, 1.0)

@test float_cond_satisfied == true

nodata_cond_satisfied = false
function stop_cond(q::QuadTreeNode{Particle})
	q.point._m <= 1.1 && return false
	global nodata_cond_satisfied = true
	@test q.point._m == 2.0
	true
end

OctTrees.map(q)

@test nodata_cond_satisfied == true

q=QuadTree(Particle; n=100)

function modify(q::QuadTreeNode{Particle}, p::Particle, i::Int64)
	@test i==1
	q.point = Particle(q.point._x, q.point._y, 7.0)
end

OctTrees.insert!(q, Particle(0.1, 0.1), 1)
OctTrees.insert!(q, Particle(0.9, 0.9), 1)
@test q.head.point._m == 7.0

N = 10000
q=QuadTree()
for i in 1:N
	OctTrees.insert!(q, Point(rand(), rand()))
end
tot=0
for i in 1:q.number_of_nodes_used
    !isfullleaf(q.nodes[i]) && continue
    global tot+=1
end
@test tot == N

##################################################################
#
#   testing OCTTREES
#
##################################################################

q = OctTree(100)

OctTrees.insert!(q, Point(0.1, 0.1, 0.1))
OctTrees.insert!(q, Point(0.9, 0.9, 0.9))

@test q.number_of_nodes_used == 9

q = OctTree(100)

OctTrees.insert!(q, Point(0.1, 0.1, 0.1))
OctTrees.insert!(q, Point(0.9, 0.9, 0.9))

tot=0
for i in 1:q.number_of_nodes_used
    !isfullleaf(q.nodes[i]) && continue
    global tot+=1
end
@test tot == 2

@test !q.head.lxlylz.is_divided
@test !q.head.lxlylz.is_empty
@test !q.head.lxlyhz.is_divided
@test q.head.lxlyhz.is_empty
@test !q.head.lxhylz.is_divided
@test q.head.lxhylz.is_empty
@test !q.head.lxhyhz.is_divided
@test q.head.lxhyhz.is_empty
@test !q.head.hxlylz.is_divided
@test q.head.hxlylz.is_empty
@test !q.head.hxlyhz.is_divided
@test q.head.hxlyhz.is_empty
@test !q.head.hxhylz.is_divided
@test q.head.hxhylz.is_empty
@test !q.head.hxhyhz.is_divided
@test !q.head.hxhyhz.is_empty
@test q.head.is_divided
@test q.head.is_empty

OctTrees.insert!(q, Point(0.55, 0.9, 0.9))

@test !q.head.hxhyhz.hxhyhz.is_divided
@test !q.head.hxhyhz.hxhyhz.is_empty
@test !q.head.hxhyhz.lxhyhz.is_divided
@test !q.head.hxhyhz.lxhyhz.is_empty
@test !q.head.hxhyhz.lxlyhz.is_divided
@test q.head.hxhyhz.lxlyhz.is_empty
@test !q.head.hxhyhz.hxlyhz.is_divided
@test q.head.hxhyhz.hxlyhz.is_empty

OctTrees.insert!(q, Point(0.9, 0.55, 0.9))

@test !q.head.hxhyhz.hxhyhz.is_divided
@test !q.head.hxhyhz.hxhyhz.is_empty
@test !q.head.hxhyhz.lxhyhz.is_divided
@test !q.head.hxhyhz.lxhyhz.is_empty
@test !q.head.hxhyhz.lxlyhz.is_divided
@test q.head.hxhyhz.lxlyhz.is_empty
@test !q.head.hxhyhz.hxlyhz.is_divided
@test !q.head.hxhyhz.hxlyhz.is_empty

##################################################################

struct Part3D <: AbstractPoint3D
	_x::Float64
	_y::Float64
	_z::Float64
	Part3D(x,y,z) = new(x,y,z)
end
Part3D() = Part3D(0., 0., 0.)
getx(p::Part3D) = p._x
gety(p::Part3D) = p._y
getz(p::Part3D) = p._z

q=OctTree(Part3D; n=4000100)

pa = [Part3D(rand(), rand(), rand()) for i in 1:1000000]
function insert_unsorted_array(pa::Array{Part3D,1}, q::OctTree)
	for p in pa
		OctTrees.insert!(q, p)
	end
end

@time insert_unsorted_array(pa,q)

##################################################################

pa = [Point(rand(), rand(), rand()) for i in 1:1000000]
function insert_unsorted_array(pa::Array{Point3D,1}, q::OctTree)
	for p in pa
		OctTrees.insert!(q, p)
	end
end
q=OctTree(4000100)
@time insert_unsorted_array(pa,q)
clear!(q)
@time insert_unsorted_array(pa,q)

# a massive particle
struct Particle3D <: AbstractPoint3D
	_x::Float64
	_y::Float64
	_z::Float64
	_m::Float64
	Particle3D(x,y,z,m) = new(x,y,z,m)
end
Particle3D(x::Float64, y::Float64, z::Float64) = Particle3D(x, y, z, 1.)
Particle3D() = Particle3D(0., 0., 0., 0.)
getx(p::Particle3D) = p._x
gety(p::Particle3D) = p._y
getz(p::Particle3D) = p._z

q=OctTree(Particle3D; n=100)

function modify(q::OctTreeNode{Particle3D}, p::Particle3D)
	total_mass = q.point._m + p._m
	newx = (q.point._x*q.point._m + p._x)/total_mass
	newy = (q.point._y*q.point._m + p._y)/total_mass
	newz = (q.point._z*q.point._m + p._z)/total_mass
	q.point = Particle3D(newx, newy, newz, total_mass)
end

@test q.head.is_empty == true

OctTrees.insert!(q, Particle3D(0.1, 0.1, 0.1), Modify)

@test q.head.is_empty == false
@test q.head.point._m == 1.0
@test q.head.point._x == 0.1
@test q.head.point._y == 0.1
@test q.head.point._z == 0.1

OctTrees.insert!(q, Particle3D(0.9, 0.9, 0.9), Modify)

@test q.head.is_empty == true
@test q.head.point._m == 2.0
@test q.head.point._x == (0.1+0.9)/2
@test q.head.point._y == (0.1+0.9)/2
@test q.head.point._z == (0.1+0.9)/2

cond_satisfied = false
function stop_cond(q::OctTreeNode{Particle3D}, cond_data::Int64)
	q.point._m <= 1.1 && return false
	global cond_satisfied = true
	@test q.point._m == 2.0
	@test cond_data==1
	true
end

OctTrees.map(q, 1)

@test cond_satisfied == true

float_cond_satisfied = false
function stop_cond(q::OctTreeNode{Particle3D}, cond_data::Float64)
	q.point._m <= 1.1 && return false
	global float_cond_satisfied = true
	@test q.point._m == 2.0
	@test cond_data==1.0
	true
end

OctTrees.map(q, 1.0)

@test float_cond_satisfied == true

nodata_cond_satisfied = false
function stop_cond(q::OctTreeNode{Particle3D})
	q.point._m <= 1.1 && return false
	global nodata_cond_satisfied = true
	@test q.point._m == 2.0
	true
end

OctTrees.map(q)

@test nodata_cond_satisfied == true

q=OctTree(Particle3D; n=100)

function modify(q::OctTreeNode{Particle3D}, p::Particle3D, i::Int64)
	@test i==1
	q.point = Particle3D(q.point._x, q.point._y, q.point._z, 7.0)
end

OctTrees.insert!(q, Particle3D(0.1, 0.1, 0.1), 1)
OctTrees.insert!(q, Particle3D(0.9, 0.9, 0.9), 1)
@test q.head.point._m == 7.0

N = 10000
q=OctTree()
for i in 1:N
	OctTrees.insert!(q, Point(rand(), rand(), rand()))
end
tot=0
for i in 1:q.number_of_nodes_used
    !isfullleaf(q.nodes[i]) && continue
    global tot+=1
end
@test tot == N

#################################################################
#
#   Testing CompiledOctTrees
#
################################################################

N = 10000
q=OctTree()
for i in 1:N
	OctTrees.insert!(q, Point(rand(), rand(), rand()))
end

c = CompiledOctTree(N, Point3D)
compile!(c, q)

for iter in 1:5
    total_number_of_particles=0
    for i in 1:c.number_of_nodes_used
        c.nodes[i].l > 0.0 && continue # not a leaf
        total_number_of_particles += 1
    end
    @test total_number_of_particles == N

    for i in 1:q.number_of_nodes_used
        n = q.nodes[i]
        if isemptyleaf(n)
            @test n.id <= 0
            continue
         end
        @test n.id > 0
        @test n.id <= c.number_of_nodes_used
        @test n.point._x == c.nodes[n.id].point._x
        @test n.point._y == c.nodes[n.id].point._y
        @test n.point._z == c.nodes[n.id].point._z
        if isleaf(n)
            @test c.nodes[n.id].l == -1.0
        else
            @test 2.0*n.r == c.nodes[n.id].l
        end
    end

    v=zeros(Int64, c.number_of_nodes_used)
    for i in 1:q.number_of_nodes_used
        if q.nodes[i].id>0
            v[q.nodes[i].id] = 1
        end
    end
    @test sum(v) == c.number_of_nodes_used

    for i in 1:q.number_of_nodes_used
        n = q.nodes[i]
        !n.is_divided && continue

        !isemptyleaf(n.lxlylz) && @test c.nodes[n.lxlylz.id].next == -1

        !isemptyleaf(n.lxlylz) && !isemptyleaf(n.lxlyhz) &&
            @test c.nodes[n.lxlyhz.id].next == n.lxlylz.id
        !isemptyleaf(n.hxlyhz) && !isemptyleaf(n.hxlylz) &&
            @test c.nodes[n.hxlyhz.id].next == n.hxlylz.id
        !isemptyleaf(n.hxhyhz) && !isemptyleaf(n.hxhylz) &&
            @test c.nodes[n.hxhyhz.id].next == n.hxhylz.id
    end

    clear!(q)
    for i in 1:N
    	OctTrees.insert!(q, Point(rand(), rand(), rand()))
    end
    compile!(c, q)
end


text:
module TestTimedata

using Base.Test
using DataArrays
using DataFrames
using Dates

using TimeData

println("Running Timedata tests")

##########
## find ##
##########

## with NAs present
##-----------------
td = TimeData.testcase(TimeData.Timenum, 2)
boolTd = td .> 130

## from boolean Timedata
@test isequal(find(boolTd), [3, 4, 5])
@test isequal(TimeData.find2sub(boolTd), ([3, 4, 5], [1, 1, 1]))

## using function
@test isequal(find(x -> x .> 130, td), [3, 4, 5])
@test isequal(TimeData.find2sub(x -> x .> 130, td), ([3, 4, 5], [1, 1, 1]))

## no true value
##--------------

td = TimeData.testcase(TimeData.Timenum, 2)
boolTd = td .> 230

## from boolean Timedata
@test isequal(find(boolTd), [])
@test isequal(TimeData.find2sub(boolTd), ([], []))

## using function
@test isequal(find(x -> x .> 230, td), [])
@test isequal(TimeData.find2sub(x -> x .> 230, td), ([], []))

## applied to non-boolean function
##--------------------------------

td = TimeData.testcase(TimeData.Timenum, 2)
@test_throws ErrorException find(x -> x .* 2, td)
@test_throws ErrorException TimeData.find2sub(x -> x .* 2, td)

## applied to numeric Timedata
td = TimeData.testcase(TimeData.Timedata, 2)
@test_throws MethodError find(td)

end


text:
using Base: AbstractCmd
using Mocking: @patch, apply
using TimeZones: TimeZone, localzone

# For mocking make sure we are actually changing the time zone
name = string(localzone()) == "Europe/Warsaw" ? "Pacific/Apia" : "Europe/Warsaw"

tzfile_path = joinpath(TZFILE_DIR, split(name, '/')...)
@assert isfile(tzfile_path) "Tests require a local tzfile present"
win_name = name == "Europe/Warsaw" ? "Central European Standard Time" : "Samoa Standard Time"
tz = TimeZone(name)

if Sys.iswindows()
    patch = @patch read(cmd::AbstractCmd, ::Type{String}) = win_name
    apply(patch) do
        @test localzone() == tz
    end

    # Dateline Standard Time -> Etc/GMT+12 -> UTC-12:00

elseif Sys.isunix()
    withenv("TZ" => nothing) do
        if Sys.islinux()
            # Determine time zone from /etc/timezone
            patches = [
                @patch isfile(f::AbstractString) = f == "/etc/timezone"
                @patch open(fn::Function, f::AbstractString) = fn(IOBuffer("$name #Works with comments\n"))
            ]
            apply(patches) do
                @test localzone() == tz
            end

            # Determine time zone from /etc/conf.d/clock
            patches = [
                @patch isfile(f::AbstractString) = f == "/etc/conf.d/clock"
                @patch open(fn::Function, f::AbstractString) = fn(IOBuffer("\n\nTIMEZONE=\"$name\""))
            ]
            apply(patches) do
                @test localzone() == tz
            end
        end

        # Determine time zone from symlink /etc/localtime
        patches = [
            @patch isfile(f::AbstractString) = false
            @patch islink(f::AbstractString) = f == "/etc/localtime"
            @patch readlink(f::AbstractString) = "/usr/share/zoneinfo/$name"
        ]
        apply(patches) do
            @test localzone() == tz
        end

        # Determine time zone from contents of /etc/localtime
        tz_from_file = open(tzfile_path) do f
            TimeZones.read_tzfile(f, "local")
        end

        patches = [
            @patch isfile(f::AbstractString) = f == "/etc/localtime"
            @patch islink(f::AbstractString) = false
            @patch open(fn::Function, f::AbstractString) = fn(open(tzfile_path))
        ]
        apply(patches) do
            @test localzone() == tz_from_file
        end

        # Unable to determine time zone
        patches = [
            @patch isfile(f::AbstractString) = false
            @patch islink(f::AbstractString) = false
        ]
        apply(patches) do
            @test_throws ErrorException localzone()
        end
    end
end

# Generates a cross-platform patched environment to run `localzone()`
function with_localzone(func::Function, name::AbstractString)
    @static if Sys.isapple()
        patches = [
            @patch islink(f::AbstractString) = f == "/etc/localtime"
            @patch readlink(::AbstractString) = "/usr/share/zoneinfo/$name"
        ]
    elseif Sys.iswindows()
        patches = [
            @patch read(cmd::AbstractCmd, ::Type{String}) = name
        ]
    elseif Sys.islinux()
        patches = [
            @patch isfile(f::AbstractString) = f == "/etc/timezone"
            @patch open(fn::Function, f::AbstractString) = fn(IOBuffer("$name\n"))
        ]
    else
        error("Unhandled OS")
    end

    withenv("TZ" => nothing) do
        apply(patches) do
            func()
        end
    end
end

# https://github.com/JuliaTime/TimeZones.jl/issues/154
@testset "legacy time zones" begin
    # "US/Pacific" is deprecated in favor of "America/Los_Angeles"
    @test istimezone("US/Pacific", Class(:LEGACY))
    name = Sys.isunix() ? "US/Pacific" : "Pacific Standard Time"
    with_localzone(name) do
        @test localzone().transitions == tz"America/Los_Angeles".transitions
    end

    # "America/Montreal" is deprecated in favor of "America/Toronto"
    @test istimezone("America/Montreal", Class(:LEGACY))
    if Sys.isunix()
        with_localzone("America/Montreal") do
            @test localzone().transitions == tz"America/Toronto".transitions
        end
    end

    # "America/Indianapolis" is deprecated in favor of "America/Indiana/Indianapolis"
    # Note: On Windows "US Eastern Standard Time" translates to "America/Indianapolis"
    @test istimezone("America/Indianapolis", Class(:LEGACY))
    name = Sys.isunix() ? "America/Indianapolis" : "US Eastern Standard Time"
    with_localzone(name) do
        @test localzone().transitions == tz"America/Indiana/Indianapolis".transitions
    end
end


text:
using ScalarRadau
using Test

#stiff test function
F(x, y, _) = 50*(cos(x) - y)
#exact solution
S(x) = (50*(-50*exp(-50*x) + sin(x) + 50*cos(x)))/2501

#test with evenly spaced output
x, y = radau(F, 0.0, 0.0, 6.0, 100, atol=1e-9, rtol=1e-9)
@test maximum(@. abs(y - S(x))) < 1e-6

#test with end-point output only
yₙ = radau(F, 0.0, 0.0, 6.0, atol=1e-9, rtol=1e-9)
@test abs(yₙ - S(6)) < 1e-6

#test for convergence failure
@test_throws ErrorException radau((x,y,p)->rand(), 0.0, 0.0, 1.0, atol=1e-9, rtol=1e-9)

#test for max steps failure
@test_throws ErrorException radau(F, 0.0, 0.0, 1e6, atol=1e-9, rtol=1e-9, maxstep=1)