text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

module Pkg

import Random
import REPL

export @pkg_str
export PackageSpec
export PackageMode, PKGMODE_MANIFEST, PKGMODE_PROJECT
export UpgradeLevel, UPLEVEL_MAJOR, UPLEVEL_MAJOR, UPLEVEL_MINOR, UPLEVEL_PATCH
export Registry, RegistrySpec

depots() = Base.DEPOT_PATH
function depots1()
    d = depots()
    isempty(d) && Pkg.Types.pkgerror("no depots found in DEPOT_PATH")
    return d[1]
end

logdir() = joinpath(depots1(), "logs")
devdir() = get(ENV, "JULIA_PKG_DEVDIR", joinpath(depots1(), "dev"))
envdir(depot = depots1()) = joinpath(depot, "environments")
const UPDATED_REGISTRY_THIS_SESSION = Ref(false)

# load snapshotted dependencies
include("../ext/TOML/src/TOML.jl")

include("GitTools.jl")
include("PlatformEngines.jl")
include("Types.jl")
include("Compress.jl")
include("Display.jl")
include("Pkg2/Pkg2.jl")
include("GraphType.jl")
include("Resolve.jl")
include("Operations.jl")
include("API.jl")
include("Registry.jl")
include("REPLMode.jl")

import .REPLMode: @pkg_str
import .Types: UPLEVEL_MAJOR, UPLEVEL_MINOR, UPLEVEL_PATCH, UPLEVEL_FIXED
import .Types: PKGMODE_MANIFEST, PKGMODE_PROJECT
# legacy CI script support
import .API: clone, dir

"""
    PackageMode

An enum with the instances

  * `PKGMODE_MANIFEST`
  * `PKGMODE_PROJECT`

Determines if operations should be made on a project or manifest level.
Used as an argument to  [`PackageSpec`](@ref) or as an argument to [`Pkg.rm`](@ref).
"""
const PackageMode = Types.PackageMode

"""
    UpgradeLevel

An enum with the instances

  * `UPLEVEL_FIXED`
  * `UPLEVEL_PATCH`
  * `UPLEVEL_MINOR`
  * `UPLEVEL_MAJOR`

Determines how much a package is allowed to be updated.
Used as an argument to  [`PackageSpec`](@ref) or as an argument to [`Pkg.update`](@ref).
"""
const UpgradeLevel = Types.UpgradeLevel

# Define new variables so tab comleting Pkg. works.
"""
    Pkg.add(pkg::Union{String, Vector{String}})
    Pkg.add(pkg::Union{PackageSpec, Vector{PackageSpec}})

Add a package to the current project. This package will be available using the
`import` and `using` keywords in the Julia REPL and if the current project is
a package, also inside that package.

# Examples
```julia
Pkg.add("Example") # Add a package from registry
Pkg.add(PackageSpec(name="Example", version="0.3")) # Specify version
Pkg.add(PackageSpec(url="https://github.com/JuliaLang/Example.jl", rev="master")) # From url
Pkg.add(PackageSpec(url="/remote/mycompany/juliapackages/OurPackage"))` # From path (has to be a gitrepo)
```

See also [`PackageSpec`](@ref).
"""
const add = API.add

"""
    Pkg.rm(pkg::Union{String, Vector{String}})
    Pkg.rm(pkg::Union{PackageSpec, Vector{PackageSpec}})

Remove a package from the current project. If the `mode` of `pkg` is
`PKGMODE_MANIFEST` also remove it from the manifest including all
recursive dependencies of `pkg`.

See also [`PackageSpec`](@ref), [`PackageMode`](@ref).
"""
const rm = API.rm

"""
    Pkg.update(; level::UpgradeLevel=UPLEVEL_MAJOR, mode::PackageMode = PKGMODE_PROJECT)
    Pkg.update(pkg::Union{String, Vector{String}})
    Pkg.update(pkg::Union{PackageSpec, Vector{PackageSpec}})

Update a package `pkg`. If no posistional argument is given, update all packages in the manifest if `mode` is `PKGMODE_MANIFEST` and packages in both manifest and project if `mode` is `PKGMODE_PROJECT`.
If no positional argument is given `level` can be used to control what how much packages are allowed to be upgraded (major, minor, patch, fixed).

See also [`PackageSpec`](@ref), [`PackageMode`](@ref), [`UpgradeLevel`](@ref).
"""
const update = API.up

"""
    Pkg.test(; coverage::Bool=true)
    Pkg.test(pkg::Union{String, Vector{String}; coverage::Bool=true)
    Pkg.test(pkgs::Union{PackageSpec, Vector{PackageSpec}}; coverage::Bool=true)

Run the tests for package `pkg` or if no positional argument is given to `test`,
the current project is tested (which thus needs to be a package).
A package is tested by running its `test/runtests.jl` file.

The tests are run by generating a temporary environment with only `pkg` and its (recursive) dependencies
(recursively) in it. If a manifest exist, the versions in that manifest is used, otherwise
a feasible set of package are resolved and installed.

During the test, test-specific dependencies are active, which are
given in the project file as e.g.

```
[extras]
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[targets]
test = ["Test"]
```

Coverage statistics for the packages may be generated by
passing `coverage=true`. The default behavior is not to run coverage.

The tests are executed in a new process with `check-bounds=yes` and by default `startup-file=no`.
If using the startup file (`~/.julia/config/startup.jl`) is desired, start julia with `--startup-file=yes`.
Inlining of functions during testing can be disabled (for better coverage accuracy)
by starting julia with `--inline=no`.
"""
const test = API.test

"""
    Pkg.gc()

Garbage collect packages that are no longer reachable from any project.
Only packages that are tracked by version are deleted, so no packages
that might contain local changes are touched.
"""
const gc = API.gc

"""
    Pkg.build(; verbose = false)
    Pkg.build(pkg::Union{String, Vector{String}}; verbose = false)
    Pkg.build(pkgs::Union{PackageSpec, Vector{PackageSpec}}; verbose = false)

Run the build script in `deps/build.jl` for `pkg` and all of the dependencies in
depth-first recursive order.
If no argument is given to `build`, the current project is built, which thus needs
to be a package.
This function is called automatically one any package that gets installed
for the first time.
`verbose = true` prints the build output to `stdout`/`stderr` instead of
redirecting to the `build.log` file.
"""
const build = API.build

# TODO: decide what to do with this
const installed = API.installed

"""
    Pkg.pin(pkg::Union{String, Vector{String}})
    Pkg.pin(pkgs::Union{Packagespec, Vector{Packagespec}})

Pin a package to the current version (or the one given in the `packagespec` or a certain
git revision. A pinned package is never updated.
"""
const pin = API.pin

"""
    Pkg.free(pkg::Union{String, Vector{String}})
    Pkg.free(pkgs::Union{Packagespec, Vector{Packagespec}})

Free a package which removes a `pin` if it exists, or if the package is tracking a path,
e.g. after [`Pkg.develop`](@ref), go back to tracking registered versions.

# Examples
```julia
Pkg.free("Package")
Pkg.free(PackageSpec("Package"))
```
"""
const free = API.free

"""
    Pkg.develop(pkg::Union{String, Vector{String}})
    Pkg.develop(pkgs::Union{Packagespec, Vector{Packagespec}})

Make a package available for development by tracking it by path.
If `pkg` is given with only a name or by a URL the packages will be downloaded
to the location by the environment variable `JULIA_PKG_DEVDIR` with
`.julia/dev` as the default.

If `pkg` is given as a local path, the package at that path will be tracked.

# Examples
```julia
# By name
Pkg.develop("Example")

# By url
Pkg.develop(PackageSpec(url="https://github.com/JuliaLang/Compat.jl"))

# By path
Pkg.develop(PackageSpec(path="MyJuliaPackages/Package.jl")
```

See also [`PackageSpec`](@ref)

"""
const develop = API.develop

#TODO: Will probably be deprecated for something in PkgDev
const generate = API.generate

"""
    Pkg.instantiate()

If a `Manifest.toml` file exist in the current project, download all
the packages declared in that manifest.
Else, resolve a set of feasible packages from the `Project.toml` files
and install them.
"""
const instantiate = API.instantiate

"""
    Pkg.resolve()

Update the current manifest with eventual changes to the dependency graph
from packages that are tracking a path.
"""
const resolve = API.resolve

"""
    Pkg.status([pkgs...]; mode::PackageMode=PKGMODE_PROJECT)

Print out the status of the project/manifest.
If `mode` is `PKGMODE_PROJECT` prints out status about only those packages
that are in the project (explicitly added). If `mode` is `PKGMODE_MANIFEST`
also print for those in the manifest (recursive dependencies). If there are
any packages listed as arguments the output will be limited to those packages.

!!! compat "Julia 1.1"
    `Pkg.status` with package arguments requires at least Julia 1.1.
"""
const status = API.status

"""
    Pkg.activate([s::String]; shared::Bool=false)

Activate the environment at `s`. The active environment is the environment
that is modified by executing package commands.
The logic for what path is activated is as follows:

  * If `shared` is `true`, the first existing environment named `s` from the depots
    in the depot stack will be activated. If no such environment exists,
    create and activate that environment in the first depot.
  * If `s` is an existing path, then activate the environment at that path.
  * If `s` is a package in the current project and `s` is tracking a path, then
    activate the environment at the tracked path.
  * Else, `s` is interpreted as a non-existing path, activate that path.

If no argument is given to `activate`, then activate the home project.
The home project is specified by either the `--project` command line option to
the julia executable, or the `JULIA_PROJECT` environment variable.

# Examples
```
Pkg.activate()
Pkg.activate("local/path")
Pkg.activate("MyDependency")
```
"""
const activate = API.activate

"""
    PackageSpec(name::String, [uuid::UUID, version::VersionNumber])
    PackageSpec(; name, url, path, rev, version, mode, level)

A `PackageSpec` is a representation of a package with various metadata.
This includes:

  * The `name` of the package.
  * The package unique `uuid`.
  * A `version` (for example when adding a package. When upgrading, can also be an instance of
   the enum [`UpgradeLevel`](@ref)
  * A `url` and an optional git `rev`ision. `rev` could be a branch name or a git commit SHA.
  * A local path `path`. This is equivalent to using the `url` argument but can be more descriptive.
  * A `mode`, which is an instance of the enum [`PackageMode`](@ref) which can be either `PKGMODE_PROJECT` or
   `PKGMODE_MANIFEST`, defaults to `PKGMODE_PROJECT`. Used in e.g. [`Pkg.rm`](@ref).

Most functions in Pkg take a `Vector` of `PackageSpec` and do the operation on all the packages
in the vector.

Below is a comparison between the REPL version and the `PackageSpec` version:

| `REPL`               | `API`                                                 |
|:---------------------|:------------------------------------------------------|
| `Package`            | `PackageSpec("Package")`                              |
| `Package@0.2`        | `PackageSpec(name="Package", version="0.2")`          |
| `Package=a67d...`    | `PackageSpec(name="Package", uuid="a67d...")`         |
| `Package#master`     | `PackageSpec(name="Package", rev="master")`           |
| `local/path#feature` | `PackageSpec(path="local/path"; rev="feature")`       |
| `www.mypkg.com`      | `PackageSpec(url="www.mypkg.com")`                    |
| `--manifest Package` | `PackageSpec(name="Package", mode=PKGSPEC_MANIFEST)`  |
| `--major Package`    | `PackageSpec(name="Package", version=PKGLEVEL_MAJOR)` |
"""
const PackageSpec = Types.PackageSpec

"""
    Pkg.setprotocol!(proto::Union{Nothing, AbstractString}=nothing)

Set the protocol used to access GitHub-hosted packages when `add`ing a url or `develop`ing a package.
Defaults to delegating the choice to the package developer (`proto == nothing`).
Other choices for `proto` are `"https"` or `"git"`.
"""
const setprotocol! = API.setprotocol!

"""
    RegistrySpec(name::String)
    RegistrySpec(; name, url, path)

A `RegistrySpec` is a representation of a registry with various metadata, much like
[`PackageSpec`](@ref).

Most registry functions in Pkg take a `Vector` of `RegistrySpec` and do the operation
on all the registries in the vector.

!!! compat "Julia 1.1"
    Pkg's registry handling requires at least Julia 1.1.

# Examples

Below is a comparison between the REPL version and the `RegistrySpec` version:

| `REPL`               | `API`                                           |
|:---------------------|:------------------------------------------------|
| `Registry`           | `RegistrySpec("Registry")`                      |
| `Registry=a67d...`   | `RegistrySpec(name="Registry", uuid="a67d..."`  |
| `local/path`         | `RegistrySpec(path="local/path")`               |
| `www.myregistry.com` | `RegistrySpec(url="www.myregistry.com")`        |
"""
const RegistrySpec = Types.RegistrySpec

function __init__()
    if isdefined(Base, :active_repl)
        REPLMode.repl_init(Base.active_repl)
    else
        atreplinit() do repl
            if isinteractive() && repl isa REPL.LineEditREPL
                isdefined(repl, :interface) || (repl.interface = REPL.setup_interface(repl))
                REPLMode.repl_init(repl)
            end
        end
    end
end

METADATA_compatible_uuid(pkg::String) = Types.uuid5(Types.uuid_package, pkg)

##################
# Precompilation #
##################

const CTRL_C = '\x03'
const precompile_script = """
    import Pkg
    tmp = mktempdir()
    cd(tmp)
    empty!(DEPOT_PATH)
    pushfirst!(DEPOT_PATH, tmp)
    # Prevent cloning registry
    mkdir("registries")
    touch("registries/blocker") # prevents default registry from cloning
    touch("Project.toml")
    ] activate .
    $CTRL_C
    Pkg.add("Test") # adding an stdlib doesn't require internet access
    ] add Te\t\t$CTRL_C
    ] st
    $CTRL_C
    rm(tmp; recursive=true)"""

end # module


text:
#!/usr/bin/julia

# Copyright (C) 2021, 2022 Matheus Fernandes Bigolin <mfrdrbigolin@disroot.org>
# SPDX-License-Identifier: MIT

"Day Six, Lanternfish."

include("Utils.jl")

using DelimitedFiles
using .Utils

"Calculate the number of fishes after `days` given a `fishpop`."
function solve(fishpop :: Array{Int}, days :: Int64)
    ages = zeros(Int, 9)

    for fish in fishpop
        ages[fish + 1] += 1
    end

    for i in 1:days
        # This cycles modularly through the ages.
        (start, final) = (mod1(i, 9), mod1(i + 7, 9))

        ages[final] += ages[start]
    end

    return sum(ages)
end

function main()
    usage_and_exit(length(ARGS) != 1)

    data = readdlm(first(ARGS), ',', Int)

    println(solve(data, 80))
    println(solve(data, 256))
end

if abspath(PROGRAM_FILE) == @__FILE__
    main()
end

ext: jl


text:
module NormalizeQuantiles

export normalizeQuantiles

export sampleRanks
export qnTiesMethods,tmMin,tmMax,tmOrder,tmReverse,tmRandom,tmAverage

if VERSION < v"0.4.0-"
	macro doc(string)
	end
	macro enum(x...)
	end
end # if VERSION < v"0.4.0-"

@enum qnTiesMethods tmMin tmMax tmOrder tmReverse tmRandom tmAverage

if VERSION >= v"0.4.0-"

@doc "
### qnmatrix::Array{Float64} function normalizeQuantiles(matrix::Array{Float64})

Method for input type Array{Float64}
" ->
function normalizeQuantiles(matrix::Array{Float64})
    damatrix=Array{Nullable{Float64}}(matrix)
    r=normalizeQuantiles(damatrix)
	convert(Array{Float64},reshape([get(r[i]) for i=1:length(r)],size(r)))
end

@doc "
### qnmatrix::Array{Float64} function normalizeQuantiles(matrix::Array{Int})

Method for input type Array{Int}
" ->
function normalizeQuantiles(matrix::Array{Int})
    dafloat=Array{Nullable{Float64}}(convert(Array{Float64},matrix))
    r=normalizeQuantiles(dafloat)
	convert(Array{Float64},reshape([get(r[i]) for i=1:length(r)],size(r)))
end

@doc "
### qnmatrix::Array{Nullable{Float64}} function normalizeQuantiles(matrix::Array{Nullable{Int}})

Method for input type Array{Nullable{Int}}
" ->
function normalizeQuantiles(matrix::Array{Nullable{Int}})
    nullable=Array{Nullable{Float64}}(matrix)
    normalizeQuantiles(nullable)
end

@doc "
### qnmatrix::Array{Nullable{Float64}} function normalizeQuantiles(matrix::Array{Nullable{Float64}})
Calculate the quantile normalized data for the input matrix

Parameter:
    matrix::Array{Nullable{Float64}}
The input data as a Array{Nullable{Float64}} of float values interpreted as Array{Nullable{Float64}}(rows,columns)

Return value:
    qnmatrix::Array{Nullable{Float64}}
The quantile normalized data as Array{Nullable{Float64}}

Examples:

    using NormalizeQuantiles

    array = [ 3.0 2.0 1.0 ; 4.0 5.0 6.0 ; 9.0 7.0 8.0 ; 5.0 2.0 8.0 ]

    a = Array{Nullable{Float64}}(array)

    qn = normalizeQuantiles(a)

    column = 2

    row = 2

    a[row,column] = Nullable{Float64}()

    qn = normalizeQuantiles(a)

" ->
function normalizeQuantiles(matrix::Array{Nullable{Float64}})
    nrows=size(matrix,1)
    ncols=size(matrix,2)
	# preparing the result matrix
    qnmatrix=Array{Nullable{Float64}}((nrows,ncols))
	if ncols>0 && nrows>0
		# foreach column: sort the values without NAs; put NAs (if any) back into sorted list
		sortColumns(matrix,qnmatrix,nrows,ncols)
		# foreach row: set all values to the mean of the row, except NAs
		meanRows(qnmatrix,nrows)
		# foreach column: equal values in original column should all be mean of normalized values
		# foreach column: reorder the values back to the original order
		equalValuesInColumnAndOrderToOriginal(matrix,qnmatrix,nrows,ncols)
	end
    qnmatrix
end

@doc "
### qnmatrix::SharedArray{Float64} function normalizeQuantiles(matrix::SharedArray{Float64})

Method for input type SharedArray{Float64}
" ->
function normalizeQuantiles(matrix::SharedArray{Float64})
    nullable=SharedArray(Nullable{Float64},(size(matrix,1),size(matrix,2)))
	nullable[:]=matrix[:]
    r=normalizeQuantiles(nullable)
	nullable=SharedArray(Nullable{Float64},(0,0))
	ra=SharedArray(Float64,(size(r,1),size(r,2)))
	ra[:]=convert(Array{Float64},reshape([get(r[i]) for i=1:length(r)],size(r)))[:]
	r=SharedArray(Nullable{Float64},(0,0))
	ra
end

@doc "
### qnmatrix::SharedArray{Float64} function normalizeQuantiles(matrix::SharedArray{Int})

Method for input type SharedArray{Int}
" ->
function normalizeQuantiles(matrix::SharedArray{Int})
    nullable=SharedArray(Nullable{Float64},(size(matrix,1),size(matrix,2)))
	nullable[:]=[ Float64(x) for x in matrix ][:]
    r=normalizeQuantiles(nullable)
	nullable=SharedArray(Nullable{Float64},(0,0))
	ra=SharedArray(Float64,(size(r,1),size(r,2)))
	ra[:]=convert(Array{Float64},reshape([get(r[i]) for i=1:length(r)],size(r)))[:]
	r=SharedArray(Nullable{Float64},(0,0))
	ra
end

@doc "
### qnmatrix::SharedArray{Nullable{Float64}} function normalizeQuantiles(matrix::SharedArray{Nullable{Int}})

Method for input type SharedArray{Int}
" ->
function normalizeQuantiles(matrix::SharedArray{Nullable{Int}})
	nullable=SharedArray(Nullable{Float64},(size(matrix,1),size(matrix,2)))
	nullable[:]=matrix[:]
	normalizeQuantiles(nullable)
end

@doc "
### qnmatrix::SharedArray{Nullable{Float64}} function normalizeQuantiles(matrix::SharedArray{Nullable{Float64}})

Quantile normalization using multiple cores (see ?normalizeQuantiles)

Example:

	addprocs()

	using NormalizeQuantiles

    array = [ 3.0 2.0 1.0 ; 4.0 5.0 6.0 ; 9.0 7.0 8.0 ; 5.0 2.0 8.0 ]

    a = Array{Nullable{Float64}}(array)

	sa = SharedArray(Nullable{Float64},(size(a,1),size(a,2)))

	sa[:]=a[:]

    qn = normalizeQuantiles(sa)

" ->
function normalizeQuantiles(matrix::SharedArray{Nullable{Float64}})
    nrows=size(matrix,1)
    ncols=size(matrix,2)
	# preparing the result matrix
    qnmatrix=SharedArray(Nullable{Float64},(nrows,ncols))
	if ncols>0 && nrows>0
		# foreach column: sort the values without NAs; put NAs (if any) back into sorted list
		multicoreSortColumns(matrix,qnmatrix,nrows,ncols)
		# foreach row: set all values to the mean of the row, except NAs
		multicoreMeanRows(qnmatrix,nrows,ncols)
		# foreach column: equal values in original column should all be mean of normalized values
		# foreach column: reorder the values back to the original order
		multicoreEqualValuesInColumnAndOrderToOriginal(matrix,qnmatrix,nrows,ncols)
	end
    qnmatrix
end

function Old_normalizeQuantiles(matrix::Array{Nullable{Float64}})
    nrows=size(matrix,1)
    ncols=size(matrix,2)
	# preparing the result matrix
    qnmatrix=Array{Nullable{Float64}}((nrows,ncols))
	if ncols>0 && nrows>0
		# foreach column: sort the values without NAs; put NAs (if any) back into sorted list
		sortColumns(matrix,qnmatrix,nrows,ncols)
		# foreach row: set all values to the mean of the row, except NAs
		meanRows(qnmatrix,nrows)
		# foreach column: equal values in original column should all be mean of normalized values
		equalValuesInColumn(matrix,qnmatrix,nrows,ncols)
		# foreach column: reorder the values back to the original order
		orderToOriginal(matrix,qnmatrix,nrows,ncols)
	end
    qnmatrix
end

function sortColumns(matrix::Array{Nullable{Float64}},qnmatrix::Array{Nullable{Float64}},nrows,ncols)
	for column = 1:ncols
		indices=[ !isnull(x) for x in vec(matrix[:,column]) ]
		sortcol=vec(matrix[:,column])[indices]
		sortcol=[Float64(get(x)) for x in sortcol]
		sort!(sortcol)
		sortcol=Array{Nullable{Float64}}(sortcol)
		naindices=(1:nrows)[!indices]
		empty=sum(indices)==0
		for napos in naindices
			if empty
				empty = false
				sortcol=Array{Nullable{Float64}}([Nullable{Float64}()])
			else
				sortcol=vcat(sortcol[1:napos-1],Nullable{Float64}(),sortcol[napos:end])
			end
		end
		qnmatrix[:,column]=sortcol
	end
end

function multicoreSortColumns(matrix::SharedArray{Nullable{Float64}},qnmatrix::SharedArray{Nullable{Float64}},nrows,ncols)
	@sync @parallel for column = 1:ncols
		indices=[ !isnull(x) for x in vec(matrix[:,column]) ]
		sortcol=vec(matrix[:,column])[indices]
		sortcol=[Float64(get(x)) for x in sortcol]
		sort!(sortcol)
		sortcol=Array{Nullable{Float64}}(sortcol)
		naindices=(1:nrows)[!indices]
		empty=sum(indices)==0
		for napos in naindices
			if empty
				empty = false
				sortcol=Array{Nullable{Float64}}([Nullable{Float64}()])
			else
				sortcol=vcat(sortcol[1:napos-1],Nullable{Float64}(),sortcol[napos:end])
			end
		end
		qnmatrix[:,column]=sortcol
	end
end

function meanRows(qnmatrix::Array{Nullable{Float64}},nrows)
	for row = 1:nrows
		indices=[ !isnull(x) for x in qnmatrix[row,:] ]
		nacount=sum(!indices)
		rowmean=mean([Float64(get(x)) for x in qnmatrix[row,indices]])
		qnmatrix[row,:]=Nullable{Float64}(rowmean)
		nacount>0?qnmatrix[row,!indices]=Nullable{Float64}():false
	end
end

function multicoreMeanRows(qnmatrix::SharedArray{Nullable{Float64}},nrows,ncols)
	@sync @parallel for row = 1:nrows
		indices=[ !isnull(x) for x in qnmatrix[row,:] ]
		nacount=sum(!indices)
		rowmean=mean([Float64(get(x)) for x in qnmatrix[row,indices]])
		qnmatrix[row,:]=Nullable{Float64}(rowmean)
		nacount>0?qnmatrix[row,!indices]=Nullable{Float64}():false
	end
end

function equalValuesInColumn(matrix::Array{Nullable{Float64}},qnmatrix::Array{Nullable{Float64}},nrows,ncols)
	for column = 1:ncols
		indices=[ !isnull(x) for x in matrix[:,column] ]
		sortp=[ Float64(get(x)) for x in vec(matrix[:,column])[indices] ]
		sortp=sortperm(sortp)
		nacount=sum(!indices)
		sortcol=matrix[:,column][indices][sortp]
		ranks=Array(Int,(length(sortcol),1))
		if nacount < nrows
			ranks[1]=1
			lastrank=1
			lastvalue=sortcol[1]
			for i in 2:length(sortcol)
				nextvalue=sortcol[i]
				get(nextvalue)==get(lastvalue)?ranks[i]=ranks[i-1]:ranks[i]=ranks[i-1]+1
				lastrank=ranks[i]
				lastvalue=sortcol[i]
			end
			indices=1:nrows
			naindices=[ isnull(x) for x in qnmatrix[:,column] ]
			indices=indices[!naindices]
			for i in 1:lastrank
				values=i.==vec(ranks)
				qnmatrix[indices[values],column]=mean([ get(x) for x in qnmatrix[indices[values],column] ])
			end
		end
	end
end

function orderToOriginal(matrix::Array{Nullable{Float64}},qnmatrix::Array{Nullable{Float64}},nrows,ncols)
	for column = 1:ncols
		indices=[ !isnull(x) for x in vec(matrix[:,column]) ]
		sortp=[ Float64(get(x)) for x in vec(matrix[:,column])[indices] ]
		sortp=sortperm(sortp)
		indices2=[ !isnull(x) for x in vec(qnmatrix[:,column]) ]
		qncol=Array{Nullable{Float64}}(nrows,1)
		fill!(qncol,Nullable{Float64}())
		qncol[(1:nrows)[indices][sortp]]=vec(qnmatrix[(1:nrows)[indices2],column])
		qnmatrix[:,column]=qncol
	end
end

function multicoreEqualValuesInColumnAndOrderToOriginal(matrix::SharedArray{Nullable{Float64}},qnmatrix::SharedArray{Nullable{Float64}},nrows,ncols)
	@sync @parallel for column = 1:ncols
		indices=[ !isnull(x) for x in vec(matrix[:,column]) ]
		sortp=[ Float64(get(x)) for x in vec(matrix[:,column])[indices] ]
		sortp=sortperm(sortp)
		indices2=[ !isnull(x) for x in vec(qnmatrix[:,column]) ]
		if length(matrix[:,column][indices][sortp])>0
			allranks=Dict{Int,Array{Int}}()
			sizehint!(allranks,nrows)
			rankColumns=getRankMatrix(matrix[:,column][indices][sortp],allranks,indices)
			for i in 1:rankColumns
				rankIndices=unique(allranks[i])
				qnmatrix[rankIndices,column]=mean([ get(x) for x in qnmatrix[rankIndices,column] ])
			end
		end
		qncol=Array{Nullable{Float64}}(nrows,1)
		fill!(qncol,Nullable{Float64}())
		qncol[(1:nrows)[indices][sortp]]=vec(qnmatrix[(1:nrows)[indices2],column])
		qnmatrix[:,column]=qncol
	end
end

function equalValuesInColumnAndOrderToOriginal(matrix::Array{Nullable{Float64}},qnmatrix::Array{Nullable{Float64}},nrows,ncols)
	for column = 1:ncols
		indices=[ !isnull(x) for x in vec(matrix[:,column]) ]
		sortp=[ Float64(get(x)) for x in vec(matrix[:,column])[indices] ]
		sortp=sortperm(sortp)
		indices2=[ !isnull(x) for x in vec(qnmatrix[:,column]) ]
		if length(matrix[:,column][indices][sortp])>0
			allranks=Dict{Int,Array{Int}}()
			sizehint!(allranks,nrows)
			rankColumns=getRankMatrix(matrix[:,column][indices][sortp],allranks,indices)
			for i in 1:rankColumns
				rankIndices=unique(allranks[i])
				qnmatrix[rankIndices,column]=mean([ get(x) for x in qnmatrix[rankIndices,column] ])
			end
		end
		qncol=Array{Nullable{Float64}}(nrows,1)
		fill!(qncol,Nullable{Float64}())
		qncol[(1:nrows)[indices][sortp]]=vec(qnmatrix[(1:nrows)[indices2],column])
		qnmatrix[:,column]=qncol
	end
end

function getRankMatrix(sortedArrayNoNAs::Array{Nullable{Float64}},allranks::Dict{Int,Array{Int}},indices::Array{Bool})
	rankColumns=0
	nrows=length(sortedArrayNoNAs)
	lastvalue=sortedArrayNoNAs[1]
	indices2=(1:length(indices))[indices]
	count=1
	for i in 2:nrows
		nextvalue=sortedArrayNoNAs[i]
		if !isnull(nextvalue) && !isnull(lastvalue) && get(nextvalue)==get(lastvalue)
			if haskey(allranks,rankColumns+1)
				allranks[rankColumns+1]=vcat(allranks[rankColumns+1],Array{Int}([indices2[i-1],indices2[i]]))
			else
				allranks[rankColumns+1]=Array{Int}([indices2[i-1],indices2[i]])
			end
			count+=1
		else
			if count>1
				rankColumns+=1
			end
			count=1
		end
		lastvalue=sortedArrayNoNAs[i]
	end
	if count>1
		rankColumns+=1
	end
	rankColumns
end

@doc "
### (Array{Int},Dict{Int,Array{Int}}) sampleRanks(array::Array{Int};tiesMethod::qnTiesMethods=tmMin,naIncreasesRank=false,resultMatrix=false)

" ->
function sampleRanks(array::Array{Int};tiesMethod::qnTiesMethods=tmMin,naIncreasesRank=false,resultMatrix=false)
	nullable=Array{Nullable{Float64}}(Array{Float64}(array))
	(rn,m)=sampleRanks(nullable,tiesMethod=tiesMethod,naIncreasesRank=naIncreasesRank,resultMatrix=resultMatrix)
	r=convert(Array{Int},reshape([get(rn[i]) for i=1:length(rn)],size(rn)))
	(r,m)
end

@doc "
### (Array{Int},Dict{Int,Array{Int}}) sampleRanks(array::Array{Float64};tiesMethod::qnTiesMethods=tmMin,naIncreasesRank=false,resultMatrix=false)

" ->
function sampleRanks(array::Array{Float64};tiesMethod::qnTiesMethods=tmMin,naIncreasesRank=false,resultMatrix=false)
	nullable=Array{Nullable{Float64}}(array)
	(rn,m)=sampleRanks(nullable,tiesMethod=tiesMethod,naIncreasesRank=naIncreasesRank,resultMatrix=resultMatrix)
	r=convert(Array{Int},reshape([get(rn[i]) for i=1:length(rn)],size(rn)))
	(r,m)
end

@doc "
### (Array{Nullable{Int}},Dict{Int,Array{Int}}) sampleRanks(array::Array{Nullable{Int}};tiesMethod::qnTiesMethods=tmMin,naIncreasesRank=false,resultMatrix=false)

" ->
function sampleRanks(array::Array{Nullable{Int}};tiesMethod::qnTiesMethods=tmMin,naIncreasesRank=false,resultMatrix=false)
	nullable=Array{Nullable{Float64}}(array)
	sampleRanks(nullable,tiesMethod=tiesMethod,naIncreasesRank=naIncreasesRank,resultMatrix=resultMatrix)
end

@doc "
### (Array{Nullable{Int}},Dict{Int,Array{Int}}) sampleRanks(array::Array{Nullable{Float64}};tiesMethod::qnTiesMethods=tmMin,naIncreasesRank=false,resultMatrix=false)

Calculate ranks of the values of a given vector.

Parameters:

	array: the input array
	tiesMethod: the method how ties (equal values) are treated
	   possible values: tmMin tmMax tmOrder tmReverse tmRandom tmAverage
	   default is tmMin
	naIncreasesRank: if true than any NA increases the following ranks by 1
	resultMatrix: if true than return a dictionary of rank keys and array of indices values

Example:

	using NormalizeQuantiles

    a = [ 5.0 2.0 4.0 3.0 1.0 ]

	n = Array{Nullable{Float64}}(a)

    (r,m)=sampleRanks(n)

    (r,m)=sampleRanks(n,tiesMethod=tmMin,naIncreasesRank=false,resultMatrix=true)

r is the vector of ranks.

m is a dictionary with rank as keys and as value the indices of all values of this rank.

" ->
function sampleRanks(array::Array{Nullable{Float64}};tiesMethod::qnTiesMethods=tmMin,naIncreasesRank=false,resultMatrix=false)
	nrows=length(array)
	indices=[ !isnull(x) for x in array ]
	reducedArray=[ Float64(get(x)) for x in array[indices] ]
	sortp=sortperm(reducedArray)
	result=Array{Nullable{Int}}(nrows)
	result[:]=Nullable{Int}()
	#resultMatrix?begin rankMatrix=Dict{Int,Array{Int}}();sizehint!(rankMatrix,nrows) end:rankMatrix=nothing
	rankMatrix=Dict{Int,Array{Int}}()
	if resultMatrix
		sizehint!(rankMatrix,nrows)
	end
	indices2=(1:nrows)[indices][sortp]
	rank=1
	narank=0
	lastvalue=reducedArray[sortp[1]]
	ties=Array{Int}(0)
	tieIndices=Array{Int}(0)
	tiesCount=0
	index=1
	for i in 1:(nrows+1)
		last=i>nrows
		if !last
			newvalue=reducedArray[sortp[index]]
		end
		if !last && !indices[i]
			if naIncreasesRank
				rank+=1
				tiesCount>0?narank+=1:false
			end
		else
			if last || newvalue != lastvalue
				if tiesMethod==tmMin
					ties[:]=minimum(ties)
					rank=ties[end]+narank+1
				elseif tiesMethod==tmMax
					ties[:]=maximum(ties)
					rank=ties[end]+narank+1
				elseif tiesMethod==tmReverse
					ties=flipdim(ties,1)
					rank=ties[1]+narank+1
				elseif tiesMethod==tmRandom
					rank=ties[end]+narank+1
					ties=ties[randperm(tiesCount)]
				elseif tiesMethod==tmAverage
					ties[:]=round(Int,mean(ties))
					rank=ties[end]+narank+1
				end
				narank=0
				for j in 1:tiesCount
					if resultMatrix
						if haskey(rankMatrix,ties[j])
							rankMatrix[ties[j]]=vcat(rankMatrix[ties[j]],tieIndices[j])
						else
							rankMatrix[ties[j]]=Array{Int}([tieIndices[j]])
						end
					end
					result[tieIndices[j]]=ties[j]
				end
				ties=Array{Int}(0)
				tieIndices=Array{Int}(0)
				tiesCount=0
			end
			if !last
				tieIndices=vcat(tieIndices,Array{Int}([indices2[index]]))
				ties=vcat(ties,Array{Int}([rank]))
				tiesCount+=1
				rank+=1
				lastvalue=newvalue
				index+=1
			end
		end
	end
	(result,rankMatrix)
end

end # if VERSION >= v"0.4.0-"

if VERSION < v"0.4.0-"

using DataArrays

"""
### qnmatrix::DataArray{Float64} function normalizeQuantiles(matrix::Array{Float64})

Method for input type Array{Float64}
"""
function normalizeQuantiles(matrix::Array{Float64})
    damatrix=DataArray(matrix)
    r=normalizeQuantiles(damatrix)
	convert(Array{Float64},reshape([r[i] for i=1:length(r)],size(r)))
end

"""
### qnmatrix::DataArray{Float64} function normalizeQuantiles(matrix::Array{Int})

Method for input type Array{Int}
"""
function normalizeQuantiles(matrix::Array{Int})
    dafloat=DataArray(convert(Array{Float64},matrix))
    r=normalizeQuantiles(dafloat)
	convert(Array{Float64},reshape([r[i] for i=1:length(r)],size(r)))
end

"""
### qnmatrix::DataArray{Float64} function normalizeQuantiles(matrix::DataArray{Int})

Method for input type DataArray{Int}
"""
function normalizeQuantiles(matrix::DataArray{Int})
    dafloat=convert(DataArray{Float64},matrix)
    normalizeQuantiles(dafloat)
end

"""
### qnmatrix::DataArray{Float64} function normalizeQuantiles(matrix::DataArray{Float64})
Calculate the quantile normalized data for the input matrix

Parameter:
    matrix::DataArray{Float64}
The input data as a DataArray of float values interpreted as DataArray(rows,columns)

Return value:
    qnmatrix::DataArray{Float64}
The quantile normalized data as DataArray{Float64}

Examples:

    using NormalizeQuantiles

    using DataArrays

    array = [ 3.0 2.0 1.0 ; 4.0 5.0 6.0 ; 9.0 7.0 8.0 ; 5.0 2.0 8.0 ]

    da = DataArray(array)

    qn = normalizeQuantiles(da)

    column = 2

    row = 2

    da[row,column] = NA

    qn = normalizeQuantiles(da)

"""
function normalizeQuantiles(matrix::DataArray{Float64})
    nrows=size(matrix,1)
    ncols=size(matrix,2)
	# preparing the result matrix
    qnmatrix=DataArray(Float64,(nrows,ncols))
	if ncols>0 && nrows>0
		# foreach column: sort the values without NAs; put NAs (if any) back into sorted list
		sortColumns(matrix,qnmatrix,nrows,ncols)
		# foreach row: set all values to the mean of the row, except NAs
		meanRows(qnmatrix,nrows,ncols)
		# foreach column: equal values in original column should all be mean of normalized values
		# foreach column: reorder the values back to the original order
		equalValuesInColumnAndOrderToOriginal(matrix,qnmatrix,nrows,ncols)
	end
    qnmatrix
end

"""
### qnmatrix::SharedArray{Float64} function normalizeQuantiles(matrix::SharedArray{Int})

Method for input type SharedArray{Int}
"""
function normalizeQuantiles(matrix::SharedArray{Int})
    sa=SharedArray(Float64,(size(matrix,1),size(matrix,2)))
	sa[:]=matrix[:]
    normalizeQuantiles(sa)
end

"""
### qnmatrix::SharedArray{Float64} function normalizeQuantiles(matrix::SharedArray{Float64})

Quantile normalization using of multiple cores (see ?normalizeQuantiles)

Example:

	addprocs()

	using NormalizeQuantiles

    a = [ 3.0 2.0 1.0 ; 4.0 5.0 6.0 ; 9.0 7.0 8.0 ; 5.0 2.0 8.0 ]

	sa = SharedArray(Float64,(size(a,1),size(a,2)))

	sa[:]=a[:]

    qn = normalizeQuantilesMultiCore(sa)

"""
function normalizeQuantiles(matrix::SharedArray{Float64})
    nrows=size(matrix,1)
    ncols=size(matrix,2)
	# preparing the result matrix
    qnmatrix=SharedArray(Float64,(nrows,ncols))
	if ncols>0 && nrows>0
		# foreach column: sort the values without NAs; put NAs (if any) back into sorted list
		multicoreSortColumns(matrix,qnmatrix,nrows,ncols)
		# foreach row: set all values to the mean of the row, except NAs
		multicoreMeanRows(qnmatrix,nrows,ncols)
		# foreach column: equal values in original column should all be mean of normalized values
		# foreach column: reorder the values back to the original order
		multicoreEqualValuesInColumnAndOrderToOriginal(matrix,qnmatrix,nrows,ncols)
	end
    qnmatrix
end

function OLD_normalizeQuantiles(matrix::DataArray{Float64})
    nrows=size(matrix,1)
    ncols=size(matrix,2)
	# preparing the result matrix
    qnmatrix=DataArray(Float64,(nrows,ncols))
	if ncols>0 && nrows>0
		# foreach column: sort the values without NAs; put NAs (if any) back into sorted list
		sortColumns(matrix,qnmatrix,nrows,ncols)
		# foreach row: set all values to the mean of the row, except NAs
		meanRows(qnmatrix,nrows,ncols)
		# foreach column: equal values in original column should all be mean of normalized values
		equalValuesInColumn(matrix,qnmatrix,nrows,ncols)
		# foreach column: reorder the values back to the original order
		orderToOriginal(matrix,qnmatrix,nrows,ncols)
	end
    qnmatrix
end

function sortColumns(matrix::DataArray{Float64},qnmatrix::DataArray{Float64},nrows,ncols)
	for column = 1:ncols
		indices=[ !isa(x,NAtype) for x in vec(matrix[:,column]) ]
		sortcol=vec(matrix[:,column])[indices]
		sort!(sortcol)
		naindices=(1:nrows)[!indices]
		empty=sum(indices)==0
		for napos in naindices
			if empty
				empty = false
				sortcol=[NA]
			else
				sortcol=vcat(sortcol[1:napos-1],NA,sortcol[napos:end])
			end
		end
		qnmatrix[:,column]=sortcol
	end
end

function multicoreSortColumns(matrix::SharedArray{Float64},qnmatrix::SharedArray{Float64},nrows,ncols)
	@sync @parallel for column = 1:ncols
		sortp=sortperm(vec(matrix[:,column]))
		sortcol=matrix[:,column][sortp]
		qnmatrix[:,column]=sortcol
	end
end

function meanRows(qnmatrix::DataArray{Float64},nrows,ncols)
	for row = 1:nrows
		naindices=[ isa(x,NAtype) for x in qnmatrix[row,:] ]
		nacount=sum(naindices)
		indices=(1:ncols)[!naindices]
		naindices=(1:ncols)[naindices]
		qnmatrix[row,:]=mean(qnmatrix[row,:][indices])
		nacount>0?qnmatrix[row,naindices]=NA:false
	end
end

function multicoreMeanRows(qnmatrix::SharedArray{Float64},nrows,ncols)
	@sync @parallel for row = 1:nrows
		indices=(1:ncols)
		qnmatrix[row,:]=mean(qnmatrix[row,:][indices])
	end
end

function equalValuesInColumn(matrix::DataArray{Float64},qnmatrix::DataArray{Float64},nrows,ncols)
	for column = 1:ncols
		sortp=sortperm(vec(matrix[:,column]))
		naindices=[ isa(x,NAtype) for x in matrix[:,column][sortp] ]
		nacount=sum(naindices)
		sortp=sortp[!naindices]
		sortcol=matrix[:,column][sortp]
		ranks=Array(Int,(length(sortcol),1))
		if nacount < nrows
			ranks[1]=1
			lastrank=1
			lastvalue=sortcol[1]
			for i in 2:length(sortcol)
				nextvalue=sortcol[i]
				nextvalue==lastvalue?ranks[i]=ranks[i-1]:ranks[i]=ranks[i-1]+1
				lastrank=ranks[i]
				lastvalue=sortcol[i]
			end
			indices=1:nrows
			naindices=[ isa(x,NAtype) for x in qnmatrix[:,column] ]
			indices=indices[!naindices]
			for i in 1:lastrank
				values=i.==vec(ranks)
				qnmatrix[indices[values],column]=mean(qnmatrix[indices[values],column])
			end
		end
	end
end

function orderToOriginal(matrix::DataArray{Float64},qnmatrix::DataArray{Float64},nrows,ncols)
	for column = 1:ncols
		sortp=sortperm(vec(matrix[:,column]))
		naindices=[ isa(x,NAtype) for x in matrix[:,column][sortp] ]
		naindices2=[ isa(x,NAtype) for x in qnmatrix[:,column] ]
		nacount=sum(naindices)
		qncol=vec(qnmatrix[:,column])
		for i in 1:length(sortp[!naindices])
			qncol[sortp[!naindices][i]]=vec(qnmatrix[:,column])[!naindices2][i]
		end
		for i in 1:length(sortp[naindices])
			qncol[sortp[naindices][i]]=NA
		end
		qnmatrix[:,column]=qncol
	end
end

function multicoreEqualValuesInColumnAndOrderToOriginal(matrix::SharedArray{Float64},qnmatrix::SharedArray{Float64},nrows,ncols)
	@sync @parallel for column = 1:ncols
		indices=[ true for x in matrix[:,column] ]
		sortp=matrix[:,column][indices]
		sortp=sortperm(sortp)
		indices2=[ true for x in qnmatrix[:,column] ]
		if length(matrix[:,column][indices][sortp])>0
			allranks=Dict{Int,Array{Int}}()
			sizehint!(allranks,nrows)
			rankColumns=getRankMatrix(matrix[:,column][indices][sortp],allranks,indices)
			for i in 1:rankColumns
				rankIndices=unique(allranks[i])
				qnmatrix[rankIndices,column]=mean(qnmatrix[rankIndices,column])
			end
		end
		qncol=vec(qnmatrix[:,column])
		qncol[(1:nrows)[indices][sortp]]=vec(qnmatrix[(1:nrows)[indices2],column])
		qnmatrix[:,column]=qncol
	end
end

function equalValuesInColumnAndOrderToOriginal(matrix::DataArray{Float64},qnmatrix::DataArray{Float64},nrows,ncols)
	for column = 1:ncols
		indices=[ !isa(x,NAtype) for x in matrix[:,column] ]
		sortp=matrix[:,column][indices]
		sortp=sortperm(sortp)
		indices2=[ !isa(x,NAtype) for x in qnmatrix[:,column] ]
		if length(matrix[:,column][indices][sortp])>0
			allranks=Dict{Int,Array{Int}}()
			sizehint!(allranks,nrows)
			rankColumns=getRankMatrix(matrix[:,column][indices][sortp],allranks,indices)
			for i in 1:rankColumns
				rankIndices=unique(allranks[i])
				qnmatrix[rankIndices,column]=mean(qnmatrix[rankIndices,column])
			end
		end
		qncol=vec(qnmatrix[:,column])
		fill!(qncol,NA)
		qncol[(1:nrows)[indices][sortp]]=vec(qnmatrix[(1:nrows)[indices2],column])
		qnmatrix[:,column]=qncol
	end
end

function getRankMatrix(sortedArrayNoNAs::Array{Float64},allranks::Dict{Int,Array{Int}},indices::Array{Bool})
	rankColumns=0
	nrows=length(sortedArrayNoNAs)
	lastvalue=sortedArrayNoNAs[1]
	indices2=(1:length(indices))[indices]
	count=1
	for i in 2:nrows
		nextvalue=sortedArrayNoNAs[i]
		if nextvalue==lastvalue
			if haskey(allranks,rankColumns+1)
				allranks[rankColumns+1]=vcat(allranks[rankColumns+1],[indices2[i-1],indices2[i]])
			else
				allranks[rankColumns+1]=[indices2[i-1],indices2[i]]
			end
			count+=1
		else
			if count>1
				rankColumns+=1
			end
			count=1
		end
		lastvalue=sortedArrayNoNAs[i]
	end
	if count>1
		rankColumns+=1
	end
	rankColumns
end

function getRankMatrix(sortedArrayNoNAs::DataArray{Float64},allranks::Dict{Int,Array{Int}},indices::Array{Bool})
	rankColumns=0
	nrows=length(sortedArrayNoNAs)
	lastvalue=sortedArrayNoNAs[1]
	indices2=(1:length(indices))[indices]
	count=1
	for i in 2:nrows
		nextvalue=sortedArrayNoNAs[i]
		if !isa(nextvalue,NAtype) && !isa(lastvalue,NAtype) && nextvalue==lastvalue
			if haskey(allranks,rankColumns+1)
				allranks[rankColumns+1]=vcat(allranks[rankColumns+1],[indices2[i-1],indices2[i]])
			else
				allranks[rankColumns+1]=[indices2[i-1],indices2[i]]
			end
			count+=1
		else
			if count>1
				rankColumns+=1
			end
			count=1
		end
		lastvalue=sortedArrayNoNAs[i]
	end
	if count>1
		rankColumns+=1
	end
	rankColumns
end

end # if VERSION < v"0.4.0-"

end # module


text:
#the bump function reaches its maximum for t = shape*scale (obtained by differentiating the function)
@inline _bumpvalmaximum(shape::AbstractFloat,scale::AbstractFloat) = shape*scale

###Calculate a bump value for a single time point
@inline function _bumpval(amplitude::AbstractFloat,shape::AbstractFloat,scale::AbstractFloat,timepoint::AbstractFloat)
    @fastmath amplitude*exp(shape*(1 - log(shape) + log(timepoint/scale)) - timepoint/scale)
end

###Calculate the bump (current) which results from a single photon
@inline function _bump!{T<:AbstractFloat}(bumpvals::Vector{T},timepoints::Vector{T},amplitude::T,shape::T,scale::T)
    @simd for i=1:length(timepoints)
        @inbounds bumpvals[i] = _bumpval(amplitude,shape,scale,timepoints[i])
    end
    bumpvals
end

bump!{T<:AbstractFloat}(bumpvals::Vector{T},timepoints::Vector{T},amplitude::T,shape::T,scale::T) = _bump!(bumpvals,timepoints,amplitude,shape,scale)
bump{T<:AbstractFloat}(timepoints::Vector{T},amplitude::T,shape::T,scale::T) = _bump!(Vector{T}(length(timepoints)),timepoints,amplitude,shape,scale)

###Function that determines how much timepoints to include in the bump calculation, based on the bump parameters
function bump(amplitude::AbstractFloat,shape::AbstractFloat,scale::AbstractFloat,cutoff::AbstractFloat)
    @assert cutoff > 0
    t = trunc(_bumpvalmaximum(shape,scale))*2+1 #starting from twice the maximum time of the bump function
    while _bumpval(amplitude,shape,scale,t) > cutoff
        t+=one(t)
    end
    timepoints = collect(one(t):one(t):t-one(t))
    _bump!(Vector{typeof(t)}(length(timepoints)),timepoints,amplitude,shape,scale)
end

###Calculate the macrocurrent generated by a sequence of photons (internal function)
@inline function _macrocurrent!{T<:AbstractFloat}(current::Vector{T},photons::Vector,bumpvals::Vector{T},nsteps::Int,bumplength::Int)
    for i = 1:nsteps-bumplength
        @inbounds p = photons[i]
        if (p>0)
            @simd for j = 1:bumplength
                @inbounds current[i+j-1] += p*bumpvals[j]
            end
        end
    end
    for i=nsteps-bumplength+1:nsteps
        @inbounds p = photons[i]
        if (p>0)
            @simd for j = 1:nsteps-i+1
                @inbounds current[i+j-1] += p*bumpvals[j]
            end
        end
    end
    current
end

function macrocurrent!{T<:AbstractFloat}(current::Vector{T},photons::Vector,bumpvals::Vector{T})
  @assert length(photons) == length(current)
  _macrocurrent!(current,photons,bumpvals,length(photons),length(bumpvals))
end

macrocurrent{T<:AbstractFloat}(photons::Vector,bumpvals::Vector{T}) = _macrocurrent!(zeros(T,length(photons)),photons,bumpvals,length(photons),length(bumpvals))

@inline function _filterphotons!{T<:AbstractFloat}(filteredphotons::Vector{Int},refractimes::Vector{Int},receptor::AbstractPhotoReceptor,
                                                   latency::Vector{T},refractory::Vector{T},bumplength::Int)
    pcounter = 1
    bl = bumplength - 1
    for t=1:receptor.numsteps
        for microvillus in microvilliwithphotons(receptor,t)
            @inbounds if refractimes[microvillus] < t
                lat = t + trunc(Int,latency[pcounter])
                if lat <= receptor.numsteps
                    @inbounds filteredphotons[lat] += 1
                    @inbounds refractimes[microvillus] = lat + bl + trunc(Int,refractory[pcounter])
                else
                    refractimes[microvillus] = receptor.numsteps
                end
                pcounter+=1
            end
        end
    end
    filteredphotons
end

function filterphotons!{T<:AbstractFloat}(filteredphotons::Vector{Int},receptor::AbstractPhotoReceptor,latency::Vector{T},refractory::Vector{T},bumplength::Int)
    _filterphotons!(filteredphotons,zeros(Int,receptor.numvilli),receptor,latency,refractory,bumplength)

end

function filterphotons{T<:AbstractFloat}(receptor::AbstractPhotoReceptor,latency::Vector{T},refractory::Vector{T},bumplength::Int)
    _filterphotons!(zeros(Int,receptor.numsteps),zeros(Int,receptor.numvilli),receptor,latency,refractory,bumplength)
end

@inline function _lightinducedcurrent!{T<:AbstractFloat}(current::Vector{T},filteredphotons::Vector{Int},refractimes::Vector{Int},receptor::AbstractPhotoReceptor,
                                                         latency::Vector{T},refractory::Vector{T},bumpvals::Vector{T})
    bumplength = length(bumpvals)
    _filterphotons!(filteredphotons,refractimes,receptor,latency,refractory,bumplength)
    _macrocurrent!(current,filteredphotons,bumpvals,receptor.numsteps,bumplength)
end

function lightinducedcurrent!{T<:AbstractFloat}(current::Vector{T},receptor::AbstractPhotoReceptor,latency::Vector{T},refractory::Vector{T},bumpvals::Vector{T})
    _lightinducedcurrent!(current,zeros(Int,receptor.numsteps),zeros(Int,receptor.numvilli),receptor,latency,refractory,bumpvals)
end

function lightinducedcurrent{T<:AbstractFloat}(receptor::AbstractPhotoReceptor,latency::Vector{T},refractory::Vector{T},bumpvals::Vector{T})
    _lightinducedcurrent!(zeros(T,receptor.numsteps),zeros(Int,receptor.numsteps),zeros(Int,receptor.numvilli),receptor,latency,refractory,bumpvals)
end


text:
##### Beginning of file

import ..package_directory

function _get_git_binary_path()::String
    deps_jl_file_path = package_directory("deps", "deps.jl")
    if !isfile(deps_jl_file_path)
        delayederror(
            string(
                "MirrorUpdater.jl is not properly installed. ",
                "Please run\nPkg.build(\"MirrorUpdater\")",
                )
            )
    end
    include(deps_jl_file_path)
    git::String = strip(string(git_cmd))
    run(`$(git) --version`)
    @debug(
        "git command: ",
        git,
        )
    return git
end

function git_version()::VersionNumber
    git::String = _get_git_binary_path()
    a::String = convert(String,read(`$(git) --version`, String))
    b::String = convert(String, strip(a))
    c::Vector{SubString{String}} = split(b, "git version")
    d::String = convert(String,last(c))
    e::String = convert(String, strip(d))
    f::VersionNumber = VersionNumber(e)
    return f
end

function clean_up_branch_name(x::String)::String
    temp::String = strip(
        strip(
            strip(
                strip(
                    strip(x,),
                    '*',
                    ),
                ),
            '*',
            ),
        )
    my_regex::Regex = r"[a-zA-Z0-9_-]*\/[a-zA-Z0-9_-]*\/([a-zA-Z0-9_-]*)"
    if occursin(my_regex, temp)
        my_match::RegexMatch = match(my_regex, temp)
        just_the_branch::String =
            clean_up_branch_name(first(my_match.captures))
        result = just_the_branch
    else
        result = temp
    end
    return result
end

clean_up_branch_name(x::AbstractString) = clean_up_branch_name(
    convert(String, x)
    )

function get_all_branches_local()
    git::String = _get_git_binary_path()
    a::String = read(`$(git) branch`, String)
    b::String = convert(String, strip(a))
    c::Vector{SubString{String}} = split(b, '\n')
    d::Vector{String} = clean_up_branch_name.(c)
    e::Vector{String} = sort(unique(d))
    return e
end

function get_all_branches_local_and_remote()::Vector{String}
    git::String = _get_git_binary_path()
    a::String = read(`$(git) branch -a`, String)
    b::String = convert(String, strip(a))
    c::Vector{SubString{String}} = split(b, '\n')
    d::Vector{String} = clean_up_branch_name.(c)
    e::Vector{String} = sort(unique(d))
    return e
end

function get_current_branch()::String
    git::String = _get_git_binary_path()
    a::String = read(`$(git) status`, String)
    b::String = convert(String, strip(a))
    c::Vector{SubString{String}} = split(b, '\n')
    d::String = convert(String, strip(first(c)))
    my_regex::Regex = r"On branch ([a-zA-Z0-9_-]*)"
    if occursin(my_regex, d)
        my_match::RegexMatch = match(my_regex, d)
        just_the_branch::String =
            clean_up_branch_name(first(my_match.captures))
        return just_the_branch
    else
        delayederror("could not determine current branch")
    end
end

function checkout_branch!(
    branch_name::AbstractString;
    create::Bool = false,
    error_on_failure::Bool = true,
    )::Nothing

    success::Bool = false

    git::String = _get_git_binary_path()

    branch_name_cleaned::String = clean_up_branch_name(
        branch_name
        )

    try
        run(`$(git) checkout $(branch_name_cleaned)`)
    catch e1
        @warn(string("ignoring exception"), e1,)
    end

    current_branch_1::String = get_current_branch()
    if strip(current_branch_1) == strip(branch_name_cleaned)
        success = true
    else
        success = false
    end

    if !success
        if create
            try
                run(`$(git) checkout --orphan $(branch_name_cleaned)`)
            catch e2
                @warn(string("ignoring exception"), e2,)
            end
        end
    end

    current_branch_2::String = get_current_branch()
    if strip(current_branch_2) == strip(branch_name_cleaned)
        success = true
    else
        success = false
    end

    if !success
        if error_on_failure
            delayederror("could not checkout the specified branch")
        else
            @warn("could not checkout the specified branch")
        end
    end

    return nothing
end

function branch_exists(branch_name::AbstractString)::Bool
    git::String = _get_git_binary_path()
    original_branch::String = get_current_branch()
    branch_name_cleaned::String = clean_up_branch_name(
        branch_name
        )
    try
        run(`$(git) checkout $(branch_name_cleaned)`)
    catch e
        @warn(string("ignoring exception"), e,)
    end
    current_branch::String = get_current_branch()
    if strip(current_branch)==strip(branch_name_cleaned)
        result = true
    else
        result = false
    end
    run(`$(git) checkout $(original_branch)`)
    return result
end

function git_add_all!()::Nothing
    git::String = _get_git_binary_path()
    try
        run(`$(git) add -A`)
    catch e
    @warn(string("ignoring exception"), e,)
    end
    return nothing
end

function git_commit!(
        ;
        message::AbstractString,
        committer_name::AbstractString,
        committer_email::AbstractString,
        allow_empty::Bool = false,
        )::Nothing
    git::String = _get_git_binary_path()
    message_stripped::String = convert(
        String,
        strip(message),
        )
    committer_name_stripped::String = convert(
        String,
        strip(committer_name),
        )
    committer_email_stripped::String = convert(
        String,
        strip(committer_email),
        )
    allow_empty_flag::String = ""
    if allow_empty
        allow_empty_flag = "--allow-empty"
    else
        allow_empty_flag = ""
    end
    run(`$(git) config user.name "$(committer_name_stripped)"`)
    run(`$(git) config user.email "$(committer_email_stripped)"`)
    run(`$(git) config commit.gpgsign false`)
    try
        run(`$(git) commit $(allow_empty_flag) -m "$(message_stripped)"`)
    catch e
        @warn(string("ignoring exception"), e,)
    end
    return nothing
end

function git_push_upstream_all!()::Nothing
    git::String = _get_git_binary_path()
    try
        run(`$(git) push -u --all`)
    catch e
    @warn(string("ignoring exception"), e,)
    end
    return nothing
end

function delete_everything_except_dot_git!(
        root_path::AbstractString,
        )::Nothing
    previous_directory::String = pwd()

    list_of_paths_to_remove::Vector{String} = Vector{String}()

    root_path_stripped::String = convert(String, strip(root_path))

    for file_or_directory in readdir(root_path_stripped)
        if strip(lowercase(file_or_directory)) != ".git"
            push!(
                list_of_paths_to_remove,
                joinpath(root_path_stripped, file_or_directory,),
                )
        end
    end

    for path_to_remove in list_of_paths_to_remove
        rm(
            path_to_remove;
            force = true,
            recursive = true,
            )
    end

    cd(previous_directory)
    return nothing
end

function delete_only_dot_git!(root_path::AbstractString)::Nothing
    previous_directory::String = pwd()

    list_of_paths_to_remove::Vector{String} = Vector{String}()

    root_path_stripped::String = convert(String, strip(root_path))

    for (rootdir, dirs, files) in walkdir(root_path_stripped)
        for dir in dirs
            if strip(lowercase(dir)) == ".git"
                push!(
                    list_of_paths_to_remove,
                    joinpath(rootdir, dir),
                    )
            end
        end
        for file in files
            if strip(lowercase(file)) == ".git"
                push!(
                    list_of_paths_to_remove,
                    joinpath(rootdir, file),
                    )
            end
        end
    end

    for path_to_remove in list_of_paths_to_remove
        rm(
            path_to_remove;
            force = true,
            recursive = true,
            )
    end

    cd(previous_directory)
    return nothing
end

function _include_branch(
        ;
        include::AbstractVector,
        branch::String,
        )::Bool
    branch_clean::String = clean_up_branch_name(branch)
    matches_any_inclusion_criteria::Vector{Bool} = Bool[
        occursin(x, branch_clean) for x in include
        ]
    result::Bool = any(matches_any_inclusion_criteria)
    return result
end

function _exclude_branch(
        ;
        exclude::AbstractVector,
        branch::String,
        )::Bool
    branch_clean::String = clean_up_branch_name(branch)
    matches_any_exclusion_criteria::Vector{Bool} = Bool[
        occursin(x, branch_clean) for x in exclude
        ]
    result::Bool = any(matches_any_exclusion_criteria)
    return result
end

function make_list_of_branches_to_snapshot(
        ;
        default_branch::String,
        include::AbstractVector,
        exclude::AbstractVector,
        )::Vector{String}

    all_branches::Vector{String} = get_all_branches_local_and_remote()

    first_pass::Vector{String} = Vector{String}()

    default_br_cl_lc_s::String = strip(
        lowercase(
            clean_up_branch_name(
                default_branch
                )
            )
        )

    for b_1 in all_branches
        if default_br_cl_lc_s == strip(lowercase(clean_up_branch_name(b_1)))
            push!(first_pass, b_1)
        end
    end

    for b_2 in all_branches
        if _include_branch(;branch=b_2,include=include)
            push!(first_pass, b_2)
        end
    end

    branches_to_snapshot::Vector{String} = Vector{String}()

    for b_3 in first_pass
        if !_exclude_branch(;branch=b_3,exclude=exclude)
            push!(branches_to_snapshot, b_3)
        end
    end

    branches_to_snapshot_cleaned::Vector{String} = clean_up_branch_name.(
        branches_to_snapshot
        )

    result::Vector{String} = sort(unique(branches_to_snapshot_cleaned))

    return result
end

##### End of file


text:
function CZ(ALPHA,BETA,EL)
      A= [0.770;   0.241;  -0.100;  -0.416;  -0.731; -1.053; -1.366;
          -1.646; -1.917; -2.120; -2.248; -2.229]
      #
      S= 0.2 * ALPHA
      K= floor(S)+3
      if(K <= 1)
          K= 2
      end
      if(K >=  12)
          K=  11
      end
      DA= S - K + 3
      L = K + fix( 1.1*sign(DA) )
      K=convert(Int64, K)
      L=convert(Int64, L)
      S= A[K] + abs(DA) * (A[L] - A[K])
      return czv = S*(1-(BETA/57.3)^2) - .19*(EL/25.0)
end


text:
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the linear advection equation

advectionvelocity = 1.0
equations = LinearScalarAdvectionEquation1D(advectionvelocity)

# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux
solver = DGSEM(polydeg=2, surface_flux=flux_lax_friedrichs)

coordinates_min = (-1.0,) # minimum coordinate
coordinates_max = (1.0,) # maximum coordinate
cells_per_dimension = (16,)

function initial_condition_shock(x, t, equation::LinearScalarAdvectionEquation1D)
    # Store translated coordinate for easy use of exact solution
    x_trans = x - equation.advectionvelocity * t
    if x[1] < 0.75 && x[1]>0.25
      scalar = 1
    else
      scalar = 0
    end
    return SVector(scalar)
  end

# Create curved mesh with 16 cells
mesh = CurvedMesh(cells_per_dimension, coordinates_min, coordinates_max)

# A semidiscretization collects data structures and functions for the spatial discretization
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_shock, solver)

###############################################################################
# ODE solvers, callbacks etc.

# Create ODE problem with time span from 0.0 to 1.0
ode = semidiscretize(semi, (0.0, 1.0));

# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup
# and resets the timers
summary_callback = SummaryCallback()

# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results
analysis_callback = AnalysisCallback(semi, interval=100)

# The SaveSolutionCallback allows to save the solution to a file in regular intervals
save_solution = SaveSolutionCallback(interval=100,
                                     solution_variables=cons2prim)

# The StepsizeCallback handles the re-calculcation of the maximum Δt after each time step
stepsize_callback = StepsizeCallback(cfl=0.5)

# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver
callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, stepsize_callback)

###############################################################################
# run the simulation

# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks
sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
            save_everystep=false, callback=callbacks);

# Print the timer summary
summary_callback()


text:
function AbstractFilter{T<:FFTW.fftwReal,G<:LineGroups}(ψmeta::ΨMeta,
        spec::Spec1D{T,FourierDomain{1},G,Morlet})
    """1. **Gaussian denominator `den = 2σ²`**
    The FWHM (full width at half maximum) bw of a Gaussian bell curve of
    variance `σ` is defined by the equation
        `g(±bw/2) = exp(- (bw/2)²/(2σ²)) = 1/sqrt(2)`
    which leads to
        `bw² = 2 log(2) * 2σ²`.
    The denominator `den = 2σ²` of the Gaussian is thus equal to
        `den = 2σ² = bw² / (2 log(2))`."""
    @inbounds log2_length = spec.log2_size
    N = 1 << log2_length
    halfN = N >> 1
    center = N * T(ψmeta.centerfrequency)
    bw = N * T(ψmeta.bandwidth)
    den = @fastmath bw * bw / T(2.0 * log(2.0))
    """2. **Angle**"""
    n_orientations = get_n_orientations(spec.pointgroup)
    θ = T(2π * ψmeta.θ / n_orientations)
    """3. **Number of periods**"""
    halfsupport = sqrt(den * log(inv(spec.ɛ)))
    firstω = max(center - halfsupport, -5N/2)
    lastω = min(center + halfsupport, 5N/2 - 1)
    nPeriods = 1 + ceil(Int, (lastω-halfN) / N)
    #"""3. **Call to morlet**"""
    aspectratio = T(ψmeta.aspectratio)
    y = morlet(FourierDomain(1), θ, aspectratio, center, den, N, nPeriods)
    #"""4. **Trimming to true support boundaries**"""
    #return AbstractFilter(y, spec)
end

function morlet{T<:FFTW.fftwReal}(::FourierDomain{2},
        θ::T, aspectratio::T, center::T, den::T, N::Int, nPeriods::Int)
    halfN = N >> 1
    pstart = - ((nPeriods-1)>>1)
    pstop = (nPeriods-1)>>1 + iseven(nPeriods)
    ωstart = - halfN + pstart * N
    ωstop = halfN + pstop * N - 1
    ωrange = ωstart:ωstop
    nΩs = length(ωrange)
    ωrows = repmat(-ωrange, 1, nΩs)
    ωcols = repmat(ωrange', nΩs, 1)
    cosθ = cos(θ)
    sinθ = sin(θ)
    ωrows_θ = ωrows * cosθ + ωcols * sinθ
    ωcols_θ = -ωrows * sinθ + ωrows * cosθ
end


text:
# test Odo, Obsv2

using IncrementalInference, KernelDensityEstimate
using Base.Test

fg = emptyFactorGraph()

N=100

doors = zeros(1,1)
pd = kde!(doors,[3.0])
pd = resample(pd,N);
bws = getBW(pd)[:,1]
doors2 = getPoints(pd);
v1 = addNode!(fg,:x1, ContinuousScalar,N=N)
f1  = addFactor!(fg,[v1],Obsv2( doors2, bws', [1.0])) #, samplefnc=getSample

tem = 2.0*randn(1,N)+getVal(v1)+50.0
v2 = addNode!(fg,:x2, ContinuousScalar, N=N)
odoc = Odo(50.0*ones(1,1),2.0*ones(1,1),[1.0])
f2 = addFactor!(fg, [:x1; :x2], odoc ) #, samplefnc=getSample

@testset "test evaluation of pose pose constraint" begin
  pts = evalFactor2(fg, f2, v2.index)
  @test norm(Base.mean(pts,2)-[50.0]) < 15.0

  ensureAllInitialized!(fg)
  tree = wipeBuildNewTree!(fg, drawpdf=false)
  inferOverTree!(fg, tree)

  @test norm(Base.mean(getVal(fg, :x2),2)-[50.0]) < 15.0
end

# using RoMEPlotting
#
# plotKDE( getVertKDE(fg,:x2) )
# plotKDE(kde!(pts))

#


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

using JuliaFEM, Test

# 1d strain

X = Dict(1 => [0.0, 0.0, 0.0], 2 => [1.0, 1.0, 1.0])
u = Dict(1 => [0.0, 0.0, 0.0], 2 => [1.0, 1.0, 1.0])
element = Element(Seg2, (1, 2))
update!(element, "geometry", X)
update!(element, "displacement", u)
xi, time = (0.0,), 0.0
detJ = element(xi, time, Val{:detJ})
J = element(xi, time, Val{:Jacobian})
# gradu = element("displacement", xi, time, Val{:Grad})
@debug("1d seg2 info", xi ,time, detJ, J)
@test isapprox(detJ, sqrt(3)/2)
@test isapprox(J, [0.5 0.5 0.5])


text:
using Pretalx
using Test

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


text:
using BinDeps

function getdata(dir::AbstractString,name::AbstractString, url::AbstractString)
    mkpath(dir)
    println(url)
    path = download(url)
    run(unpack_cmd(path,dir,".gz",".tar"))
    dirnext=readdir(dir)[1]
    dirnext=joinpath(dir,dirnext)
    if(isdir(dirnext))
        dir=dirnext
    end
    dictall = Dict()
    cd(joinpath(Pkg.dir("DataDeps"),"datasets"))
    try
        f=open("downloads.json", "r")
        dicttxt = readall(f)
        close(f)
        dictall=JSON.parse(dicttxt)
    end
    dictall[name]=dir
    stringdata = JSON.json(dictall)
    f=open("downloads.json", "w")
    write(f, stringdata)
    close(f)
end


text:
primitive type DiseaseState 8 end

DiseaseState(x::UInt8) = reinterpret(DiseaseState, x)
UInt8(x::DiseaseState) = reinterpret(UInt8, x)

# `DieaseState` const representations
const State_S = DiseaseState(0b0001)
const State_E = DiseaseState(0b0010)
const State_I = DiseaseState(0b0100)
const State_R = DiseaseState(0b1000)

const _DiseaseStateVector = [State_S; State_E; State_I; State_R]
const _DiseaseStateCharVector = ['S'; 'E'; 'I'; 'R']

# `DiseaseState` conversion to/from `Char`
function Base.convert(::Type{Char}, x::DiseaseState)
  return _DiseaseStateCharVector[_DiseaseStateVector .== Ref(x)][1]
end

function Base.convert(::Type{DiseaseState}, x::Char)
  return _DiseaseStateVector[_DiseaseStateCharVector .== Ref(x)][1]
end

# `DiseaseState` conversion to/from `BitArray`
function Base.convert(::Type{BitArray{1}}, x::DiseaseState)
  return BitArray(_DiseaseStateVector .== Ref(x))
end

function Base.convert(::Type{DiseaseState}, x::BitArray{1})
  return _DiseaseStateVector[x]
end

# `DiseaseState` conversion to/from `Int64`
function Base.convert(::Type{Int64}, x::DiseaseState)
  return trailing_zeros(UInt8(x)) + 1
end

function Base.convert(::Type{DiseaseState}, x::Int64)
  return _DiseaseStateVector[x]
end

function Base.show(io::IO, x::DiseaseState)
  return print(io, convert(Char, x))
end

function Base.length(x::DiseaseState)
  return 1
end

const _state_progressions = Dict{DataType, Vector{DiseaseState}}()
_state_progressions[SEIR] = [State_S; State_E; State_I; State_R]
_state_progressions[SEI] = [State_S; State_E; State_I]
_state_progressions[SIR] = [State_S; State_I; State_R]
_state_progressions[SI] = [State_S; State_I]

function advance(x::DiseaseState, ::Type{T}) where T <: EpidemicModel
  current_index = findfirst(Ref(x) .== _state_progressions[T])
  return _state_progressions[T][current_index + 1]
end

function regress(x::DiseaseState, ::Type{T}) where T <: EpidemicModel
  current_index = findfirst(Ref(x) .== _state_progressions[T])
  return _state_progressions[T][current_index - 1]
end

function advance!(x::DiseaseState, ::Type{T}) where T <: EpidemicModel
  x = advance(x, T)
end

function regress!(x::DiseaseState, ::Type{T}) where T <: EpidemicModel
  x = regress(x, T)
end


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl
=#
# TODO add hess_prod, generalize for complex

"""
$(TYPEDEF)

Epigraph of matrix relative entropy cone of dimension `dim` in svec format.

    $(FUNCTIONNAME){T}(dim::Int, use_dual::Bool = false)
"""
mutable struct EpiTrRelEntropyTri{T <: Real} <: Cone{T}
    use_dual_barrier::Bool
    dim::Int
    d::Int

    point::Vector{T}
    dual_point::Vector{T}
    grad::Vector{T}
    dder3::Vector{T}
    vec1::Vector{T}
    vec2::Vector{T}
    feas_updated::Bool
    grad_updated::Bool
    hess_updated::Bool
    inv_hess_updated::Bool
    hess_fact_updated::Bool
    dder3_aux_updated::Bool
    is_feas::Bool
    hess::Symmetric{T, Matrix{T}}
    inv_hess::Symmetric{T, Matrix{T}}
    hess_fact_mat::Symmetric{T, Matrix{T}}
    hess_fact::Factorization{T}

    rt2::T
    vw_dim::Int
    V_idxs::UnitRange{Int}
    W_idxs::UnitRange{Int}
    V::Matrix{T}
    W::Matrix{T}
    V_fact::Eigen{T}
    W_fact::Eigen{T}
    Vi::Matrix{T}
    Wi::Matrix{T}
    W_sim::Matrix{T}
    V_λ_log::Vector{T}
    W_λ_log::Vector{T}
    V_log::Matrix{T}
    W_log::Matrix{T}
    WV_log::Matrix{T}
    z::T
    Δ2_V::Matrix{T}
    Δ2_W::Matrix{T}
    Δ3_V::Array{T, 3}
    Δ3_W::Array{T, 3}
    dzdV::Vector{T}
    dzdW::Vector{T}
    d2zdV2::Matrix{T}
    d2zdW2::Matrix{T}
    d2zdVW::Matrix{T}
    mat::Matrix{T}
    mat2::Matrix{T}
    mat3::Matrix{T}
    mat4::Matrix{T}
    ten3d::Array{T, 3}
    matd2::Matrix{T}

    function EpiTrRelEntropyTri{T}(
        dim::Int;
        use_dual::Bool = false,
        ) where {T <: Real}
        @assert dim > 2
        @assert isodd(dim)
        cone = new{T}()
        cone.use_dual_barrier = use_dual
        cone.dim = dim
        cone.vw_dim = div(dim - 1, 2)
        cone.d = svec_side(cone.vw_dim)
        return cone
    end
end

reset_data(cone::EpiTrRelEntropyTri) = (cone.feas_updated = cone.grad_updated =
    cone.hess_updated = cone.inv_hess_updated = cone.hess_fact_updated =
    cone.dder3_aux_updated = false)

function setup_extra_data!(cone::EpiTrRelEntropyTri{T}) where {T <: Real}
    vw_dim = cone.vw_dim
    d = cone.d
    cone.rt2 = sqrt(T(2))
    cone.V_idxs = 2:(vw_dim + 1)
    cone.W_idxs = (vw_dim + 2):cone.dim
    cone.V = zeros(T, d, d)
    cone.W = zeros(T, d, d)
    cone.Vi = zeros(T, d, d)
    cone.Wi = zeros(T, d, d)
    cone.W_sim = zeros(T, d, d)
    cone.V_λ_log = zeros(T, d)
    cone.W_λ_log = zeros(T, d)
    cone.V_log = zeros(T, d, d)
    cone.W_log = zeros(T, d, d)
    cone.WV_log = zeros(T, d, d)
    cone.Δ2_V = zeros(T, d, d)
    cone.Δ2_W = zeros(T, d, d)
    cone.Δ3_V = zeros(T, d, d, d)
    cone.Δ3_W = zeros(T, d, d, d)
    cone.dzdV = zeros(T, vw_dim)
    cone.dzdW = zeros(T, vw_dim)
    cone.d2zdV2 = zeros(T, vw_dim, vw_dim)
    cone.d2zdW2 = zeros(T, vw_dim, vw_dim)
    cone.d2zdVW = zeros(T, vw_dim, vw_dim)
    cone.mat = zeros(T, d, d)
    cone.mat2 = zeros(T, d, d)
    cone.mat3 = zeros(T, d, d)
    cone.mat4 = zeros(T, d, d)
    cone.ten3d = zeros(T, d, d, d)
    cone.matd2 = zeros(T, d^2, d^2)
    return
end

get_nu(cone::EpiTrRelEntropyTri) = 2 * cone.d + 1

function set_initial_point!(
    arr::AbstractVector{T},
    cone::EpiTrRelEntropyTri{T},
    ) where {T <: Real}
    arr .= 0
    # at the initial point V and W are diagonal, equivalent to epirelentropy
    (arr[1], v, w) = get_central_ray_epirelentropy(cone.d)
    k = 1
    for i in 1:cone.d
        arr[1 + k] = v
        arr[cone.vw_dim + 1 + k] = w
        k += i + 1
    end
    return arr
end

function update_feas(cone::EpiTrRelEntropyTri{T}) where {T <: Real}
    @assert !cone.feas_updated
    point = cone.point

    cone.is_feas = false
    for (X, idxs) in zip((cone.V, cone.W), (cone.V_idxs, cone.W_idxs))
        @views svec_to_smat!(X, point[idxs], cone.rt2)
    end
    VH = Hermitian(cone.V, :U)
    WH = Hermitian(cone.W, :U)
    if isposdef(VH) && isposdef(WH)
        # TODO use LAPACK syev! instead of syevr! for efficiency
        V_fact = cone.V_fact = eigen(VH)
        W_fact = cone.W_fact = eigen(WH)
        if isposdef(V_fact) && isposdef(W_fact)
            for (fact, λ_log, X_log) in zip((V_fact, W_fact),
                (cone.V_λ_log, cone.W_λ_log), (cone.V_log, cone.W_log))
                (λ, vecs) = fact
                @. λ_log = log(λ)
                spectral_outer!(X_log, vecs, λ_log, cone.mat)
            end
            @. cone.WV_log = cone.W_log - cone.V_log
            cone.z = point[1] - dot(WH, Hermitian(cone.WV_log, :U))
            cone.is_feas = (cone.z > 0)
        end
    end

    cone.feas_updated = true
    return cone.is_feas
end

function update_grad(cone::EpiTrRelEntropyTri{T}) where {T <: Real}
    @assert cone.is_feas
    rt2 = cone.rt2
    zi = inv(cone.z)
    (V_λ, V_vecs) = cone.V_fact
    (W_λ, W_vecs) = cone.W_fact
    Vi = cone.Vi
    Wi = cone.Wi
    W_sim = cone.W_sim
    Δ2_V = cone.Δ2_V
    dzdW = cone.dzdW
    mat = cone.mat
    mat2 = cone.mat2
    mat3 = cone.mat3
    g = cone.grad

    spectral_outer!(Vi, V_vecs, inv.(V_λ), cone.mat)
    spectral_outer!(Wi, W_vecs, inv.(W_λ), cone.mat)

    g[1] = -zi

    @views g_W = g[cone.W_idxs]
    smat_to_svec!(g_W, Wi, rt2)
    copyto!(mat, -zi * I)
    @. mat -= zi * cone.WV_log
    smat_to_svec!(dzdW, mat, rt2)
    axpby!(-1, dzdW, -1, g_W)

    Δ2!(Δ2_V, V_λ, cone.V_λ_log)

    spectral_outer!(W_sim, V_vecs', Symmetric(cone.W, :U), mat)
    @. mat = W_sim * Δ2_V
    spectral_outer!(mat2, V_vecs, Symmetric(mat, :U), mat3)
    @views smat_to_svec!(cone.dzdV, mat2, rt2)

    axpby!(-1, Vi, -zi, mat2)
    @views smat_to_svec!(g[cone.V_idxs], mat2, rt2)

    cone.grad_updated = true
    return cone.grad
end

function update_hess(cone::EpiTrRelEntropyTri{T}) where {T <: Real}
    @assert cone.grad_updated
    isdefined(cone, :hess) || alloc_hess!(cone)
    rt2 = cone.rt2
    V_idxs = cone.V_idxs
    W_idxs = cone.W_idxs
    zi = inv(cone.z)
    (V_λ, V_vecs) = cone.V_fact
    (W_λ, W_vecs) = cone.W_fact
    Δ2_V = cone.Δ2_V
    Δ2_W = cone.Δ2_W
    Δ3_V = cone.Δ3_V
    dzdW = cone.dzdW
    d2zdV2 = cone.d2zdV2
    d2zdW2 = cone.d2zdW2
    d2zdVW = cone.d2zdVW
    mat = cone.mat
    mat2 = cone.mat2
    mat3 = cone.mat3
    mat4 = cone.mat4
    H = cone.hess.data

    # u
    @views dzdVzi = H[V_idxs, 1]
    @. dzdVzi = zi * cone.dzdV

    H[1, 1] = abs2(zi)
    @. H[1, V_idxs] = zi * dzdVzi
    @. H[1, W_idxs] = zi * dzdW

    # vv
    Δ3!(Δ3_V, Δ2_V, V_λ)
    d2zdV2!(d2zdV2, V_vecs, cone.W_sim, Δ3_V, cone.ten3d, cone.matd2,
        mat, mat2, rt2)

    @. d2zdV2 *= -1
    @views Hvv = H[V_idxs, V_idxs]
    symm_kron!(Hvv, cone.Vi, rt2)
    mul!(Hvv, dzdVzi, dzdVzi', true, true)
    @. Hvv += zi * d2zdV2

    # vw
    eig_dot_kron!(d2zdVW, Δ2_V, V_vecs, mat, mat2, mat3, mat4, rt2)
    @views Hvw = H[V_idxs, W_idxs]
    @. Hvw = -zi * d2zdVW
    mul!(Hvw, dzdVzi, dzdW', true, true)

    # ww
    Δ2!(Δ2_W, W_λ, cone.W_λ_log)
    eig_dot_kron!(d2zdW2, Δ2_W, W_vecs, mat, mat2, mat3, mat4, rt2)
    @views Hww = H[W_idxs, W_idxs]
    symm_kron!(Hww, cone.Wi, rt2)
    mul!(Hww, dzdW, dzdW', true, true)
    @. Hww += zi * d2zdW2

    cone.hess_updated = true
    return cone.hess
end

function update_dder3_aux(cone::EpiTrRelEntropyTri)
    @assert !cone.dder3_aux_updated
    @assert cone.hess_updated
    Δ3!(cone.Δ3_W, cone.Δ2_W, cone.W_fact.values)
    cone.dder3_aux_updated = true
    return
end

function dder3(
    cone::EpiTrRelEntropyTri{T},
    dir::AbstractVector{T},
    ) where {T <: Real}
    cone.dder3_aux_updated || update_dder3_aux(cone)
    rt2 = cone.rt2
    V_idxs = cone.V_idxs
    W_idxs = cone.W_idxs
    zi = inv(cone.z)
    (V_λ, V_vecs) = cone.V_fact
    (W_λ, W_vecs) = cone.W_fact
    Δ3_V = cone.Δ3_V
    Δ3_W = cone.Δ3_W
    dzdV = cone.dzdV
    dzdW = cone.dzdW
    mat = cone.mat
    mat2 = cone.mat2
    dder3 = cone.dder3

    u_dir = dir[1]
    @views v_dir = dir[V_idxs]
    @views w_dir = dir[W_idxs]

    # TODO in-place
    Vvd = Symmetric(cone.d2zdV2, :U) * v_dir
    Wwd = Symmetric(cone.d2zdW2, :U) * w_dir
    d2zdVW = Symmetric(cone.d2zdVW, :U)
    VWwd = d2zdVW * w_dir

    const0 = zi * (u_dir + dot(v_dir, dzdV)) + dot(w_dir, dzdW)
    const1 = abs2(const0) + zi * (-dot(v_dir, VWwd) +
        (dot(v_dir, Vvd) + dot(w_dir, Wwd)) / 2)

    V_part_1a = const0 * (Vvd - VWwd)
    W_part_1a = Wwd - d2zdVW * v_dir

    # u
    ziconst1 = dder3[1] = zi * const1

    # v, w
    # TODO prealloc
    V_dir = Symmetric(zero(mat), :U)
    W_dir = Symmetric(zero(mat), :U)
    V_dir_sim = zero(mat)
    W_dir_sim = zero(mat)
    VW_dir_sim = zero(mat)
    W_part_1 = zero(mat)
    V_part_1 = zero(mat)
    d3WlogVdV = zero(mat)
    diff_dot_V_VV = zero(mat)
    diff_dot_V_VW = zero(mat)
    diff_dot_W_WW = zero(mat)

    svec_to_smat!(V_dir.data, v_dir, rt2)
    svec_to_smat!(W_dir.data, w_dir, rt2)
    spectral_outer!(V_dir_sim, V_vecs', V_dir, mat)
    spectral_outer!(W_dir_sim, W_vecs', W_dir, mat)
    spectral_outer!(VW_dir_sim, V_vecs', W_dir, mat)

    @inbounds @views for j in 1:cone.d
        Vds_j = V_dir_sim[:, j]
        Wds_j = W_dir_sim[:, j]
        for i in 1:j
            Vds_i = V_dir_sim[:, i]
            Wds_i = W_dir_sim[:, i]
            DΔ3_V_ij = Diagonal(Δ3_V[:, i, j])
            DΔ3_W_ij = Diagonal(Δ3_W[:, i, j])
            diff_dot_V_VV[i, j] = dot(Vds_j, DΔ3_V_ij, Vds_i)
            diff_dot_W_WW[i, j] = dot(Wds_j, DΔ3_W_ij, Wds_i)
        end
        for i in 1:cone.d
            VWds_i = VW_dir_sim[:, i]
            DΔ3_V_ij = Diagonal(Δ3_V[:, i, j])
            diff_dot_V_VW[i, j] = dot(Vds_j, DΔ3_V_ij, VWds_i)
        end
    end

    # v
    d3WlogVdV!(d3WlogVdV, Δ3_V, V_λ, V_dir_sim, cone.W_sim, mat)
    svec_to_smat!(V_part_1, V_part_1a, rt2)
    rdiv!(V_dir_sim, Diagonal(sqrt.(V_λ)))
    ldiv!(Diagonal(V_λ), V_dir_sim)
    V_part_2 = d3WlogVdV
    @. V_part_2 += diff_dot_V_VW + diff_dot_V_VW'
    mul!(V_part_2, V_dir_sim, V_dir_sim', true, zi)
    mul!(mat, Symmetric(V_part_2, :U), V_vecs')
    mul!(V_part_1, V_vecs, mat, true, zi)
    @views dder3_V = dder3[V_idxs]
    smat_to_svec!(dder3_V, V_part_1, rt2)
    @. dder3_V += ziconst1 * dzdV

    # w
    svec_to_smat!(W_part_1, W_part_1a, rt2)
    spectral_outer!(mat2, V_vecs, Symmetric(diff_dot_V_VV, :U), mat)
    axpby!(true, mat2, const0, W_part_1)
    rdiv!(W_dir_sim, Diagonal(sqrt.(W_λ)))
    ldiv!(Diagonal(W_λ), W_dir_sim)
    W_part_2 = diff_dot_W_WW
    mul!(W_part_2, W_dir_sim, W_dir_sim', true, -zi)
    mul!(mat, Symmetric(W_part_2, :U), W_vecs')
    mul!(W_part_1, W_vecs, mat, true, zi)
    @views dder3_W = dder3[W_idxs]
    smat_to_svec!(dder3_W, W_part_1, rt2)
    @. dder3_W += const1 * dzdW

    return dder3
end

function Δ2!(
    Δ2::Matrix{T},
    λ::Vector{T},
    log_λ::Vector{T},
    ) where {T <: Real}
    rteps = sqrt(eps(T))
    d = length(λ)

    @inbounds for j in 1:d
        λ_j = λ[j]
        lλ_j = log_λ[j]
        for i in 1:(j - 1)
            λ_i = λ[i]
            λ_ij = λ_i - λ_j
            if abs(λ_ij) < rteps
                Δ2[i, j] = 2 / (λ_i + λ_j)
            else
                Δ2[i, j] = (log_λ[i] - lλ_j) / λ_ij
            end
        end
        Δ2[j, j] = inv(λ_j)
    end

    # make symmetric
    copytri!(Δ2, 'U')
    return Δ2
end

function Δ3!(
    Δ3::Array{T, 3},
    Δ2::Matrix{T},
    λ::Vector{T},
    ) where {T <: Real}
    @assert issymmetric(Δ2) # must be symmetric (wrapper is less efficient)
    rteps = sqrt(eps(T))
    d = length(λ)

    @inbounds for k in 1:d, j in 1:k, i in 1:j
        λ_j = λ[j]
        λ_k = λ[k]
        λ_jk = λ_j - λ_k
        if abs(λ_jk) < rteps
            λ_i = λ[i]
            λ_ij = λ_i - λ_j
            if abs(λ_ij) < rteps
                t = abs2(3 / (λ_i + λ_j + λ_k)) / -2
            else
                t = (Δ2[i, j] - Δ2[j, k]) / λ_ij
            end
        else
            t = (Δ2[i, j] - Δ2[i, k]) / λ_jk
        end

        Δ3[i, j, k] = Δ3[i, k, j] = Δ3[j, i, k] =
            Δ3[j, k, i] = Δ3[k, i, j] = Δ3[k, j, i] = t
    end

    return Δ3
end

function d2zdV2!(
    d2zdV2::Matrix{T},
    vecs::Matrix{T},
    inner::Matrix{T},
    Δ3::Array{T, 3},
    ten3d::Array{T, 3}, # temp
    matd2::Matrix{T}, # temp
    mat::Matrix{T}, # temp
    mat2::Matrix{T}, # temp
    rt2::T,
    ) where {T <: Real}
    d = size(vecs, 1)

    @inbounds for i in 1:d
        @views ten3d_i = ten3d[i, :, :]
        @views Δ3_i = Δ3[i, :, :]
        @. mat2 = inner * Δ3_i
        spectral_outer!(ten3d_i, vecs, Symmetric(mat2, :U), mat)
    end
    @inbounds for j in 1:d, i in 1:j
        @views temp_ij = matd2[block_idxs(d, i), block_idxs(d, j)]
        @views D_ij = ten3d[:, i, j]
        spectral_outer!(temp_ij, vecs, D_ij, mat)
    end
    copytri!(matd2, 'U')

    rt2i = inv(rt2)
    row_idx = 1
    @inbounds for j in 1:d, i in 1:j
        col_idx = 1
        ijeq = (i == j)
        di = d * (i - 1)
        dj = d * (j - 1)
        for l in 1:d, k in 1:l
            rho = (k == l ? (ijeq ? T(0.5) : rt2i) : (ijeq ? rt2i : one(T)))
            dlk = d * (l - 1) + k
            dkl = d * (k - 1) + l
            dji = dj + i
            dij = di + j
            d2zdV2[row_idx, col_idx] = rho * (
                matd2[dji, dlk] + matd2[dij, dlk] +
                matd2[dji, dkl] + matd2[dij, dkl])
            col_idx += 1
        end
        row_idx += 1
    end

    return d2zdV2
end

function VdWs_element(
    i::Int,
    j::Int,
    k::Int,
    l::Int,
    Vds::Matrix{T},
    Ws::Matrix{T},
    ) where {T <: Real}
    @inbounds begin
        a = Vds[l, j] * Ws[i, k] + Vds[l, i] * Ws[j, k]
        b = Vds[k, i] * Vds[l, j] * Ws[k, l]
        c = Vds[k, l] * a + b
    end
    return c
end

function d3WlogVdV!(
    d3WlogVdV::Matrix{T},
    Δ3::Array{T, 3},
    λ::Vector{T},
    Vds::Matrix{T},
    Ws::Matrix{T},
    Δ4_ij::Matrix{T}, # temp
    ) where {T <: Real}
    d = length(λ)

    @inbounds for j in 1:d, i in 1:j
        Δ4_ij!(Δ4_ij, i, j, Δ3, λ)
        t = zero(T)
        for l in 1:d
            for k in 1:(l - 1)
                t += Δ4_ij[k, l] * (VdWs_element(i, j, k, l, Vds, Ws) +
                    VdWs_element(i, j, l, k, Vds, Ws))
            end
            t += Δ4_ij[l, l] * VdWs_element(i, j, l, l, Vds, Ws)
        end
        d3WlogVdV[i, j] = t
    end

    copytri!(d3WlogVdV, 'U')
    return d3WlogVdV
end

function Δ4_ij!(
    Δ4_ij::Matrix,
    i::Int,
    j::Int,
    Δ3::Array{T, 3},
    λ::Vector{T},
    ) where {T <: Real}
    rteps = sqrt(eps(T))
    d = length(λ)
    λ_i = λ[i]
    λ_j = λ[j]

    @inbounds for l in 1:d, k in 1:l
        λ_k = λ[k]
        λ_l = λ[l]
        λ_ij = λ_i - λ_j
        λ_ik = λ_i - λ_k
        λ_il = λ_i - λ_l
        B_ik = (abs(λ_ik) < rteps)
        B_il = (abs(λ_il) < rteps)

        if (abs(λ_ij) < rteps) && B_ik && B_il
            t = λ_i^-3 / 3
        elseif B_ik && B_il
            t = (Δ3[i, i, i] - Δ3[i, i, j]) / λ_ij
        elseif B_il
            t = (Δ3[i, i, j] - Δ3[i, j, k]) / λ_ik
        else
            t = (Δ3[i, j, k] - Δ3[j, k, l]) / λ_il
        end

        Δ4_ij[k, l] = t
    end

    return Δ4_ij
end


text:
__precompile__(true)

module Somoclu

using BinDeps
using MultivariateStats: PCA, fit, principalvars, projection
using Base: AsyncCondition

export Som, train!, distance

if isfile(joinpath(dirname(@__FILE__),"..","deps","deps.jl"))
    include("../deps/deps.jl")
else
    error("Somoclu not properly installed. Please run Pkg.build(\"Somoclu\")")
end

"""
    Som(ncolumns, nrows; <keyword arguments>)

Self-organizing map of size `ncolumns`x`nrows`.

# Arguments
* `compactsupport::Bool=true`: Cut off map updates beyond the training radius
                                with the Gaussian neighborhood.
* `gridtype::String="rectangular"`: Specify the grid form of the nodes:
                                    `"rectangular"` or `"hexagonal"`
* `initialcodebook::Array{32, 2}=nothing`: Specify an initial codebook.
* `initialization::String="random"`: Specify the codebook initialization:
                                     `"random"` or `"pca"`.
* `kerneltype::Integer=0`: Specify which kernel to use: 0 for dense CPU kernel
                          or 1 for dense GPU kernel (if compiled with it).
* `maptype::String="planar"`: Specify the map topology: `"planar"` or `"toroid"`
* `neighborhood::String="gaussian"`: Specify the neighborhood function:
                                     `"gaussian"` or `"bubble"`.
* `stdcoeff::Float32=0.5`: Coefficient in the Gaussian neighborhood function
                           exp(-||x-y||^2/(2\*(coeff\*radius)^2))
* `verbose::Integer=0p`: Specify verbosity: 0, 1, or 2.
"""
mutable struct Som
    ncolumns::Int
    nrows::Int
    kerneltype::Int
    maptype::String
    gridtype::String
    compactsupport::Bool
    neighborhood::String
    stdcoeff::Float32
    verbose::Int
    useCustomDistance::Bool
    initialization::String
    codebook::Array{Float32, 2}
    bmus::Array{Cint, 2}
    umatrix::Array{Float32, 2}

    function Som(ncolumns, nrows; kerneltype=0, maptype="planar",
        gridtype="square", compactsupport=true, neighborhood="gaussian",
        stdcoeff=0.5, verbose=0, useCustomDistance=false, initialization="random",
        initialcodebook=nothing)

        maptype != "planar"        && maptype != "toroid" &&
            error("Unknown map type!")
        gridtype != "square"       && gridtype != "hexagonal" &&
            error("Unknown grid type!")
        neighborhood != "gaussian" && neighborhood != "bubble" &&
            error("Unknown neighborhood function!")
        initialization != "random" && initialization != "pca" &&
            error("Unknown initialization!")
        kerneltype != 0            && kerneltype != 1 &&
            error("Unsupported kernel type!")
        verbose < 0                && verbose > 2 &&
            error("Unsupported verbosity level!")

        if initialcodebook != nothing
            new(ncolumns, nrows, kerneltype, maptype, gridtype, compactsupport,
            neighborhood, stdcoeff, verbose, useCustomDistance, initialization,
            initialcodebook)
        else
            new(ncolumns, nrows, kerneltype, maptype, gridtype, compactsupport,
            neighborhood, stdcoeff, verbose, useCustomDistance, initialization)
        end
    end
end

distance(p1::Ptr{Cfloat}, p2::Ptr{Cfloat}, d::Cuint) = Cfloat(0.0)::Cfloat

function get_parameters(context::Ptr{Void})
    pr = Ptr{Cfloat}(context)::Ptr{Cfloat}
    dim = unsafe_load(Ptr{Cuint}(context + Core.sizeof(Cfloat)))::Cuint
    p1a = Ptr{Ptr{Cfloat}}(pr + Core.sizeof(Cfloat) + Core.sizeof(Cuint))::Ptr{Ptr{Cfloat}}
    p2a = (p1a + Core.sizeof(Ptr{Cfloat}))::Ptr{Ptr{Cfloat}}
    p1 = unsafe_load(p1a)
    p2 = unsafe_load(p2a)
    rem = Ptr{Void}(p2a + Core.sizeof(Ptr{Cfloat}))::Ptr{Void}
    return dim, p1, p2, pr, rem
end

function nothreads_distance(context::Ptr{Void})
    dim, p1, p2, pr, ignore = get_parameters(context)
    r = distance(p1, p2, dim)
    unsafe_store!(pr, r)
    return r::Cfloat
end

function __init__()
    global const fdist_nt_c = cfunction(nothreads_distance, Cfloat, (Ptr{Void}, ))
    global mt_data_map = Dict{AsyncCondition, Ref{Ptr{Void}}}()
end

"""
    train!(som::Som, data::Array{Float32, 2}; <keyword arguments>)

Train a self-organizing map `som` on `data`.

The som will be updated during training.

# Arguments
* `epochs::Integer=10`: The number of epochs to train the map for.
* `radius0::Float32=0`: The initial radius on the map where the update happens
                       around a best matching unit. Default value of 0 will
                       trigger a value of min(n_columns, n_rows)/2.
* `radiusN::Float32=1`: The radius on the map where the update happens around a
                       best matching unit in the final epoch.
* `radiuscooling::String="linear"`: The cooling strategy between radius0 and
                                    radiusN: `"linear"` or `"exponential"`.
* `scale0::Float32=0.1`: The initial learning scale.
* `scaleN::Float32=0.01`: The learning scale in the final epoch.
* `scalecooling::String="linear"`: The cooling strategy between scale0 and
                                   scaleN: `"linear"` or `"exponential"`.

"""
function train!(som::Som, data::Array{Float32, 2}; epochs=10, radius0=0, radiusN=1, radiuscooling="linear", scale0=0.1, scaleN=0.01, scalecooling="linear")
    nDimensions, nVectors = size(data)
    if som.initialization == "random"
        som.codebook = Array{Float32}(nDimensions, som.ncolumns*som.nrows);
        # These two lines trigger the C++ code to randomly initialize the codebook
        som.codebook[1, 1] = 1000.0
        som.codebook[2, 1] = 2000.0
    elseif som.initialization == "pca"
        coord = zeros(Float32, som.ncolumns*som.nrows, 2);
        for i = 1:som.ncolumns*som.nrows
            coord[i, 1] = div(i-1, som.ncolumns)
            coord[i, 2] = rem(i-1, som.ncolumns)
        end
        coord = coord ./ [som.nrows-1 som.ncolumns-1];
        coord = 2*(coord - .5);
        me = mean(data, 2);
        M = fit(PCA, data.-me; maxoutdim=2);
        eigval = principalvars(M);
        eigvec = projection(M);
        norms = [norm(eigvec[:, i]) for i in 1:2];
        eigvec = ((eigvec' ./ norms) .* eigval)'
        som.codebook = repeat(me, outer=[1, som.ncolumns*som.nrows])
        for j = 1:som.ncolumns*som.nrows
            for i = 1:2
                som.codebook[:, j] = som.codebook[:, j] + coord[j, i] * eigvec[:, i]
            end
        end
    else
        error("Unknown initialization method")
    end
    _radiuscooling = radiuscooling == "linear"      ? 0 :
                     radiuscooling == "exponential" ? 1 :
                     error("Unknown radius cooling")

    _scalecooling  = scalecooling == "linear"       ? 0 :
                     scalecooling == "exponential"  ? 1 :
                     error("Unknown scale cooling")

    _maptype       = som.maptype == "planar"        ? 0 :
                     som.maptype == "toroid"        ? 1 :
                     error("Unknown map type!")

    _gridtype      = som.gridtype == "square"       ? 0 :
                     som.gridtype == "hexagonal"    ? 1 :
                     error("Unknown grid type!")

    bmus = Array{Cint}(nVectors*2);
    umatrix = Array{Float32}(som.ncolumns*som.nrows);

    global fdist_nt_c

    fp = !som.useCustomDistance ? C_NULL : fdist_nt_c

    # Note that som.ncolumns and som.nrows are swapped because Julia is column-first
    ccall((:julia_train, libsomoclu), Void,
          (Ptr{Float32}, Cint, Cuint, Cuint, Cuint, Cuint, Cuint, Float32, Float32, Cuint,
           Float32, Float32, Cuint, Cuint, Cuint, Cuint, Bool, Bool, Float32, Cuint,
           Ptr{Float32}, Cint, Ptr{Cint}, Cint, Ptr{Float32}, Cint, Ptr{Void}),
          reshape(data, length(data)), length(data), epochs, som.nrows, som.ncolumns,
          nDimensions, nVectors, radius0, radiusN, _radiuscooling, scale0, scaleN,
          _scalecooling, som.kerneltype, _maptype, _gridtype, som.compactsupport,
          som.neighborhood=="gaussian", som.stdcoeff, som.verbose, reshape(som.codebook,
          length(som.codebook)), length(som.codebook), bmus, length(bmus), umatrix,
          length(umatrix), fp)

    som.umatrix = reshape(umatrix, som.nrows, som.ncolumns);
    som.bmus = reshape(bmus, 2, nVectors)
    som.bmus[1, :], som.bmus[2, :] = som.bmus[2, :]+1, som.bmus[1, :]+1;
    return nothing
end

end


text:
using RCall
using JellyMe4

using DataFrames
using MixedModels
using Test

const datasets = MixedModels.datasets
dataset(x) = DataFrame(MixedModels.dataset(x))

# this should only occur on CIs
# the Linux CI installs lme4 via apt beforehand
# and there are binary builds for Mac and Windows
reval("""
if(!require(lme4)){
    # use tmp for tests if lme4 isn't available
    .libPaths("/tmp")
    lib <- .libPaths()[1L]
    warning("installing lme4")
    warning(lib)
    install.packages("lme4",repos="https://cloud.r-project.org", libs=lib, type="binary")
    library(lme4)
}
if(!require(afex)){
    .libPaths("/tmp")
    lib <- .libPaths()[1L]
    install.packages("afex",repos="https://cloud.r-project.org", libs=lib, , type="binary")
}
""")

include("merMod.jl")
include("lmerMod.jl")
include("glmerMod.jl")


text:
struct PiecewiseUniform <: Distribution{Float64} end

"""
    piecewise_uniform(bounds, probs)

Samples a `Float64` value from a piecewise uniform continuous distribution.

There are `n` bins where `n = length(probs)` and `n + 1 = length(bounds)`.
Bounds must satisfy `bounds[i] < bounds[i+1]` for all `i`.
The probability density at `x` is zero if `x <= bounds[1]` or `x >= bounds[end]` and is otherwise `probs[bin] / (bounds[bin] - bounds[bin+1])` where `bounds[bin] < x <= bounds[bin+1]`.
"""
const piecewise_uniform = PiecewiseUniform()

function check_dims(::PiecewiseUniform, bounds, probs)
    if length(bounds) != length(probs) + 1
        error("Dimension mismatch")
    end
end

function get_bin(bounds, x)
    @assert x <= bounds[end]
    bin = 1
    while x > bounds[bin+1]
        bin += 1
    end
    @assert x > bounds[bin] && x <= bounds[bin+1]
    bin
end

function logpdf(::PiecewiseUniform, x::Real, bounds::AbstractVector{T},
                    probs::AbstractVector{U}) where {T <: Real, U <: Real}
    check_dims(piecewise_uniform, bounds, probs)

    # bounds[1]      bounds[2]           bounds[3]      bounds[4]
    # ^              ^                   ^              ^
    # |    probs[1]  |  probs[2]         | probs[3]     |
    if x <= bounds[1] || x >= bounds[end]
        -Inf
    else
        bin = get_bin(bounds, x)
        log(probs[bin]) - log(bounds[bin+1] - bounds[bin])
    end
end

function random(::PiecewiseUniform, bounds::Vector{T},
                    probs::Vector{U}) where {T <: Real, U <: Real}
    bin = categorical(probs)
    uniform_continuous(bounds[bin], bounds[bin+1])
end

(::PiecewiseUniform)(bounds, probs) = random(PiecewiseUniform(), bounds, probs)

function logpdf_grad(::PiecewiseUniform, x::Real, bounds, probs)
    check_dims(piecewise_uniform, bounds, probs)
    if x <= bounds[1] || x >= bounds[end]
        error("Out of bounds")
    end
    bin = get_bin(bounds, x)
    bounds_grad = fill(0., length(bounds))
    bin_length = bounds[bin+1] - bounds[bin]
    bounds_grad[bin] = 1. / bin_length
    bounds_grad[bin+1] = - 1. / bin_length
    probs_grad = fill(0., length(probs))
    probs_grad[bin] = 1. / probs[bin]
    (0., bounds_grad, probs_grad)
end

is_discrete(::PiecewiseUniform) = false

has_output_grad(::PiecewiseUniform) = true
has_argument_grads(::PiecewiseUniform) = (true, true)

export piecewise_uniform


text:
## ExoplanetsSysSim/src/io.jl
## (c) 2015 Eric B. Ford

module SysSimIO
using ExoplanetsSysSim
using HDF5, JLD

#if VERSION >= v"0.5-"
#  import Compat: UTF8String, ASCIIString
#end

export save_sim_param, save_sim_results, load_sim_param, load_distances, load_summary_stats

function save_sim_param(filename::String, p::SimParam)
 local file
 try
   file = JLD.jldopen(filename,"w")
   write_sim_param(file,p)
 catch
   println("# Problem writing parameters to jld file: ", filename)
   return false
 finally
   close(file)
 end
 return true
end

function save_sim_results(filename::String, p::SimParam; distances::Vector{Float64}=Array{Float64}(0), summary_stats::CatalogSummaryStatistics = CatalogSummaryStatistics() )
 local file
 try
   file = JLD.jldopen(filename,"w")
   write_sim_param(file,p)
   write_sim_summary_stats(file,summary_stats)
   write_sim_distances(file,distances)
 catch
   println("# Problem writing data to jld file: ", filename)
 finally
   close(file)
  end
end

function write_sim_distances(file::JLD.JldFile, d::Vector{Float64} )
  JLD.write(file,"distances",d)
end

function write_sim_summary_stats(file::JLD.JldFile, ss::CatalogSummaryStatistics )
  JLD.write(file,"summary_stats",ss.stat)
end

function write_sim_param(file::JLD.JldFile, p::SimParam)
 sim_param_bool = Dict{String,Bool}()
 sim_param_int = Dict{String,Integer}()
 sim_param_real = Dict{String,Real}()
 sim_param_function = Dict{String,String}()
 sim_param_string = Dict{String,String}()
 for k in keys(p.param)
   #println("# k=",k,".  v=",p.param[k])
   if typeof(p.param[k]) <: Bool
      sim_param_bool[k] = p.param[k]
   elseif typeof(p.param[k]) <: Integer
      sim_param_int[k] = p.param[k]
   elseif typeof(p.param[k]) <: Real
      sim_param_real[k] = p.param[k]
   elseif typeof(p.param[k]) <: Function
      sim_param_function[k] = string(p.param[k])
   elseif typeof(p.param[k]) <: AbstractString
      sim_param_string[k] = convert(String,p.param[k])
   else
	  warn(string("Can't store value of >",k,"< due to type ", typeof(p.param[k])))
   end
 end
   JLD.write(file,"sim_param_int",sim_param_int)
   JLD.write(file,"sim_param_real",sim_param_real)
   JLD.write(file,"sim_param_function",sim_param_function)
   JLD.write(file,"sim_param_string",sim_param_string)
   JLD.write(file,"sim_param_bool",sim_param_bool)
   JLD.write(file,"sim_param_active",p.active)
end

function load_sim_param(filename::String)
  local jld_data
  try
    jld_data = load(filename)
  catch
   println("# Problem reading parameters from jld file: ", filename)
  end
  p = SimParam()
  merge!(p.active, jld_data["sim_param_active"])
  merge!(p.param, jld_data["sim_param_int"])
  merge!(p.param, jld_data["sim_param_string"])
  merge!(p.param, jld_data["sim_param_real"])
  #=
  merge!(p.param, jld_data["sim_param_function"])
  df = jld_data["sim_param_function"]
  for k in keys(df)
    p.param[k]::Function = symbol(df[k])
  end
  =#
  merge!(p.param, jld_data["sim_param_bool"])
  return p
end

function load_distances(filename::String)
  local jld_data
  try
    jld_data = load(filename)
  catch
   println("# Problem reading distances from jld file: ", filename)
  end
  d::Array{Float64,1} = jld_data["distances"]
  return d
end

function load_summary_stats(filename::String)
  local jld_data
  try
    jld_data = load(filename)
  catch
   println("# Problem reading parameters from jld file: ", filename)
  end
  s = CatalogSummaryStatistics()
  merge!(s.stat, jld_data["summary_stats"])
  return s
end

function test_io()
  sim_param = setup_sim_param_demo()
  save_sim_param("test.jld",sim_param)
  spd = load_sim_param("test.jld")
  rm("test.jld")
end

end # module SysSimIO


text:
export gaussian_kernel, sigmoid_kernel, cauchy_kernel, linear_kernel, tanh_kernel

"""
Creates a Gaussian kernel with the given distance function and `sigma` value
"""
function gaussian_kernel(dist, sigma=1.0)
    sigma2 = sigma * 2
    function fun(obj, ref)::Float64
        d = dist(obj, ref)
        (d == 0 || sigma == 0) && return 1.0
        exp(-d / sigma2)
    end

    fun
end

"""
Creates a sigmoid kernel with the given `sigma` value and distance function
"""
function sigmoid_kernel(dist, sigma=1.0)
    sqrtsigma = sqrt(sigma)
    function fun(obj, ref)::Float64
        x = dist(obj, ref)
        2 * sqrtsigma / (1 + exp(-x))
    end

    fun
end

"""
Creates a Cauchy's kernel with the given `sigma` value and distance function
"""
function cauchy_kernel(dist, sigma=1.0)
    sqsigma = sigma^2
    function fun(obj, ref)::Float64
        x = dist(obj, ref)
        (x == 0 || sqsigma == 0) && return 1.0
        1 / (1 + x^2 / sqsigma)
    end

    fun
end

"""
Creates a tanh kernel with the given `sigma` value and distance function
"""
function tanh_kernel(dist, sigma=1.0)
    function fun(obj, ref)::Float64
        x = dist(obj, ref)
        (exp(x-sigma) - exp(-x+sigma)) / (exp(x-sigma) + exp(-x+sigma))
    end

    fun
end

"""
Creates a linear kernel with the given distance function and `sigma` slope
"""
function linear_kernel(dist, sigma=1.0)
    function fun(obj, ref)::Float64
        dist(obj, ref) * sigma
    end

    fun
end