text:
# constants
const ATTR_LAT   = Dict("description" => "Latitude", "unit" => "°");
const ATTR_LON   = Dict("description" => "Longitude", "unit" => "°");
const ATTR_CYC   = Dict("description" => "Cycle index", "unit" => "-");
const ATTR_ABOUT = Dict("about" => "This is a file generated using PkgUtility.jl",
                        "notes" => "PkgUtility.jl uses NCDatasets.jl to create NC files");

# size_nc Wrapper using NCDatasets
"""
    size_nc(file::String, var::String)

Return the dimensions and size of a NetCDF dataset, given
- `file` Dataset path
- `var` Variable name

---
# Examples
```julia
ndims,sizes = read_nc("test.nc", "test");
```
"""
function size_nc(file::String, var::String)
    _dset = Dataset(file, "r");
    _dvar = _dset[var];
    _ndim = ndims(_dvar);
    _size = size(_dvar);
    close(_dset);

    return _ndim, _size
end

# read_nc Wrapper using NCDatasets
"""
NCDatasets.jl and NetCDF.jl both provide function to read data out from NC dataset. However, while NetCDF.jl is more convenient to use (less lines of code to read data), NCDatasets.jl is better to
    read a subset from the dataset and is able to detect the scale factor and offset. Here, we used a wrapper function to read NC dataset using NCDatasets.jl:

$(METHODLIST)

"""
function read_nc end

"""
When only file name and variable label are provided, `read_nc` function reads out all the data:

    read_nc(file::String, var::String)

Read data from NC file, given
- `file` Dataset path
- `var` Variable to read

Note that the missing data will be labeled as NaN.

---
Examples
```julia
# read data labeled as test from test.nc
data = read_nc("test.nc", "test");
```
"""
read_nc(file::String, var::String) = (
    _dset = Dataset(file, "r");
    _dvar = _dset[var][:,:];
    _nvar = replace(_dvar, missing=>NaN);
    _dvar = nothing;
    close(_dset);

    return _nvar
);

"""
If a float type is given, the data will be converted to FT, namely the output will be an array of float numbers:

    read_nc(FT, file::String, var::String)

Read data from nc file, given
- `FT` Float number type
- `file` Dataset path
- `var` Variable name

---
Examples
```julia
# read data labeled as test from test.nc as Float32
data = read_nc(Float32, "test.nc", "test");
```
"""
read_nc(FT, file::String, var::String) = FT.(read_nc(file, var));

"""
In many cases, the NC dataset can be very huge, and reading all the data points into one array could be time and memory consuming. In this case, reading a subset of data would be the best option:

    read_nc(file::String, var::String, indz::Int)

Read a subset from nc file, given
- `file` Dataset path
- `var` Variable name
- `indz` The 3rd index of subset data to read

Note that the dataset must be a 3D array to use this method.

---
Examples
```julia
# read 1st layer data labeled as test from test.nc
data = read_nc("test.nc", "test", 1);
```
"""
read_nc(file::String, var::String, indz::Int) = (
    _dset = Dataset(file, "r");
    _dvar = _dset[var][:,:,indz];
    _data = replace(_dvar, missing=>NaN);
    _dvar = nothing;
    close(_dset);

    return _data
);

"""
Similarly, one may want to read the subset as a certain float type using

    read_nc(FT, file::String, var::String, indz::Int)

Read a subset from nc file, given
- `FT` Float number type
- `file` Dataset path
- `var` Variable name
- `indz` The 3rd index of subset data to read

---
Examples
```julia
# read 1st layer data labeled as test from test.nc as Float32
data = read_nc(Float32, "test.nc", "test", 1);
```
"""
read_nc(FT, file::String, var::String, indz::Int) = FT.(read_nc(file, var, indz));

"""
Another convenient wrapper is to read all the data for given index in x and y, for example, if one wants to read the time series of data at a given site:

    read_nc(file::String, var::String, indx::Int, indy::Int)

Read the time series of data for a site, given
- `file` Dataset path
- `var` Variable name
- `indx` The 1st index of subset data to read, typically longitude
- `indy` The 2nd index of subset data to read, typically latitude

---
Examples
```julia
data = read_nc("test.nc", "test", 1, 1);
```
"""
read_nc(file::String, var::String, indx::Int, indy::Int) = (
    _dset = Dataset(file, "r");
    _dvar = _dset[var][indx,indy,:];
    _data = replace(_dvar, missing=>NaN);
    _dvar = nothing;
    close(_dset);

    return _data
);

"""
Similarly, one may want to read the subset as a certain float type using

    read_nc(FT, file::String, var::String, indx::Int, indy::Int)

Read the time series of data for a site, given
- `FT` Float number type
- `file` Dataset path
- `var` Variable name
- `indx` The 1st index of subset data to read, typically longitude
- `indy` The 2nd index of subset data to read, typically latitude

---
Examples
```julia
data = read_nc(Float32, "test.nc", "test", 1, 1);
```
"""
read_nc(FT, file::String, var::String, indx::Int, indy::Int) = FT.(read_nc(file, var, indx, indy));

"""
Another convenient wrapper is to read the data for given index in x, y, and z, for example, if one wants to read the time series of data at a given site:

    read_nc(file::String, var::String, indx::Int, indy::Int, indz::Int)

Read the time series of data for a site, given
- `file` Dataset path
- `var` Variable name
- `indx` The 1st index of subset data to read, typically longitude
- `indy` The 2nd index of subset data to read, typically latitude
- `indz` The 3rd index of subset data to read, typically time

---
Examples
```julia
data = read_nc("test.nc", "test", 1, 1, 1);
```
"""
read_nc(file::String, var::String, indx::Int, indy::Int, indz::Int) = (
    _dset = Dataset(file, "r");
    _dvar = _dset[var][indx,indy,indz];
    _data = ismissing(_dvar) ? NaN : _dvar;
    _dvar = nothing;
    close(_dset);

    return _data
);

"""
Similarly, one may want to read the data as a certain float type using

    read_nc(FT, file::String, var::String, indx::Int, indy::Int, indz::Int)

Read the time series of data for a site, given
- `FT` Float number type
- `file` Dataset path
- `var` Variable name
- `indx` The 1st index of subset data to read, typically longitude
- `indy` The 2nd index of subset data to read, typically latitude
- `indz` The 3rd index of subset data to read, typically time

---
Examples
```julia
data = read_nc(Float32, "test.nc", "test", 1, 1, 1);
```
"""
read_nc(FT, file::String, var::String, indx::Int, indy::Int, indz::Int) = FT.(read_nc(file, var, indx, indy, indz));

# save Arrays to NC file
"""
NCDatasets.jl does not have a convenient function (1 line command) to save dataset as a file. Thus, we provide a few methods as supplements:

$(METHODLIST)

"""
function save_nc! end

"""
This method is a case if one wants to save both variable and attributes into the target file. This method support saving multiple (N) dimension arrays:

    save_nc!(file::String,
             var_name::String,
             var_attr::Dict,
             var_data::Array{FT,N},
             atts_name::Vector{String},
             atts_attr::Vector{Dict},
             atts_data::Vector{Vector{FT}},
             notes::Dict{String,String}
    ) where {FT<:AbstractFloat,N}

Save dataset as NC file, given
- `file` Path to save the dataset
- `var_name` Variable name for the data in the NC file
- `var_attr` Variable attributes for the data, such as unit and long name
- `var_data` Data to save
- `atts_name` vector of supporting attribute labels, such as `lat` and `lon`
- `atts_attr` Vector of attributes for the supporting attributes, such as unit
- `atts_data` Vector of attributes data, such as the latitude range
- `notes` Global attributes (notes)
- `compress` Compression level fro NetCDF, default is 4

---
Examples
```julia
# generate data to write into NC file
lats = collect(Float64, -85:10:85);
lons = collect(Float64, -175:10:175);
inds = collect(Int, 1:12);
data1 = rand(18) .+ 273.15;
data2 = rand(36,18) .+ 273.15;
data3 = rand(36,18,12) .+ 273.15;

# define the attributes of the dimensions and data
attrn = Dict("description" => "Random temperature", "unit" => "K");
latat = Dict("description" => "Latitude", "unit" => "°");
lonat = Dict("description" => "Longitude", "unit" => "°");
indat = Dict("description" => "Cycle index", "unit" => "-");

# define attributes names, information, and data
atts_name1 = ["lat"];
atts_name2 = ["lon", "lat"];
atts_name3 = ["lon", "lat", "ind"];
atts_attr1 = [latat];
atts_attr2 = [lonat, latat];
atts_attr3 = [lonat, latat, indat];
atts_data1 = Any[lats];
atts_data2 = Any[lons, lats];
atts_data3 = Any[lons, lats, inds];
notes = Dict("description" => "This is a file generated using PkgUtility.jl", "notes" => "PkgUtility.jl uses NCDatasets.jl to create NC files");

# save data as NC files (1D, 2D, and 3D)
save_nc!("data1.nc", "data1", attrn, data1, atts_name1, atts_attr1, atts_data1, notes);
save_nc!("data2.nc", "data2", attrn, data2, atts_name2, atts_attr2, atts_data2, notes);
save_nc!("data3.nc", "data3", attrn, data3, atts_name3, atts_attr3, atts_data3, notes);
```
"""
save_nc!(file::String,
         var_name::String,
         var_attr::Dict{String,String},
         var_data::Array{FT,N},
         atts_name::Vector{String},
         atts_attr::Vector{Dict{String,String}},
         atts_data::Vector,
         notes::Dict{String,String};
         compress::Int = 4
) where {FT<:AbstractFloat,N} = (
    # make sure the data provided match in dimensions and ranges
    @assert length(atts_attr) == length(atts_data) == length(atts_name) == N;
    @assert 0 <= compress <= 9;

    # create a dataset using "c" mode
    _dset = Dataset(file, "c");

    # global title attribute
    for (_title,_notes) in notes
        _dset.attrib[_title] = _notes;
    end;

    # dimensions for each attribute with their own sizes
    for _i in 1:N
        defDim(_dset, atts_name[_i], length(atts_data[_i]));
        _var = defVar(_dset, atts_name[_i], eltype(atts_data[_i]), atts_name[_i:_i]; attrib=atts_attr[_i], deflatelevel=compress);
        _var[:,:] = atts_data[_i];
    end;

    # define variable with attribute units and copy data into it
    _data = defVar(_dset, var_name, FT, atts_name; attrib=var_attr, deflatelevel=compress);
    _data[:,:] = var_data;

    # close dataset file
    close(_dset);

    return nothing
);

"""
To save the code and effort to redefine the common attributes like latitude, longitude, and cycle index, we provide a shortcut method that handles these within the function:

    save_nc!(file::String,
             var_name::String,
             var_attr::Dict{String,String},
             var_data::Array{FT,N};
             notes::Dict{String,String} = DEFAULT_DICT,
             compress::Int = 4
    ) where {FT<:AbstractFloat,N}

Save the 2D or 3D data as NC file, given
- `file` Path to save the dataset
- `var_name` Variable name for the data in the NC file
- `var_attr` Variable attributes for the data, such as unit and long name
- `var_data` Data to save
- `notes` Global attributes (notes)
- `compress` Compression level fro NetCDF, default is 4

---
Examples
```julia
# generate data to write into NC file
data2 = rand(36,18) .+ 273.15;
data3 = rand(36,18,12) .+ 273.15;

# define the attributes and notes
attrn = Dict("description" => "Random temperature", "unit" => "K");
notes = Dict("description" => "This is a file generated using PkgUtility.jl", "notes" => "PkgUtility.jl uses NCDatasets.jl to create NC files");

# save data as NC files (2D and 3D)
save_nc!("data2.nc", "data2", attrn, data2);
save_nc!("data2.nc", "data2", attrn, data2; notes=notes);
save_nc!("data3.nc", "data3", attrn, data3);
save_nc!("data3.nc", "data3", attrn, data3; notes=notes);
```
"""
save_nc!(file::String,
         var_name::String,
         var_attr::Dict{String,String},
         var_data::Array{FT,N};
         notes::Dict{String,String} = ATTR_ABOUT,
         compress::Int = 4
) where {FT<:AbstractFloat,N} = (
    @assert 2 <= N <= 3;
    @assert 0 <= compress <= 9;

    # generate lat and lon information based on the dimensions of the data
    _N_lat = size(var_data, 2);
    _N_lon = size(var_data, 1);
    _res_lat = FT(180) / _N_lat;
    _res_lon = FT(360) / _N_lon;
    _lats = collect(FT, _res_lat/2:_res_lat:180) .- 90;
    _lons = collect(FT, _res_lon/2:_res_lon:360) .- 180;
    if N==3
        _inds = collect(Int,1:size(var_data,3));
    end;

    # define the attributes of the dimensions and data
    if N==2
        _atts_name = ["lon", "lat"];
        _atts_attr = [ATTR_LON, ATTR_LAT];
        _atts_data = Any[_lons, _lats];
    else
        _atts_name = ["lon", "lat", "ind"];
        _atts_attr = [ATTR_LON, ATTR_LAT, ATTR_CYC];
        _atts_data = Any[_lons, _lats, _inds];
    end;

    # save the data to NC file
    save_nc!(file, var_name, var_attr, var_data, _atts_name, _atts_attr, _atts_data, notes; compress=compress);

    return nothing
);

# append Arrays to NC file
"""
NCDatasets.jl does not have a convenient function (1 line command) to append dataset into a file. Thus, we provide a few methods as supplements:

$(METHODLIST)

"""
function append_nc! end

"""
This method append data to an exisiting NC file. If the attributes exist already, then only save the data:

    append_nc!(file::String,
               var_name::String,
               var_attr::Dict{String,String},
               var_data::Array{FT,N},
               atts_name::Vector{String},
               atts_attr::Vector{Dict{String,String}},
               atts_data::Vector;
               compress::Int = 4
    ) where {FT<:AbstractFloat,N}

Append data to existing file, given
- `file` Path to save the dataset
- `var_name` Variable name for the data in the NC file
- `var_attr` Variable attributes for the data, such as unit and long name
- `var_data` Data to save
- `atts_name` vector of supporting attribute labels, such as `lat` and `lon`
- `atts_attr` Vector of attributes for the supporting attributes, such as unit
- `atts_data` Vector of attributes data, such as the latitude range
- `compress` Compression level fro NetCDF, default is 4

---
Examples
```julia
# generate data to write into NC file
lats = collect(Float64, -85:10:85);
lons = collect(Float64, -175:10:175);
inds = collect(Int, 1:12);
data1 = rand(18) .+ 273.15;
data2 = rand(36,18) .+ 273.15;
data3 = rand(36,18,12) .+ 273.15;

# define the attributes of the dimensions and data
attrn = Dict("description" => "Random temperature", "unit" => "K");
latat = Dict("description" => "Latitude", "unit" => "°");
lonat = Dict("description" => "Longitude", "unit" => "°");
indat = Dict("description" => "Cycle index", "unit" => "-");

# define attributes names, information, and data
atts_name1 = ["lat"];
atts_name2 = ["lon", "lat"];
atts_name3 = ["lon", "lat", "ind"];
atts_attr1 = [latat];
atts_attr2 = [lonat, latat];
atts_attr3 = [lonat, latat, indat];
atts_data1 = Any[lats];
atts_data2 = Any[lons, lats];
atts_data3 = Any[lons, lats, inds];
notes = Dict("description" => "This is a file generated using PkgUtility.jl", "notes" => "PkgUtility.jl uses NCDatasets.jl to create NC files");

# save data as NC files (1D, 2D, and 3D)
append_nc!("data1.nc", "datax", attrn, data1, atts_name1, atts_attr1, atts_data1);
append_nc!("data2.nc", "datax", attrn, data2, atts_name2, atts_attr2, atts_data2);
append_nc!("data3.nc", "datax", attrn, data3, atts_name3, atts_attr3, atts_data3);
```
"""
append_nc!(file::String,
           var_name::String,
           var_attr::Dict{String,String},
           var_data::Array{FT,N},
           atts_name::Vector{String},
           atts_attr::Vector{Dict{String,String}},
           atts_data::Vector;
           compress::Int = 4
) where {FT<:AbstractFloat,N} = (
    # make sure the data provided match in dimensions and ranges
    @assert length(atts_attr) == length(atts_data) == length(atts_name) == N;
    @assert 0 <= compress <= 9;

    # read a dataset using "a" mode
    _dset = Dataset(file, "a");

    # dimensions for each attribute with their own sizes if not exisiting
    for _i in 1:N
        if !(atts_name[_i] in listVar(_dset.ncid))
            defDim(_dset, atts_name[_i], length(atts_data[_i]));
            _var = defVar(_dset, atts_name[_i], eltype(atts_data[_i]), atts_name[_i:_i]; attrib=atts_attr[_i], deflatelevel=compress);
            _var[:,:] = atts_data[_i];
        end;
    end;

    # define variable with attribute units and copy data into it
    _data = defVar(_dset, var_name, FT, atts_name; attrib=var_attr, deflatelevel=compress);
    _data[:,:] = var_data;

    # close dataset file
    close(_dset);

    return nothing
);


text:
include("./dedupe.jl")
using Test
using Combinatorics

@testset "Deduplicated vote distributions" begin
    area::Int = 36
    side::Int = 6
    n_voters::Int = 4
    @testset "Rotations" begin
        last_idx = 0
        for row in 1:side
            for col in 1:side
                idx = coord_index(row, col, side)
                @test idx == last_idx + 1  # Indices wrap around
                @test rot(idx, side, 0) == idx
                @test rot(idx, side, 180) == rot(rot(idx, side, 90), side, 90)
                @test rot(idx, side, 270) == rot(rot(idx, side, 90), side, 180)
                @test rot(rot(idx, side, 180), side, 180) == idx
                @test rot(rot(idx, side, 270), side, 90) == idx
                @test_throws DomainError rot(idx, side, 80)
                last_idx = idx
            end
        end
    end

    @testset "Transpose" begin
        for row in 1:side
            for col in 1:side
                idx = coord_index(row, col, side)
                t_row, t_col = coords(transpose(idx, side), side)
                @test t_row == col && t_col == row
                @test transpose(transpose(idx, side), side) == idx
            end
        end
    end

    @testset "Reflections" begin
        for row in 1:side
            for col in 1:side
                idx = coord_index(row, col, side)
                horiz_idx = reflect(idx, side, true)
                vert_idx = reflect(idx, side, false)
                horiz_row, horiz_col = coords(horiz_idx, side)
                vert_row, vert_col = coords(vert_idx, side)
                @test horiz_row == row && horiz_col == side + 1 - col
                @test vert_row == side + 1 - row && vert_col == col
                @test reflect(horiz_idx, side, true) == idx
                @test reflect(vert_idx, side, false) == idx
            end
        end
    end

    @testset "Enumeration of vote distributions" begin
        # We dedupe the vote distributions (that is, we remove vote distributions
        # that differ from one we've seen only by orientation) and then rehydrate
        # the deduped vote distributions by applying 90º, 180º, and 270º rotations,
        # the transpose, horizontal and vertical flips, and the transpose of the 180º
        # rotation.
        all_combos = [c for c in combinations(1:area, n_voters)]
        combos_from_deduped = []
        @test length(all_combos) == binomial(area, n_voters)
        for c in all_combos
            if !combo_duplicated(c, minimum(c), side)
                push!(combos_from_deduped, c)
                for angle in (90, 180, 270)
                    push!(combos_from_deduped, sort([rot(v, side, angle) for v in c]))
                end
                push!(combos_from_deduped, sort([reflect(v, side, true) for v in c]))
                push!(combos_from_deduped, sort([reflect(v, side, false) for v in c]))
                push!(combos_from_deduped, sort([transpose(v, side) for v in c]))
                push!(combos_from_deduped, sort([rot(transpose(v, side), side, 180) for v in c]))
            end
        end
        @test Set(all_combos) == Set(combos_from_deduped)
    end
end

ext: jl


text:
#!/usr/bin/env julia --color=yes

using Documenter, MLToolkit

makedocs(modules=[MLToolkit])

max_stars_repo_name: coreywoodfield/mltoolkit-julia

max_issues_repo_name: coreywoodfield/mltoolkit-julia

max_forks_repo_name: coreywoodfield/mltoolkit-julia


text:
### Methods ###

function CortexModel(config::Dict)
  numCol = get(config, "numCol", 4)
  szL23e = get(config, "szL23e", 8)
  szL4e = get(config, "szL4e", 8)
  szL5e = get(config, "szL5e", 8)
  szL6e = get(config, "szL6e", 8)

  agent = Agent()

  # Create neurons
  for c = 1:numCol
    agent_Col = Agent()

    # Layer 2/3
    addnode!(agent_Col, GenericNeurons(szL23e); name="L23e")

    # Layer 4
    addnode!(agent_Col, GenericNeurons(szL4e); name="L4e")

    # Layer 5
    addnode!(agent_Col, GenericNeurons(szL5e); name="L5e")

    # Layer 6
    addnode!(agent_Col, GenericNeurons(szL6e); name="L6e")

    agent["$c"] = agent_Col
  end

  # Connect synapses
  for c = 1:numCol
    ## Intra-column intra-layer synapses

    # L4e -> L4e
    addnode!(agent, GenericSynapses(szL4e, szL4e); name="L4e_$(c)_L4e_$(c)")
    addedge!(agent, agent["L4e_$(c)_L4e_$(c)"], (
      (agent["$c"]["L4e"], :input),
      (agent["$c"]["L4e"], :output)
    ))

    ## Intra-column inter-layer synapses

    # L4e -> L23e
    addnode!(agent, GenericSynapses(szL4e, szL23e); name="L4e_$(c)_L23e_$(c)")
    addedge!(agent, agent["L4e_$(c)_L23e_$(c)"], (
      (agent["$c"]["L4e"], :input),
      (agent["$c"]["L23e"], :output)
    ))

    # L23e -> L5e
    addnode!(agent, GenericSynapses(szL23e, szL5e); name="L23e_$(c)_L5e_$(c)")
    addedge!(agent, agent["L23e_$(c)_L5e_$(c)"], (
      (agent["$c"]["L23e"], :input),
      (agent["$c"]["L5e"], :output)
    ))

    # L4e -> L6e
    addnode!(agent, GenericSynapses(szL4e, szL6e); name="L4e_$(c)_L6e_$(c)")
    addedge!(agent, agent["L4e_$(c)_L6e_$(c)"], (
      (agent["$c"]["L4e"], :input),
      (agent["$c"]["L6e"], :output)
    ))

    ## Inter-column synapses
    for d = 1:numCol
      # TODO
    end
  end

  return agent
end


text:
immutable logNormal <: ContinuousUnivariateDistribution
    meanlog::Float64
    sdlog::Float64
    function logNormal(ml::Real, sdl::Real)
    	if sdl > 0.0
    		new(float64(ml), float64(sdl))
    	else
    		error("sdlog must be positive")
    	end
	end
end

logNormal(ml::Real) = logNormal(ml, 1.0)
logNormal() = logNormal(0.0, 1.0)

@_jl_dist_2p logNormal lnorm

entropy(d::logNormal) = 1.0 / 2.0 + (1.0 / 2.0) *
                        log(2.0 * pi * d.sdlog^2) + d.meanlog

insupport(d::logNormal, x::Number) = isreal(x) && isfinite(x) && 0 < x

mean(d::logNormal) = exp(d.meanlog + d.sdlog^2 / 2)

function var(d::logNormal)
	sigsq = d.sdlog^2
	return (exp(sigsq) - 1) * exp(2d.meanlog + sigsq)
end


text:
export Tiling, encode

import Base: length, -

"""
    Tiling(ranges::NTuple{N, Tr}) where {N, Tr}

Using a tuple of `ranges` to simulate a tiling.
The length of `ranges` indicates the dimension of tilling.

# Example
```julia
julia> t = Tiling((1:2:5, 10:5:20))
Tiling{2,StepRange{Int64,Int64}}((1:2:5, 10:5:20), [1 3; 2 4])

julia> encode(t, (2, 12))  # encode into an Int
1

julia> encode(t, (2, 18))
3

julia> t2 = t - (1, 3)  # shift a little to get a new Tiling
Tiling{2,StepRange{Int64,Int64}}((0:2:4, 7:5:17), [1 3; 2 4])
```
"""
struct Tiling{N,Tr<:AbstractRange}
    ranges::NTuple{N,Tr}
    inds::LinearIndices{N,NTuple{N,Base.OneTo{Int}}}
    Tiling(ranges::NTuple{N,Tr}) where {N,Tr} =
        new{N,Tr}(ranges, LinearIndices(Tuple(length(r) - 1 for r in ranges)))
end

"""
    (-)(t::Tiling, xs)

Shift `t` along each dimension by each element in `xs`.
"""
function Base.:-(t::Tiling, xs)
    Tiling(Tuple(r .- x for (r, x) in zip(t.ranges, xs)))
end

Base.length(t::Tiling) = reduce(*, (length(r) - 1 for r in t.ranges))

encode(range::AbstractRange, x) = floor(Int, div(x - range[1], step(range)) + 1)

# TODO: use @generator here!
encode(t::Tiling{1}, x::Number) = encode(t.ranges[1], x)
encode(t::Tiling{1}, xs) = encode(t.ranges[1], xs[1])
encode(t::Tiling{2}, xs) =
    t.inds[CartesianIndex(encode(t.ranges[1], xs[1]), encode(t.ranges[2], xs[2]))]
encode(t::Tiling{3}, xs) =
    t.inds[CartesianIndex(
        encode(t.ranges[1], xs[1]),
        encode(t.ranges[2], xs[2]),
        encode(t.ranges[3], xs[3]),
    )]


text:
using Test
using StatsPlots
using MCMCChains

import Logging

unicodeplots()

n_iter = 500
n_name = 3
n_chain = 2

val = randn(n_iter, n_name, n_chain) .+ [1, 2, 3]'
val = hcat(val, rand(1:2, n_iter, 1, n_chain))

chn = Chains(val)

# Silence all warnings.
level = Logging.min_enabled_level(Logging.current_logger())
Logging.disable_logging(Logging.Warn)

@testset "Plotting tests" begin
    # plotting singe plotting types
    println("traceplot")
    display(traceplot(chn, 1))
    println()

    println("meanplot")
    display(meanplot(chn, 1))
    println()

    println("density")
    display(density(chn, 1))
    display(density(chn, 1, append_chains=true))
    println()

    println("autocorplot")
    display(autocorplot(chn, 1))
    println()

    println("ridgelineplot")
    display(ridgelineplot(chn,chn.name_map[:parameters]))
    println()

    println("forestplot")
    display(forestplot(chn,chn.name_map[:parameters]))
    println()

    #ps_contour = plot(chn, :contour)

    println("histogram")
    display(histogram(chn, 1))
    println()

    println("\nmixeddensity")
    display(mixeddensity(chn, 1))

    # plotting combinations
    display(plot(chn))
    display(plot(chn, append_chains=true))
    display(plot(chn, seriestype = (:mixeddensity, :autocorplot)))

    # Test plotting using colordim keyword
    display(plot(chn, colordim = :parameter))

    # Test if plotting a sub-set work.s
    display(plot(chn, 2))
    display(plot(chn, 2, colordim = :parameter))
    println()
end

# Reset log level.
Logging.disable_logging(level)


text:
function addkey!(membernames, nam)
    if !haskey(membernames, nam)
        membernames[nam] = gensym()
    end
    membernames[nam]
end

onearg(e::Expr, f) = e.head == :call && length(e.args) == 2 && e.args[1] == f
onearg(e, f) = false

is_column_expr(x) = false
is_column_expr(e::Expr) = e.head == :$

mapexpr(f, e) = Expr(e.head, map(f, e.args)...)

replace_syms!(x, membernames) = x
replace_syms!(q::QuoteNode, membernames) =
    replace_syms!(Meta.quot(q.value), membernames)
function replace_syms!(e::Expr, membernames)
    if onearg(e, :^)
        e.args[2]
    elseif onearg(e, :_I_)
        @warn "_I_() for escaping variables is deprecated, use cols() instead"
        addkey!(membernames, :($(e.args[2])))
    elseif onearg(e, :cols)
        @warn "cols(x) for escaping variables is deprecated, use \$x instead"
        addkey!(membernames, :($(e.args[2])))
    elseif is_column_expr(e)
        addkey!(membernames, :($(e.args[1])))
    elseif e.head == :quote
        addkey!(membernames, Meta.quot(e.args[1]) )
    elseif e.head == :.
        replace_dotted!(e, membernames)
    else
        mapexpr(x -> replace_syms!(x, membernames), e)
    end
end

is_simple_non_broadcast_call(x) = false
function is_simple_non_broadcast_call(expr::Expr)
    expr.head == :call &&
        length(expr.args) >= 2 &&
        expr.args[1] isa Symbol &&
        all(x -> x isa QuoteNode || onearg(x, :cols) || is_column_expr(x), expr.args[2:end])
end

is_simple_broadcast_call(x) = false
function is_simple_broadcast_call(expr::Expr)
    expr.head == :. &&
        length(expr.args) == 2 &&
        expr.args[1] isa Symbol &&
        expr.args[2] isa Expr &&
        expr.args[2].head == :tuple &&
        all(x -> x isa QuoteNode || onearg(x, :cols) || is_column_expr(x), expr.args[2].args)
end

function args_to_selectors(v)
    t = map(v) do arg
        if arg isa QuoteNode
            arg
        elseif onearg(arg, :cols)
            @warn "cols(x) is deprecated, use \$x instead"
            arg.args[2]
        elseif is_column_expr(arg)
            arg.args[1]
        else
            throw(ArgumentError("This path should not be reached, arg: $(arg)"))
        end
    end

    :(DataFramesMeta.make_source_concrete($(Expr(:vect, t...))))
end

is_macro_head(ex, name) = false
is_macro_head(ex::Expr, name) = ex.head == :macrocall && ex.args[1] == Symbol(name)

const BYROW_SYM = Symbol("@byrow")
const PASSMISSING_SYM = Symbol("@passmissing")
const DEFAULT_FLAGS = (;BYROW_SYM => Ref(false), PASSMISSING_SYM => Ref(false))

extract_macro_flags(ex, exprflags = deepcopy(DEFAULT_FLAGS)) = (ex, exprflags)
function extract_macro_flags(ex::Expr, exprflags = deepcopy(DEFAULT_FLAGS))
    if ex.head == :macrocall
        macroname = ex.args[1]
        if macroname in keys(exprflags)
            exprflag = exprflags[macroname]
            if exprflag[] == true
                throw(ArgumentError("Redundant flag $macroname used."))
            end
            exprflag[] = true
            return extract_macro_flags(MacroTools.unblock(ex.args[3]), exprflags)
        else
            return (ex, exprflags)
        end
    end

    return (ex, exprflags)
end

"""
    check_macro_flags_consistency(exprflags)

Check that the macro flags are consistent with
one another. For now this only checks that
`@passmissing` is only called when `@byrow` is
also called. In the future we may expand
this function or eliminate it all together.
"""
function check_macro_flags_consistency(exprflags)
    if exprflags[PASSMISSING_SYM][]
        if !exprflags[BYROW_SYM][]
            s = "The `@passmissing` flag is currently only allowed with the `@byrow` flag"
            throw(ArgumentError(s))
        end
    end
end

"""
    get_source_fun(function_expr; wrap_byrow::Bool=false)

Given an expression that may contain `QuoteNode`s (`:x`)
and items wrapped in `cols`, return a function
that is equivalent to that expression where the
`QuoteNode`s and `cols` items are the inputs
to the function.

For fast compilation `get_source_fun` returns
the name of a called function where possible.

* `f(:x, :y)` will return `f`
* `f.(:x, :y)` will return `ByRow(f)`
* `:x .+ :y` will return `.+`

`get_source_fun` also returns an expression
representing the vector of inputs that will be
used as the `src` in the `src => fun => dest`
call later on.

If `wrap_byrow=true` then the function gets wrapped
in `ByRow`. If the expression begins with `@byrow`,
then `get_source_fun` is recurively called on the
expression that `@byrow` acts on, with `wrap_byrow=true`.

### Examples

julia> using MacroTools

julia> ex = :(:x + :y);

julia> DataFramesMeta.get_source_fun(ex)
(:(DataFramesMeta.make_source_concrete([:x, :y])), :+)

julia> ex = quote
           :x .+ 1 .* :y
       end |> MacroTools.prettify

julia> src, fun = DataFramesMeta.get_source_fun(ex);

julia> MacroTools.prettify(fun)
:((mammoth, goat)->mammoth .+ 1 .* goat)

julia> ex = :(@byrow :x * :y);

julia> src, fun = DataFramesMeta.get_source_fun(ex);

julia> MacroTools.prettify(fun)
:(ByRow(*))
```

"""
function get_source_fun(function_expr; exprflags = deepcopy(DEFAULT_FLAGS))
    function_expr = MacroTools.unblock(function_expr)

    if is_simple_non_broadcast_call(function_expr)
        source = args_to_selectors(function_expr.args[2:end])
        fun_t = function_expr.args[1]

        # .+ to +
        if startswith(string(fun_t), '.')
            f_sym_without_dot = Symbol(chop(string(fun_t), head = 1, tail = 0))
            fun = :(DataFrames.ByRow($f_sym_without_dot))
        else
            fun = fun_t
        end
    elseif is_simple_broadcast_call(function_expr)
        # extract source symbols from quotenodes
        source = args_to_selectors(function_expr.args[2].args)
        fun_t = function_expr.args[1]
        fun = :(DataFrames.ByRow($fun_t))
    else
        membernames = Dict{Any, Symbol}()

        body = replace_syms!(function_expr, membernames)

        source = :(DataFramesMeta.make_source_concrete($(Expr(:vect, keys(membernames)...))))
        inputargs = Expr(:tuple, values(membernames)...)

        fun = quote
            $inputargs -> begin
                $body
            end
        end
    end

    if exprflags[BYROW_SYM][]
        if exprflags[PASSMISSING_SYM][]
            fun = :(ByRow(Missings.passmissing($fun)))
        else
            fun = :(ByRow($fun))
        end
    end

    return source, fun
end

# `nolhs` needs to be `true` when we have syntax of the form
# `@combine(gd, fun(:x, :y))` where `fun` returns a `table` object.
# We don't create the "new name" pair because new names are
# given by the table.
# We need wrap_byrow as a keyword argument here in case someone
# uses `@transform df @byrow begin ... end`, which we
# deal with outside of this function.
function fun_to_vec(ex::Expr;
                    gensym_names::Bool=false,
                    outer_flags::NamedTuple=deepcopy(DEFAULT_FLAGS),
                    no_dest::Bool=false)
    # classify the type of expression
    # :x # handled via dispatch
    # $:x # handled as though above
    # f(:x) # requires no_dest, for `@with` and `@subset` in future
    # :y = :x # Simple pair
    # :y = $:x # Extract and return simple pair (no function)
    # $:y = :x # Simple pair
    # $:y = $:x # Simple pair
    # :y = f(:x) # re-write as simple call
    # :y = f($:x) # re-write as simple call, interpolation elsewhere
    # :y = :x + 1 # re-write as complicated call
    # :y = $:x + 1 # re-write as complicated call, interpolation elsewhere
    # $:y = f(:x) # re-write as simple call, unblock extract function
    # $:y = f($:x) # re-write as simple call, unblock, interpolation elsewhere
    # $y = :x + 1 # re-write as complicated col, unblock
    # $:y = $:x + 1 # re-write as complicated call, unblock, interpolation elsewhere
    # `@byrow` before any of the above
    ex, final_flags = extract_macro_flags(MacroTools.unblock(ex), deepcopy(outer_flags))
    check_macro_flags_consistency(final_flags)

    if gensym_names
        ex = Expr(:kw, gensym(), ex)
    end

    # :x
    # handled below via dispatch on ::QuoteNode

    # Fix any references to `cols` and replace them
    # with $
    if onearg(ex, :cols)
        ex = Expr(:$, ex.args[2])
    end

    # $:x
    if is_column_expr(ex)
        return ex.args[1]
    end

    if no_dest
        source, fun = get_source_fun(ex, exprflags = final_flags)
        return quote
            $source => $fun
        end
    end

    @assert ex.head == :kw || ex.head == :(=)
    lhs = ex.args[1]

    # fix cols
    if onearg(lhs, :cols)
        lhs = Expr(:$, lhs.args[2])
    end

    rhs = MacroTools.unblock(ex.args[2])

    if onearg(rhs, :cols)
        @warn "cols(x) is deprecated, use \$x instead"
        rhs = Expr(:$, rhs.args[2])
    end

    if is_macro_head(rhs, "@byrow")
        s = "In keyword argument inputs, `@byrow` must be on the left hand side. " *
        "Did you write `y = @byrow f(:x)` instead of `@byrow y = f(:x)`?"
        throw(ArgumentError(s))
    end

    # y = ...
    if lhs isa Symbol
        msg = "Using an un-quoted Symbol on the LHS is deprecated. " *
              "Write $(QuoteNode(lhs)) = ... instead."
        Base.depwarn(msg, "")
        lhs = QuoteNode(lhs)
    end

    # :y = :x
    if lhs isa QuoteNode && rhs isa QuoteNode
        source = rhs
        dest = lhs

        return quote
            $source => $dest
        end
    end

    # :y = $:x
    if lhs isa QuoteNode && is_column_expr(rhs)
        source = rhs.args[1]
        dest = lhs

        return quote
            $source => $dest
        end
    end

    # $:y = :x
    if is_column_expr(lhs) && rhs isa QuoteNode
        source = rhs
        dest = lhs.args[1]

        return quote
            $source => $dest
        end
    end

    # $:y = $:x
    if is_column_expr(lhs) && is_column_expr(rhs)
        source = rhs.args[1]
        dest = lhs.args[1]
        return quote
            $source => $dest
        end
    end

    # :y = f(:x)
    # :y = f($:x)
    # :y = :x + 1
    # :y = $:x + 1
    source, fun = get_source_fun(rhs; exprflags = final_flags)
    if lhs isa QuoteNode
        dest = lhs
        return quote
            $source => $fun => $dest
        end
    end

    # $:y = f(:x)
    if is_column_expr(lhs)
        dest = lhs.args[1]

        return quote
            $source => $fun => $dest
        end
    end

    throw(ArgumentError("This path should not be reached"))
end
fun_to_vec(ex::QuoteNode;
           no_dest::Bool=false,
           gensym_names::Bool=false,
           outer_flags::Union{NamedTuple, Nothing}=nothing) = ex

function make_source_concrete(x::AbstractVector)
    if isempty(x) || isconcretetype(eltype(x))
        return x
    elseif all(t -> t isa Union{AbstractString, Symbol}, x)
        return Symbol.(x)
    else
        throw(ArgumentError("Column references must be either all the same " *
                            "type or a a combination of `Symbol`s and strings"))
    end
end

protect_replace_syms!(e, membernames) = e
function protect_replace_syms!(e::Expr, membernames)
    if e.head == :quote
        e
    else
        replace_syms!(e, membernames)
    end
end

function replace_dotted!(e, membernames)
    x_new = replace_syms!(e.args[1], membernames)
    y_new = protect_replace_syms!(e.args[2], membernames)
    Expr(:., x_new, y_new)
end

function create_args_vector(args...; wrap_byrow::Bool=false)
    create_args_vector(Expr(:block, args...); wrap_byrow = wrap_byrow)
end

"""
   create_args_vector(arg) -> vec, outer_flags

Given an expression return a vector of operations
and a `NamedTuple` of the macro-flags that appear
in the expression.

If a `:block` expression, return the `args` of
the block as an array. If a simple expression,
wrap the expression in a one-element vector.
"""
function create_args_vector(arg; wrap_byrow::Bool=false)
    arg, outer_flags = extract_macro_flags(MacroTools.unblock(arg))

    if wrap_byrow
        if outer_flags[BYROW_SYM][]
            throw(ArgumentError("Redundant @byrow calls"))
        end

        outer_flags[BYROW_SYM][] = true
    end

    if arg isa Expr && arg.head == :block
        x = MacroTools.rmlines(arg).args
    else
        x = Any[arg]
    end

    return x, outer_flags
end


text:
#
# errors.jl -
#
# Management of errors for the Julia interface to the Spinnaker SDK.
#
#------------------------------------------------------------------------------

# Throws a `CallError` exception if `err` indicates an error in function `func`.
function _check(err::Cint, func::Symbol)
    if err != 0
        throw_call_error(err, func)
    end
    return nothing
end

const _error_symbols = Dict{Cint,Symbol}()

@noinline throw_call_error(err::Integer, func::Symbol) =
    throw(CallError(err, func))

show(io::IO, ::MIME"text/plain", err::CallError) =
    print(io, " Fail to execute", err.func)


text:
# This example still fails, not sure where it goes wrong

using Flux, Zygote
using Flux: @nograd, @adjoint, @epochs, onehotbatch, throttle, params
using Distributions
using Distributions: Gumbel
using Statistics: mean
using Plots

add_dim(x::Array) = reshape(x, (1,size(x)...))

function qammod(M)
    r = 1:sqrt(M)
    r = 2 * (r .- mean(r))
    r = [i for i in r, j in r]
    constellation = vcat(complex.(r, r')...)
    norm = sqrt(mean(abs2.(r)))
    constellation / norm
end

p_norm(p, x, fun) = sum(p .* fun(x))

function straight_through_estimator(x)
    M = size(x)[1]
    min_idx = argmin.(eachcol(x))
    Flux.onehotbatch(min_idx, 1:M)
end
# forward function, backward identity
@adjoint straight_through_estimator(x) = straight_through_estimator(x), x -> (x,)

TR = Float32
TC = ComplexF32

M = 64
constellation_dim = 2
N = 32 * M
temperature = 1.
SNR = 15
SNRlin = 10^(SNR/10) |> TR

nHidden = 128

constellation = qammod(M)
g_dist = Gumbel()

encoder = Chain(Dense(1, nHidden, Flux.relu),
                Dense(nHidden, M));
decoder = Chain(Dense(constellation_dim, nHidden, Flux.relu),
                Dense(nHidden, nHidden, Flux.relu),
                Dense(nHidden, M));

function gumbel_sample(M, N)
    rand(g_dist, (M, N))
end

function model(X)
    # sample from discrete distribution
    s_logits = encoder(X)
    g = gumbel_sample(M, N)
    s_bar = Flux.softmax((g .+ s_logits) / temperature)
    s = straight_through_estimator(s_bar)
    p_s = Flux.softmax(s_logits)

    # modulation
    norm_factor = sqrt(p_norm(p_s, constellation, x -> abs2.(x)))
    norm_constellation = constellation / complex.(norm_factor, 0)
    x = add_dim(norm_constellation) * complex.(s, 0)

    # Channel
    𝜎 = sqrt(1 / SNRlin) |> TR
    r = x + 𝜎 * randn(TC, 1, N)
    r = [real(r); imag(r)]

    # decoder
    Y = decoder(r)
    return p_s, s, Y, x, norm_constellation
end

stop_gradient(x) = x
@nograd stop_gradient

function loss(X)
    p_s, s, Y, x, norm_constellation = model(X)
    logit_loss = Flux.logitcrossentropy(Y, stop_gradient(s))
    entropy_x = -p_norm(p_s, p_s, x -> log2.(x))
    logit_loss - entropy_x
end

X = ones(1, 1)

@show loss(X)

ps = params(encoder, decoder);

gs = gradient(ps) do
    loss(X)
end

@show gs[encoder[1].W]

# opt = ADAM(0.001);
# data = [[X]]
# evalcb() = @show(loss(X));
# @epochs 500 Flux.train!(loss, ps, data, opt, cb = throttle(evalcb, 5));
#
# p_s, s, Y, x, norm_constellation = model(X)
#
# scatter(real(norm_constellation), imag(norm_constellation), aspect_ratio = :equal,
#                                                             markershape = :circle,
#                                                             markersize = 500*p_s,
#                                                             markerstrokealpha = 0)
# lim_ = 1.5
# ylims!((-lim_,lim_))
# xlims!((-lim_,lim_))
#
# @show p_s

max_stars_repo_name: Rassibassi/jaude

max_issues_repo_name: Rassibassi/jaude

max_forks_repo_name: Rassibassi/jaude


text:
"""
heqtest(n=100, c=.5, printh="yes"; jprecision=Float64, jmaxit=10)

Generates the examples for the SIREV paper.
"""
function heqtest(n = 100, c = 0.5, printh = "yes";
         jprecision = Float64, jmaxit = 10, hsham=1)
    FS = ones(n)
    x0 = ones(n)
    hdata = heqinit(x0, c);
    FPS = ones(jprecision, n, n)
    bargs=(atol = 1.e-20, rtol = 1.e-20, sham = hsham, resdec = .1, dx=1.e-6,
        pdata = hdata, maxit = jmaxit, stagnationok=true,
        armmax=0, printerr=false)
    nsoloutfd = nsol( heqf!, x0, FS, FPS; bargs...)
    nsoloutex = nsol( heqf!, x0, FS, FPS, heqJ!; bargs...)
    if printh == "yes"
        sol = nsoloutex.solution
        cpout = chandprint(sol, hdata)
        for i = 1:21
            @printf("%1.2e      %1.5e \n", cpout[i, 1], cpout[i, 2])
        end
        @printf("\n \n \n")
        println(nsoloutex.history)
    end
    return (exactout = nsoloutex, fdout = nsoloutfd)
end

"""
chandprint(x)

Print the table on page 125 (Dover edition) of Chandresekhar's book.
"""

function chandprint(x, pdata)
    c=pdata.cval[1]
    muc = collect(0:0.05:1)
    mu = pdata.mu
    n = length(mu)
    nx = length(x)
    LC = zeros(21, n)
    for j = 1:n
        for i = 1:21
            LC[i, j] = muc[i] / (muc[i] + mu[j])
        end
    end
    p = c * 0.5 / n
    LC = p * LC
    hout = LC * x
    onex = ones(size(muc))
    hout = onex ./ (onex - hout)
    return [muc hout]
end

function htime(n)
(x0, FS, FPS, FPS32, hdata) = hset(n);
println("analytic, Newton, Float64");
@time htest(x0, FS, FPS, hdata; analytic=true, hsham=1);
println("analytic, Newton, Float32");
@time htest(x0, FS, FPS32, hdata; analytic=true, hsham=1)
end

function htest(x0, FS, FPS, hdata; analytic=false, hsham=5)
    n=length(FS)
    #
    # I've preallocaed x0, FS, and FPS. But they may have been changed by previous runs.
    # The cost of resetting their entries to 1.0 is insignificant.
    #
    FS.=1.0
    FPS.=1.0
    bargs=(atol = 1.e-10, rtol = 1.e-10, sham = hsham, resdec = .1, pdata=hdata)
    if analytic
        nout=nsol( heqf!, x0, FS, FPS, heqJ!; bargs...)
    else
        nout=nsol( heqf!, x0, FS, FPS; bargs...)
    end
    return nout
end

function hset(n)
FS=ones(n,);
FPS=ones(n,n);
FPS32=ones(Float32,n,n);
x0=ones(n,);
c=.5;
hdata = heqinit(x0, c);
return(x0, FS, FPS, FPS32, hdata)
end

max_stars_repo_name: ctkelley/MPResults


text:
module TestSample
  using Test
  using EcologicalNetworks

  # Generate some data

  A = [true true false; true true true; false true false]
  N = UnipartiteNetwork(A)

  @test size(sample(N, 2)) == (2,2)
  @test size(sample(N, (2,2))) == (2,2)
  @test size(sample(N, (2,))) == (2,2)

  A = [true true false; true true true; false true false]
  N = BipartiteNetwork(A)

  @test size(sample(N, 2)) == (2,2)
  @test size(sample(N, (2,))) == (2,2)
  @test size(sample(N, (2,1))) == (2,1)

end


text:
using MarketTechnicals, MarketData, FactCheck

facts("Momentum") do

    context("rsi") do
        @fact rsi(cl, 14).values[end]              --> roughly(55.84921708898833) # TTR value is 55.84922
        @fact rsi(cl, 14, wilder=true).values[end] --> roughly(55.95931797240996) # TTR value is 55.95932
        @fact rsi(cl, 14).timestamp[end]           --> Date(2001,12,31)
    end

    context("macd") do
        @fact macd(cl).values[end, 1] --> roughly(0.42117515229271874) # TTR value with percent=FALSE is 0.421175152
        @fact macd(cl).values[end, 2] --> roughly(0.44142752613313185) # TTR value with percent=FALSE is 0.4414275
        @fact macd(cl).timestamp[end] --> Date(2001,12,31)
    end

#     context("cci") do
#         @fact cci(ohlc).values[1]      --> roughly(-38.931614)      # TTR::CCI value is -38.931614
#         @fact cci(ohlc).values[end]    --> roughly(46.3511339)      # TTR::CCI value is 46.3511339
#         @fact cci(ohlc).timestamp[end] --> Date(2001,12,31)
#     end
end


text:
using StochasticDiffEq, DiffEqNoiseProcess, Test, Random, LinearAlgebra

function true_general_function(dt,dW,C,m)
  # Choose p
  sum_dW² = dot(dW,dW)
  M = div(m*(m-1),2)
  p = Int(floor(sqrt(M/(12*dt*C))*sqrt(m + 4*sum_dW²/dt)/π + 1))
  # Alternative choice of p which on average results in larger values:
  #p = Int(floor(sqrt(5*m*M/(12*dt*C))/π+1))
  #@show p
  # a_p below (22) for tail approx.
  ap = pi^2/6
  for r=1:p
    ap -= 1/r^2
  end

  # define permutation operator
  PermOp = zeros(Int,m^2,m^2)
  for i=1:m^2
    j = 1 + m*((i-1)%m) + div(i-1,m)
    PermOp[i,j] = 1
  end

  # define K operator
  KOp = zeros(Int,M,m^2)
  row = 1
  for i=1:m-1
    col = (i-1)*m + i + 1
    s = m-i
    KOp[row:(row+s)-1,col:(col+s)-1] = Matrix{eltype(dW)}(I, s, s)
    row += s
  end

  # define identiy matrices
  Idm = Matrix{eltype(dW)}(I, m, m)
  IdM = Matrix{eltype(dW)}(I, M, M)
  Idm2 = Matrix{eltype(dW)}(I, m^2, m^2)

  # compute Σinf
  op1 = 2/dt * KOp*(Idm2-PermOp)
  op2 = kron(Idm, dW .* dW') #kron(I+zeros(m,m),W.dW .* W.dW')
  op3 = (Idm2-PermOp)*KOp'
  Σinf = 2*IdM + op1*op2*op3

  # compute sqrt( Σinf ) from Cholesky decomposition
  α = sqrt(1 + sum_dW²/dt)
  SqΣinf = (Σinf + 2*α*IdM)/(sqrt(2)*(1+α))

  # initialze Gp from distribution N(0, I_M) to compute tail approx
  Gp₁ = randn(M)
  prefac = (dt/(2*π))*sqrt(ap)
  Atail = prefac*SqΣinf*Gp₁ # check influence of op3 here!!

  #@show Atail

  # compute Atilde
  Ap = zeros(M)

  # for r=1:p
  #   sd = sqrt(dt/(2*pi*r))
  #   ζ = randn(m)*sd
  #   η = randn(m)*sd
  #   aᵢ₀ -= (2/sqrt(π*r))*ζ
  #   for i=1:m
  #     for j=1:m
  #       Ap[i,j] += ζ[i]*η[j]-η[i]*ζ[j]
  #     end
  #   end
  # end
  # for i=1:m
  #   for j=1:m
  #     Jp[i,j] += Ap[i,j] - 1//2*(aᵢ₀[j]*dW[i]-dW[j]*aᵢ₀[i])
  #   end
  # end

  op4 = KOp*(PermOp-Idm2)

  for k in 1:p
    ζ = randn(eltype(dW),m)
    η = randn(eltype(dW),m)
    #@show ζ , η
    #@show ζ * sqrt(dt/(2*π*k)), η * sqrt(dt/(2*π*k))
    #kronprod1 = vec(ζ) .* vec(η+convert(eltype(dW),sqrt(2/dt))*dW)'
    kronprod2 = vec(η+convert(eltype(dW),sqrt(2/dt))*dW) .* vec(ζ)'

    #@show size(vec(kronprod2)), size(op4), size(Ap)
    Ap += op4*vec(kronprod2)/k #(kronprod1 - kronprod2)/k
  end

  Ap *= dt/(2*π)

  # compute J_{ij}^p
  Jp_com = vec(1//2 .* vec(dW) .* vec(dW)')

  Jp =  Jp_com + op3*(Ap + Atail)

  #@show op3*Atail

  return reshape(Jp, (m, m)), reshape((op3*Ap), (m, m)), reshape(op3*Atail, (m, m))
end

seed = 10
Random.seed!(seed)

m = 10
W = WienerProcess(0.0,zeros(m),nothing)

dt = 0.1
calculate_step!(W,dt,nothing,nothing)

for i in 1:10
  accept_step!(W,dt,nothing,nothing)
end

# Diagonal noise tests

true_diag = 1//2 .* W.dW .* W.dW

Wikdiag = StochasticDiffEq.WikJDiagonal_iip(W.dW)
StochasticDiffEq.get_iterated_I!(dt, W.dW, W.dZ, Wikdiag, nothing, 1)
Wikdiagoop = StochasticDiffEq.WikJDiagonal_oop()

@test StochasticDiffEq.get_iterated_I(dt, W.dW, W.dZ, Wikdiagoop, nothing, 1) == Wikdiag.WikJ
@test Wikdiag.WikJ == true_diag

KPWdiagiip = StochasticDiffEq.KPWJ_iip(W.dW)
Random.seed!(seed)
StochasticDiffEq.get_iterated_I!(dt, W.dW, W.dZ, KPWdiagiip, Int(1e3), 1, 1//1)
KPWdiagoop = StochasticDiffEq.KPWJ_oop()
Random.seed!(seed)
@test isapprox(StochasticDiffEq.get_iterated_I(dt, W.dW, W.dZ, KPWdiagoop, Int(1e3), 1, 1//1), KPWdiagiip.WikJ, atol=1e-14)

StochasticDiffEq.get_iterated_I(dt, W.dW, W.dZ, KPWdiagoop, Int(1e3), 1, 1//1)

KPWdiagonly = [KPWdiagiip.WikJ[i, i] for i in 1:m]

@test KPWdiagonly == true_diag
@test KPWdiagonly == Wikdiag.WikJ

# Commutative noise tests

true_commute = 1//2 .* W.dW .* W.dW'

Wikcommute = StochasticDiffEq.WikJCommute_iip(W.dW)
StochasticDiffEq.get_iterated_I!(dt, W.dW, W.dZ, Wikcommute, nothing, 1)
Wikcommuteoop = StochasticDiffEq.WikJCommute_oop()

@test StochasticDiffEq.get_iterated_I(dt, W.dW, W.dZ, Wikcommuteoop, nothing, 1) == Wikcommute.WikJ
@test Wikcommute.WikJ == true_commute

# General noise tests

Random.seed!(seed)
Wikgeneral = StochasticDiffEq.WikJGeneral_iip(W.dW)
StochasticDiffEq.get_iterated_I!(dt, W.dW, W.dZ, Wikgeneral, nothing, 1)

Random.seed!(seed)
Wikgeneraloop  = StochasticDiffEq.WikJGeneral_oop(W.dW)
@test StochasticDiffEq.get_iterated_I(dt, W.dW, W.dZ, Wikgeneraloop, nothing, 1) == Wikgeneral.WikJ

Random.seed!(seed)
true_noncom = true_general_function(dt, W.dW, 1.0, m)

@test isapprox(Wikgeneral.WikJ, true_noncom[1], atol=1e-15)
@test isapprox(Wikgeneral.WikJ2, true_noncom[3], atol=1e-15)

# Test the relations given in Wiktorsson Eq.(2.1)
Aij = true_noncom[2]+true_noncom[3]
Aii = [Aij[i, i] for i in 1:m]
Jii = [true_noncom[1][i, i] for i in 1:m]

@test isapprox(true_noncom[1] + true_noncom[1]', 2*true_commute)
@test isapprox(true_noncom[1], true_commute + Aij)
@test isapprox(Aij, -Aij', atol=1e-15)
@test isapprox(Jii, 1//2 .* W.dW .* W.dW, atol=1e-15)
@test isapprox(Aii, zeros(m), atol=1e-15)

# .. for KPW scheme
Aii = [KPWdiagiip.WikA[i, i] for i in 1:m]
@test isapprox(KPWdiagiip.WikJ + KPWdiagiip.WikJ', 2*true_commute)
@test !isapprox(KPWdiagiip.WikJ, true_commute + KPWdiagiip.WikA) # ! because of correction with \mu term
@test isapprox(KPWdiagiip.WikA, -KPWdiagiip.WikA', atol=1e-15)
@test isapprox(KPWdiagonly, 1//2 .* W.dW .* W.dW, atol=1e-15)
@test isapprox(Aii, zeros(m), atol=1e-15)

# moment conditions
# E1(I_{j1} I_{j2}) = Δ δ_{j1,j2}
# E2(I_{j1, j2} I_{j1, j2}) = 1/2 Δ^2
# E3(I_{j1} I_{j2} I_{j3, j4}) = {Δ^2 if j1=j2=j3=j4, 1/2 Δ^2 if j3!=j4 and j1=j3, j2=j4 or j1=j4, j2=j3, 0 otherwise}

function moments!(tmp, E1, E2, E3, W, Wik, Δ, samples, p=nothing)
  for _ in 1:samples
    calculate_step!(W,Δ,nothing,nothing)
    accept_step!(W,Δ,nothing,nothing)
    mul!(tmp,vec(W.dW),vec(W.dW)')

    StochasticDiffEq.get_iterated_I!(Δ, W.dW, W.dZ, Wik, p, 1)
    @. E1 = E1 + tmp
    @. E2 = E2 + Wik.WikJ * Wik.WikJ
    @. E3 = E3 + tmp * Wik.WikJ
  end
  @. E1 = E1/samples
  @. E2 = E2/samples
  @. E3 = E3/samples

  return nothing
end

samples = Int(5e5)
E1 = false .* vec(W.dW) .* vec(W.dW)'
E2 = zero(E1)
E3 = zero(E1)
tmp = zero(E1)
# generate new random dW
W = WienerProcess!(0.0,zeros(m),nothing)
Wikmom = StochasticDiffEq.WikJGeneral_iip(W.dW)
@time moments!(tmp, E1, E2, E3, W, Wikmom, dt, samples)

@test maximum(abs.(E1 - dt*I)) < 1e-3
@test maximum(abs.(E2- dt^2*(zero(E2).+1)/2)) < 4e-3
@test maximum(abs.(E3 - dt^2*(one(E2).+1)/2)) < 6e-3

"""
Problem 2.3.3 from
Kloeden, P. E., Platen, E., & Schurz, H. Numerical solution of SDE through computer
experiments. Springer Science & Business Media. (2012)
"""

Random.seed!(seed)
dt = 1.0
ps = [1,Int(1e2),Int(5e3),Int(1e4)]
m = 2
W2 = WienerProcess(0.0,zeros(m),nothing)
calculate_step!(W2,dt,nothing,nothing)
for i in 1:10
  accept_step!(W2,dt,nothing,nothing)
end

function path_convergence(dt, ps, W, Wik)
    sample_path = []
    for (i, p) in enumerate(ps)
      Random.seed!(seed)
      @show p
      StochasticDiffEq.get_iterated_I!(dt, W.dW, W.dZ, Wik, p, 1)
      #@show Wik.WikJ
      push!(sample_path,Wik.WikJ[1,2])
    end
    return sample_path
end

Wikpath = StochasticDiffEq.KPWJ_iip(W2.dW)
path = path_convergence(dt, ps, W2, Wikpath)

v = abs.(path .- path[end])
@test sort(v, rev = true) == v

"""
Exercise 2.3.1 from
Kloeden, P. E., Platen, E., & Schurz, H. Numerical solution of SDE through computer
experiments. Springer Science & Business Media. (2012)
"""

mutable struct StatsJ{AType}
  mean::Vector{AType}
  var::Vector{AType}
  meanold::Vector{AType}
  varold::Vector{AType}
  tmp::AType
  x::AType
end

function StatsJ(u)
  mean = Vector{typeof(u)}()
  var = Vector{typeof(u)}()
  meanold = Vector{typeof(u)}()
  varold = Vector{typeof(u)}()

  for k=1:2
    push!(mean,zero(u))
    push!(var,zero(u))
    push!(meanold,zero(u))
    push!(varold,zero(u))
  end

  tmp = zero(u)
  x = zero(u)
  StatsJ(mean,var,meanold,varold,tmp,x)
end

function Welford!(cache::StatsJ, W, Wik, Δ, samples, p=nothing)
  for k in 1:2
    fill!(cache.var[k], zero(eltype(cache.tmp)))
  end
  for i in 1:samples
    calculate_step!(W,Δ,nothing,nothing)
    accept_step!(W,Δ,nothing,nothing)
    mul!(cache.tmp,1//2*vec(W.dW),vec(W.dW)')

    StochasticDiffEq.get_iterated_I!(Δ, W.dW, W.dZ, Wik, p, 1)
    for k in 1:2
      if k==1
        copyto!(cache.x,cache.tmp)
      else
        copyto!(cache.x,Wik.WikJ)
      end
      if i == 1
        copyto!(cache.meanold[k], cache.x)
      else
        @. cache.mean[k] = cache.meanold[k] + (cache.x - cache.meanold[k]) / i
        @. cache.var[k] = cache.varold[k] + (cache.x - cache.meanold[k]) * (cache.x - cache.mean[k])

        copyto!(cache.meanold[k],cache.mean[k])
        copyto!(cache.varold[k],cache.var[k])
      end
    end
  end

  for k in 1:2
    @. cache.mean[k] = cache.mean[k]/samples
    @. cache.var[k] = cache.var[k]/(samples-1)
  end

  return nothing
end

cache = StatsJ(false .* vec(W2.dW) .* vec(W2.dW)')
@time Welford!(cache, W2, Wikpath, dt, Int(1e5), ps[2])

@test maximum(abs.(cache.mean[1])) < 1e-5
@test maximum(abs.(cache.mean[2])) < 1e-5
@test maximum(abs.(cache.var[1]-cache.var[2])) > 1e-1


text:
printstyled(stdout,"  show\n", color=:light_green)

# SeisChannel show
S = SeisChannel()
nx = (1, 2, 3, 4, 5, 10, 100, 10000)
redirect_stdout(out) do
  for i in nx
    S.t = [1 0; i 0]
    S.x = randn(i)
    show(S)
  end
end

redirect_stdout(out) do
  # show
  show(breaking_seis())
  show(randSeisData(1))
  show(SeisChannel())
  show(SeisData())
  show(randSeisChannel())

  # summary
  summary(randSeisChannel())
  summary(randSeisData())

  # invoke help-only functions
  @test seed_support() == nothing
  @test chanspec() == nothing
  @test mseed_support() == nothing
  @test timespec() == nothing
  @test RESP_wont_read() == nothing
end


text:
abstract type AbstractIntegrationMethod end
struct AutoQuadrature end
struct Quadrature end
struct TrapezoidRule end

"""
    SFMDistribution(pdf::Function, dq::DiscretisedFluidQueue{<:Mesh}, fun_evals::Int = 6)

Return a discretised version of the `pdf`. The method of discretisation depends on the discretisation method
of `dq`.

# Arguments:
- `pdf::Function`: a function `f(x::Float64,i::Int)` where `f(x,i)dx=P(X(0)∈dx,φ(0)=i)` is the initial
    distribution of a fluid queue.
- `dq::DiscretisedFluidQueue{<:Mesh}`:
- `fun_evals`: the number of function evaluations used to approximate `f(x,i)`
"""
function SFMDistribution(pdf::Function,dq::DiscretisedFluidQueue{<:Mesh},method::Type{AbstractIntegrationMethod})
    throw(DomainError("No such integration method"))
end

"""
    SFMDistribution(pdf::Function, dq::DiscretisedFluidQueue{DGMesh{T}})

Approximates `pdf` via polynomials using interpolation.
"""
function SFMDistribution(pdf::Function,dq::DiscretisedFluidQueue{DGMesh{T}},method::Type{AutoQuadrature}=AutoQuadrature) where T
    cellnodes = cell_nodes(dq)
    n₋ = N₋(dq)
    coeffs = zeros(n_bases_per_cell(dq),n_phases(dq),n_intervals(dq))
    for i in phases(dq)
        for cell in 1:n_intervals(dq)
            nodes = cellnodes[:,cell]
            if length(nodes)==1
                weights = Δ(dq,cell)
            else
                weights = gauss_lobatto_weights(nodes[1],nodes[end],length(nodes))
            end
            coeffs[:,i,cell] = pdf.(nodes,i).*weights
        end
    end
    coeffs = [zeros(n₋); coeffs[:]; zeros(N₊(dq))]
    return SFMDistribution(coeffs,dq)
end

function SFMDistribution(pdf::Function,dq::DiscretisedFluidQueue{DGMesh{T}},method::Type{TrapezoidRule};fun_evals::Int=2001) where T
    n₋ = N₋(dq)
    coeffs = zeros(n_bases_per_cell(dq),n_phases(dq),n_intervals(dq))
    V = vandermonde(n_bases_per_cell(dq))
    for i in phases(dq)
        for cell in 1:n_intervals(dq)
            nodes = dq.mesh.nodes[cell:cell+1]
            if n_bases_per_cell(dq)==1
                weights = fill(Δ(dq,cell),1,1)
            else
                weights = gauss_lobatto_weights(nodes[1],nodes[end],n_bases_per_cell(dq))
            end
            lagrange_polynomials(x) = DiscretisedFluidQueues.lagrange_polynomials(
                DiscretisedFluidQueues.gauss_lobatto_points(nodes[1],nodes[end],n_bases_per_cell(dq)),x)./weights

            eval_pts = range(nodes[1],nodes[end],length=fun_evals)
            h = Δ(dq,cell)/(fun_evals-1)
            vals = zeros(n_bases_per_cell(dq),fun_evals)
            for (m,x) in enumerate(eval_pts)
                vals[:,m] = lagrange_polynomials(x).*pdf(x,i)*h
            end
            lagrange_coeffs = (sum(vals[:,2:end-1],dims=2)+vals[:,1]/2+vals[:,2]/2)'*(weights[:].*V.V*V.V'.*weights[:]')*2.0/Δ(dq,cell)
            coeffs[:,i,cell] = lagrange_coeffs
        end
    end
    coeffs = [zeros(n₋); coeffs[:]; zeros(N₊(dq))]
    return SFMDistribution(coeffs,dq)
end

"""
    interior_point_mass(x::Float64, i::Int, dq::DiscretisedFluidQueue{<:Mesh})

Construct an approximation to a distribution of fluid queue with a point mass at (x,i).
Returns a SFMDistribution. The method of approximation is defined by the type of <:Mesh.

# Arguments:
- `x::Float64`: The position of the point mass
- `i::Int`: The phase of the point mass
- `dq::DiscretisedFluidQueue`:
"""
function interior_point_mass(x::Float64,i::Int,dq::DiscretisedFluidQueue{<:Mesh}) end
"""
    interior_point_mass(x::Float64, i::Int, dq::DiscretisedFluidQueue{DGMesh{T}})

Constructs a polynomial approximation to the point mass at (x,i)
"""
function interior_point_mass(x::Float64,i::Int,dq::DiscretisedFluidQueue{DGMesh{T}}) where T
    (x<dq.mesh.nodes[1])&&throw(DomainError("x is not in interior"))
    (x>dq.mesh.nodes[end])&&throw(DomainError("x is not in interior"))
    if _has_right_boundary(dq.model.S,i)
        cell_idx, cellnodes, coeff_idx = _get_coeff_index_pos(x,i,dq)
    elseif _has_left_boundary(dq.model.S,i)
        cell_idx, cellnodes, coeff_idx = _get_coeff_index_neg(x,i,dq)
    end
    n₊ = N₊(dq)
    n₋ = N₋(dq)
    coeffs = zeros(n₊+n₋+total_n_bases(dq))
    nodes = cell_nodes(dq.mesh)[:,cell_idx]
    V = vandermonde(n_bases_per_cell(dq))
    # we solve a∫ψ(x)ψ(x)'dx = ∫δ(x-x₀)ψ'(x)dx where ψ(x) is a column vector of
    # lagrange polynomials on the current cell where each basis function integrates to 1.
    # The right-hand integral is ψ'(x₀). Let W = diag(V.w) and note W=W' and W^-1=(W^-1)'
    # Transform to a interpolating basis via ϕ(x) = (Δ/2×W)*ψ(x) where Δ/2V.w are the
    # integrals of the interpolating basis functions ϕ(x).
    # The inner product on the left-hand side
    # becomes (2/Δ)W^-1 * ∫ϕ(x)ϕ(x)'dx * (2/Δ)W^-1.
    # Maths tells us that ∫ϕ(x)ϕ(x)'dx = (V*V')^-1 Δ/2.
    # Thus a∫ψ(x)ψ(x)'dx = ∫δ(x-x₀)ψ'(x)dx
    # becomes a(2/Δ)W^-1 * ∫ϕ(x)ϕ(x)'dx * (2/Δ)W^-1 = ϕ'(x₀)(2/Δ)W^-1
    # solving for a
    # a = ϕ'(x₀) (2/Δ)W^-1 * ((2/Δ)W^-1*∫ϕ(x)ϕ(x)'dx*(2/Δ)W^-1)^-1
    #   = ϕ'(x₀) (∫ϕ(x)ϕ(x)'dx)^-1 * (Δ/2)W
    #   = ϕ'(x₀) (V.V*V.V')2/Δ * (Δ/2)W
    #   = ϕ'(x₀) (V.V*V.V') W
    coeffs[coeff_idx] = lagrange_polynomials(nodes,x)'*(V.V*V.V').*V.w'
    return SFMDistribution(coeffs,dq)
end

"""
    left_point_mass(i::Int, dq::DiscretisedFluidQueue)

Construct a SFMDistribution with a point mass at the left boundary of `dq` in phase `i` and 0 elsewhere.
"""
function left_point_mass(i::Int,dq::DiscretisedFluidQueue)
    !_has_left_boundary(dq.model.S,i)&&throw(DomainError("only phases with lpm=true have left point masses"))
    n₊ = N₊(dq)
    n₋ = N₋(dq)
    coeffs = zeros(n₊+n₋+total_n_bases(dq))
    nᵢ = N₋(dq.model.S[1:i])
    coeffs[nᵢ] = 1.0
    return SFMDistribution(coeffs,dq)
end

"""
    right_point_mass(i::Int, dq::DiscretisedFluidQueue)

Construct a SFMDistribution with a point mass at the left boundary of `dq` in phase `i` and 0 elsewhere.
"""
function right_point_mass(i::Int,dq::DiscretisedFluidQueue)
    !_has_right_boundary(dq.model.S,i)&&throw(DomainError("only phases with rpm=false have left point masses"))
    n₊ = N₊(dq)
    n₋ = N₋(dq)
    coeffs = zeros(n₊+n₋+total_n_bases(dq))
    nᵢ = N₊(dq.model.S[1:i])
    coeffs[end-n₊+nᵢ] = 1.0
    return SFMDistribution(coeffs,dq)
end

"""
    SFMDistribution(pdf::Function, dq::DiscretisedFluidQueue{FVMesh}, fun_evals::Int = 6)

Construct an approximation to `pdf` as the average of `pdf` of each cell.

Uses quadrature with `fun_evals` function evaluations to approximate cell averages.
"""
function SFMDistribution(pdf::Function,dq::DiscretisedFluidQueue{FVMesh{T}},method::Type{Quadrature}=Quadrature;fun_evals::Int=6) where T
    n₋ = N₋(dq)
    coeffs = zeros(n_bases_per_cell(dq),n_phases(dq),n_intervals(dq))
    for i in phases(dq)
        for cell in 1:n_intervals(dq)
            a,b = dq.mesh.nodes[cell:cell+1]
            quad = gauss_lobatto_quadrature(x->pdf(x,i),a,b,fun_evals)
            coeffs[1,i,cell] = quad./Δ(dq,cell)
        end
    end
    coeffs = [zeros(n₋); coeffs[:]; zeros(N₊(dq))]
    return SFMDistribution(coeffs,dq)
end

"""
    interior_point_mass(x::Float64, i::Int, dq::DiscretisedFluidQueue{FVMesh})

Construct an approximation to a point mass at (x,i). Basically, just puts mass 1 in the cell containing (x,i).
"""
function interior_point_mass(x::Float64,i::Int,dq::DiscretisedFluidQueue{FVMesh{T}}) where T
    if _has_right_boundary(dq.model.S,i)
        cell_idx, ~, ~ = _get_coeff_index_pos(x,i,dq)
    elseif _has_left_boundary(dq.model.S,i)
        cell_idx, ~, ~ = _get_coeff_index_neg(x,i,dq)
    end
    n₊ = N₊(dq)
    n₋ = N₋(dq)
    coeffs = zeros(n₊+n₋+total_n_bases(dq))
    coeffs[cell_idx] = 1.0./Δ(dq,cell_idx)
    return SFMDistribution(coeffs,dq)
end

## constructors for FRAPMesh
function SFMDistribution_from_cdf(cdf::Function,dq::DiscretisedFluidQueue{FRAPMesh{T}},method::Type{TrapezoidRule}=TrapezoidRule;fun_evals::Int=2001) where T
    n₋ = N₋(dq)
    coeffs = zeros(n_bases_per_cell(dq),n_phases(dq),n_intervals(dq))
    for i in phases(dq)
        for cell in 1:n_intervals(dq)
            a,b = dq.mesh.nodes[cell:cell+1]
            if _has_right_boundary(dq.model.S,i)
                o = expected_orbit_from_cdf(x->cdf(Δ(dq,cell)-x,i),dq.mesh.me,a,b,fun_evals)
            elseif _has_left_boundary(dq.model.S,i)
                o = expected_orbit_from_cdf(x->cdf(x,i),dq.mesh.me,a,b,fun_evals)
            end
            coeffs[:,i,cell] = o
        end
    end
    coeffs = [zeros(n₋); coeffs[:]; zeros(N₊(dq.model.S))]
    return SFMDistribution(coeffs,dq)
end

"""
    SFMDistribution(pdf::Function,dq::DiscretisedFluidQueue{FRAPMesh{T}},method::Type{TrapezoidRule}=TrapezoidRule;fun_evals=2001) where T

Return the appropriate initial condition to approximate the initial distribution `pdf`
for the numerical discretisation scheme defined by the `dq` DiscretisedFluidQueue. `method`
specifies the integration rule to approximate the integral of the pdf and the orbit function.

i.e. for each cell compute the expected orbit position
``∫pdf(x,i) ⋅ a(x) dx ``
where the integral is over each cell and
``a(x) = a exp(S(x-yₖ))``
if the membership of `i` is `-1` and
``a(x) = a exp(S(yₖ₊₁-x))``
if the membership of `i` is `1` where
`yₖ` and `yₖ₊₁` are the left and right cell edges, respectively, and
`a` and `S` are defined by `dq.mesh.me` are parameters of a MatrixExponential.
"""
function SFMDistribution(pdf::Function,dq::DiscretisedFluidQueue{FRAPMesh{T}},method::Type{TrapezoidRule}=TrapezoidRule;fun_evals=2001) where T
    n₋ = N₋(dq)
    coeffs = zeros(n_bases_per_cell(dq),n_phases(dq),n_intervals(dq))
    for i in phases(dq)
        for cell in 1:n_intervals(dq)
            a,b = dq.mesh.nodes[cell:cell+1]
            if _has_right_boundary(dq.model.S,i)
                o = expected_orbit_from_pdf(x->(b-a)*pdf(a+x*(b-a),i),dq.mesh.me,0.0,1.0,fun_evals)
            else
                o = expected_orbit_from_pdf(x->(b-a)*pdf(b-x*(b-a),i),dq.mesh.me,0.0,1.0,fun_evals)
            end
            coeffs[:,i,cell] = o
        end
    end
    coeffs = [zeros(n₋); coeffs[:]; zeros(N₊(dq))]
    return SFMDistribution(coeffs,dq)
end

"""
    interior_point_mass(x::Float64, i::Int, dq::DiscretisedFluidQueue{FRAPMesh})

Return the appropriate initial condition to approximate the initial distribution which
is a point mass at (x,i) for the numerical discretisation scheme defined by the
`dq` DiscretisedFluidQueue.

i.e. Compute the orbit position
``a(x) = a exp(S(x-yₖ))``
if the membership of `i` is `-1` and
``a(x) = a exp(S(yₖ₊₁-x))``
if the membership of `i` is `1` where
`yₖ` and `yₖ₊₁` are the left and right cell edges and `x∈[yₖ,yₖ₊₁]`, respectively, and
`a` and `S` are defined by `dq.mesh.me` are parameters of a MatrixExponential.
"""
function interior_point_mass(x::Float64,i::Int,dq::DiscretisedFluidQueue{FRAPMesh{T}}) where T
    if _has_right_boundary(dq.model.S,i)
        cell_idx, ~, coeff_idx = _get_coeff_index_pos(x,i,dq)
        yₖ = dq.mesh.nodes[cell_idx]
        d = x-yₖ
    elseif _has_left_boundary(dq.model.S,i)
        cell_idx, ~, coeff_idx = _get_coeff_index_neg(x,i,dq)
        yₖ₊₁ = dq.mesh.nodes[cell_idx+1]
        d = yₖ₊₁-x
    end
    o = orbit(dq.mesh.me,d/Δ(dq.mesh,cell_idx))
    coeffs = zeros(N₊(dq)+N₋(dq)+total_n_bases(dq))
    coeffs[coeff_idx] = o
    return SFMDistribution(coeffs,dq)
end

max_stars_repo_name: angus-lewis/StochasticFluidQueues

max_issues_repo_name: angus-lewis/StochasticFluidQueues

max_forks_repo_name: angus-lewis/StochasticFluidQueues


text:
import .MonteCarloMeasurements: Particles
# From simulation result to particles

vv2m(x) = copy(reduce(hcat, x)')
function simulate(pf, n, du, Npart)
    sims = map(_->vv2m.(simulate(pf, n, du)), 1:Npart)
    ns = length(sims[1]) # number of time series, usually x,u,y
    ntuple(ns) do s
        xes = getindex.(sims,s)
        M = reduce((x,y)->cat(x,y,dims=3), xes)
        M = permutedims(M, (3,2,1))
        copy(dropdims(mapslices(Particles, M, dims=(1,2)), dims=2)')
    end
end

# From filter result to particles
function MonteCarloMeasurements.Particles(x::AbstractMatrix{<:AbstractVector},we=nothing) # Helper function
    xp = copy(vv2m([Particles(vv2m(c)) for c in eachcol(x)]))
    we === nothing && return xp
    ## Perform resampling so that all particles have the same weight
    choices = map(LowLevelParticleFilters.resample, eachcol(we))
    for j in 1:size(xp,2), i in 1:size(xp,1)
        xp[i,j].particles .= xp[i,j].particles[choices[j]]
    end
    xp
end


text:
module WindingNumbers

include("Points.jl")
include("Vertices.jl")
include("Implementation.jl")

export windingnumber, windingangle

end

max_stars_repo_name: jwscook/WindingNumbers

max_issues_repo_name: jwscook/WindingNumbers

max_forks_repo_name: jwscook/WindingNumbers


text:
todisplay = VERSION <= VersionNumber("1.2") ? true : false

using LinearAlgebra
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using CAGD
if todisplay
    using ViewerGL
    GL = ViewerGL
	include("views.jl")
end

# Rod Generation

# npts = 4
npts = 16
V,_ = Lar.rod()([npts, 1])
EV = [
    [[2*i, (2*i+1)%2npts+1] for i = 1 : npts];             # horizontal upper edges
    [[2*i-1, (2*i)%2npts+1] for i = 1 : npts];           # horizontal lower edges
    [[2*i+1, 2*i+2] for i = 0 : npts-1];                # vertical edges
]
FE = [
    [i for i = 1 : npts],                                       # upper face
    [npts + i for i = 1 : npts],                                # lower face
    [[i, npts+i, 2*npts+i, 2*npts+i+1] for i = 1:npts-1]...,    # vertical faces
    [npts, 2*npts, 3*npts, 2*npts+1]
]
FV = Lar.cop2lar(Lar.lar2cop(FE) * Lar.lar2cop(EV))
CV = [collect(1 : size(V,2))]

if todisplay
    GL.VIEW([
        GL.GLGrid(V,EV,GL.COLORS[1],1)
        GL.GLFrame
    ]);
end

# Three Rods Generation

cyl = Lar.Struct([
    Lar.r(0,0,0),
    #Lar.s(0.6, 0.6, 2.0),
    Lar.s(1.0, 1.0, 1.5),
    Lar.t(0,0,-1.5),
    (V,CV,FV,EV)
])
cyl = Lar.struct2lar(cyl)
tris = Lar.Struct([
    cyl,
    Lar.Struct([Lar.r(pi/2,0,0), cyl ]) ,
    Lar.Struct([Lar.r(0,pi/2,0), cyl ])
])
V, CV, FV, EV = Lar.struct2lar(tris)
cyls = (V, CV, FV, EV)

# Cube Generation

V, (_, EV, FV, CV) = Lar.cuboid([1,1,1],true,[-1,-1,-1])
cube = (V, CV, FV, EV)

# Sphere Generation

# V, EV, FV = catmullclark(cube[1], cube[3], cube[2], 4)
# sphere = (V, FV, EV)

V, FV = Lar.sphere(1.0)([15,25])
# V, FV = Lar.sphere(1.0)([5,5])
FV = sort(sort.(FV))
function getFaceEdges(fV::Array{Int,1})
    return [fV[1], fV[2]], [fV[1], fV[3]], [fV[2], fV[3]]
end
EVs = unique([(map(f -> getFaceEdges(f), FV)...)...])
CV = [collect(1:size(V, 2))]
sphere = (V, CV, FV, EVs)

# Object Creation

carry = Lar.Struct([
    cyls,
    Lar.s(1.3,1.3,1.3),
    cube,
    Lar.s(1.4,1.4,1.4),
    #Lar.s(1.5,1.5,1.5)
    sphere
])
V,CV,FV,EV = Lar.struct2lar(carry)

if todisplay
    GL.VIEW([
        GL.GLGrid(V,EV,GL.COLORS[1],1)
        GL.GLFrame
    ]);
end

# Model Generation

model = CAGD.Model(V)
cFE = convert(Lar.ChainOp, Lar.coboundary_1(model.G, FV, EV))
CAGD.addModelCells!(model, 1, EV, signed = true)
CAGD.addModelCells!(model, 2, cFE)
# Since CF is only used in boolean evaluation, the sign is not needed:
cFV = Lar.lar2cop(FV)
cCV = Lar.lar2cop(CV)
cCF = convert(Lar.ChainOp, ((cCV*cFV').>0))
CAGD.addModelCells!(model, 3, cCF)

atol = 1e-9;
if todisplay  displayModel(model)  end

split_model = CAGD.facesplitting(model, atol = atol)

if todisplay  displayModel(split_model)  end

congr_model = CAGD.mergeModelVertices(split_model, err=atol, signed_merge=true)
#split_model = CAGD.facesplitting(congr_model, atol = atol)
#congr_model = CAGD.mergeModelVertices(split_model, err=atol, signed_merge=true)

if todisplay  displayModel(congr_model)  end

gift_model = deepcopy(congr_model);
FC, bicon_comps = CAGD.tgw(congr_model, 3)
CAGD.addModelCells!(gift_model, 3, convert(Lar.ChainOp, FC'))

if todisplay  viewExplode(gift_model)  end

##==============================================================================
##  Boolean Decomposition
##==============================================================================

arranged_model, boolean_matrix = CAGD.bool3(model)

bXRod = boolean_matrix[:, 2]
bYRod = boolean_matrix[:, 3]
bZRod = boolean_matrix[:, 4]
bCube = boolean_matrix[:, 5]
bSphe = boolean_matrix[:, 6]


text:
module C

value = 3

end


text:
import Base: isnan
import DataArrays: DataArray, isna, dropna, NA, NAtype, array
import DataFrames: DataFrame, ncol, convert
export GLRM

# TODO: identify categoricals automatically from PooledDataArray columns

function GLRM(df::DataFrame, k::Int;
              losses = Loss[], rx = QuadReg(.01), ry = QuadReg(.01),
              offset = true, scale = false,
              prob_scale = true, NaNs_to_NAs = true)
    if NaNs_to_NAs
        df = copy(df)
        NaNs_to_NAs!(df)
    end
    if losses == Loss[] # if losses not specified, identify ordinal, boolean and real columns
        reals, real_losses = get_reals(df)
        bools, bool_losses = get_bools(df)
        ordinals, ordinal_losses = get_ordinals(df)
        A = [df[reals] df[bools] df[ordinals]]
        labels = [names(df)[reals]; names(df)[bools]; names(df)[ordinals]]
        losses = [real_losses; bool_losses; ordinal_losses]
    else # otherwise require one loss function per column
        A = df
        ncol(df)==length(losses) ? labels = names(df) : error("please input one loss per column of dataframe")
    end
    # identify which entries in data frame have been observed (ie are not N/A)
    obs = observations(A)
    # initialize X and Y
    X = randn(k,size(A,1))
    Y = randn(k,size(A,2))
    # form model
    glrm = GLRM(A, losses, rx, ry, k, obs=obs, X=X, Y=Y, offset=offset, scale=scale)

    # scale model so it really computes the MAP estimator of the parameters
    if prob_scale
        prob_scale!(glrm)
    end
    return glrm, labels
end

function get_reals(df::DataFrame)
    m,n = size(df)
    reals = [typeof(df[i])<:AbstractArray{Float64,1} for i in 1:n]
    n1 = sum(reals)
    losses = Array(Loss,n1)
    for i=1:n1
        losses[i] = QuadLoss()
    end
    return reals, losses
end

function get_bools(df::DataFrame)
    m,n = size(df)
    bools = [isa(df[i], AbstractArray{Bool,1}) for i in 1:n]
    n1 = sum(bools)
    losses = Array(Loss,n1)
    for i=1:n1
        losses[i] = HingeLoss()
    end
    return bools, losses
end

function get_ordinals(df::DataFrame)
    m,n = size(df)
    # there must be a better way to check types...
    ordinals = [(isa(df[i], AbstractArray{Int,1}) ||
                 isa(df[i], AbstractArray{Int32,1}) ||
                 isa(df[i], AbstractArray{Int64,1})) for i in 1:n]
    nord = sum(ordinals)
    ord_idx = (1:size(df,2))[ordinals]
    maxs = zeros(nord,1)
    mins = zeros(nord,1)
    for i in 1:nord
        col = df[ord_idx[i]]
        try
            maxs[i] = maximum(dropna(col))
            mins[i] = minimum(dropna(col))
        end
    end

    # set losses and regularizers
    losses = Array(Loss,nord)
    for i=1:nord
        losses[i] = OrdinalHinge(mins[i],maxs[i])
    end
    return ordinals, losses
end


text:
module RocketTimerObservableTest

using Test
using Rocket

include("../test_helpers.jl")

@testset "TimerObservable" begin

    println("Testing: timer")

    @testset begin
        @test timer(100, 100) == Rocket.TimerObservable(100, 100)
    end

    @testset begin
        source = timer(42, 13)

        io = IOBuffer()

        show(io, source)

        printed = String(take!(io))

        @test occursin("TimerObservable", printed)
        @test occursin("42", printed)
        @test occursin("13", printed)

        subscription = subscribe!(source, void())

        show(io, subscription)

        printed = String(take!(io))

        @test occursin("TimerSubscription", printed)

        unsubscribe!(subscription)
    end

    run_testset([
        (
            source = timer(100, 30) |> take(3),
            values = @ts(100 ~ [ 0 ] ~ 30 ~ [ 1 ] ~ 30 ~ [ 2, c ] )
        ),
        (
            source = timer(100),
            values = @ts(100 ~ [ 0, c ])
        )
    ])

end

end


text:
using CSV

function get_compression(fh)
    if endswith(fh, "gz") return "gzip"
    elseif endswith(fh, "bz2") return "bz2"
    else return nothing end
end

function parse_sumstats(fh; alleles=false, dropna=true)
    @info "Reading from:" fh
    x = CSV.File(fh)

    if dropna x = x.dropna(how="any") end # TODO double check Julian equivalent

    @info x
    return x
end

#= TODO
function series_eq(x, y) end
function read_csv(fh, **kwargs) end
function sub_chr(s, chrom) end
function get_present_chrs(fh, num) end
function which_compression(fh) end
function read_cts(fh, match_snps) end
function ldscore_fromlist(flist, num=nothing) end
function l2_parser(fh, compression) end
function annot_parser(
    fh, compression, frqfile_full=nothing, compression_frq=nothing,
) end
function frq_parser(fh, compression) end
function ldscore(fh, num=nothing) end
function M(fh; num=nothing, N=2, common=False) end
function M_fromlist(flist; num=nothing, N=2, common=False) end
function annot(fh_list; num=nothing, frqfile=nothing) end
=#

#= TODO
function __ID_List_Factory__(
    colnames, keepcol, fname_end; header=nothing, usecols=nothing,
) end
    struct IDContainer end
        function __init__(container::IDContainer, fname) end
        function __read__(container::IDContainer, fname) end
        function loj(container::IDContainer, externalDf) end
=#


text:
using AbstractPlotting: PointTrans, xyz_boundingbox

@testset "Basic transforms" begin
    function point2(x::Point2)
        return Point2f0(x[1] + 10, x[2] - 77)
    end

    function point3(x::Point3)
        return Point3f0(x[1] + 10, x[2] - 77, x[3] /  4)
    end
    trans2 = PointTrans{2}(point2)
    trans3 = PointTrans{3}(point3)
    points = [Point2f0(0, 0), Point2f0(0, 1)]
    bb = xyz_boundingbox(trans2, points)
    @test bb == Rect(Vec3f0(10, -77, 0), Vec3f0(0, 1, 0))
    bb = xyz_boundingbox(trans3, points)
    @test bb == Rect(Vec3f0(10, -77, 0), Vec3f0(0, 1, 0))

    points = [Point3f0(0, 0, 4), Point3f0(0, 1, -8)]
    bb = xyz_boundingbox(trans2, points)
    @test bb == Rect(Vec3f0(10, -77, -8), Vec3f0(0, 1, 12))
    bb = xyz_boundingbox(trans3, points)
    @test bb == Rect(Vec3f0(10, -77, -2.0), Vec3f0(0, 1, 3.0))
end


text:
# Function to fit data using easyfit
"""
function to fit (using exponentials curves) time-correlation functions using EasyFit.
fitting("Directory", "file.dat", 2, 1)

path_to_file : Directory of the .dat file
file : name of the .dat file
comp : componente that will be analised. 2 for anion and 3 for cation / 2 - anion1; 3 - anion2; 4- cation
order : number of exponentials

"""
function fitting(path_to_file::String, file::String, comp::Int64 , order::Int64)
      archive  = readdlm("$path_to_file/$file", comments=true, comment_char='#')
      time     = archive[:,1]  # dt values
      data     = archive[:,comp]
      # fittings
      fit      = fitexp(time, data, n=order);
      return fit
end

##EasyFit
#function adjust(il::String, c)
#    file1 = readdlm("DATA_TIMECORR_4A/timecorr-$il-$c-ions.dat", comments=true, comment_char='#')
#    # time for each concentration
#    time1 = file1[1:500, 1]
#    # data for each system compoent
#    if length(il) != 10
#      data1 = file1[1:500, 2] # anion
#      data2 = file1[1:500, 3] # cation
#    else
#      data1 = file1[1:500, 2]  # anion 1
#      data2 = file1[1:500, 3]  # anion 2
#      data3 = file1[1:500, 4]  # cation
#    end
#
#    if length(il) != 10
#      println("Exponetials fits")
#      println("*************************************")
#      println("$il")
#      println("       ")
#      println("Cation : $(il[1:4])")
#      fit2 = fitexp(time1, data2, n=3)
#      println(fit2)
#      println("  ")
#      println("Anion  : $(il[5:7])")
#      fit1 = fitexp(time1, data1, n=3)
#      println(fit1)
#      println("*************************************")
#    end
#  end
#end
#
#
#
#
#
#
#
#
##using Polynomials
##using LsqFit
##using Plots
#
#```
#Function to perform a biexponential fit
#
#```
#function fitting(xdata, ydata)
#  p0 = [0.5, 0.5, 0.5, 0.5, 0.5, 0.5]   # guess points -
#  model(t, p) = p[1] * exp.(-p[2] * t) +  p[3] * exp.(-p[4] * t)  + p[5] * exp.(-p[6] * t)  # model trying to fit
#  nlin_fit(model, xdata, ydata, p0)
#end
#
#function fitting_md2(xdata, ydata)
#  p0 = [0.5, 0.5, 0.5]   # guess points -
#  model(t, p) = (1 - p[1]) * exp(-p[2] * t) + p[1] * exp(-p[3] * t)
#  nlin_fit(model, xdata, ydata, p0)
#end
#
#function nlin_fit(model, xdata, ydata, p0)
#
#  nlinfit = curve_fit(model, xdata, ydata, p0)
#  pfit = nlinfit.param
#  println(pfit)
#  xlin = range(xdata[1], xdata[end], length=200)
#
##  scatter(xdata, ydata, markersize=3, legend=:topright, label="data")
##  plot!(xlin, model(xlin, [p0[1], p0[2], p0[3], p0[4]]), label="initial model")
##  plot!(xlin, model(xlin, [pfit[1], pfit[2], pfit[3], pfit[4]]), linestyle=:dash, label="fitted model", dpi=200)
##
##  xaxis!("x")
##  yaxis!("y")
##  title!("nonlinear fit")
##
##  savefig("fitting_example.png")
#
#end
#
#


text:
const _small_numbers = String[
    "zero", "one", "two", "three", "four", "five", "six", "seven", "eight", "nine", "ten",
    "eleven", "twelve", "thirteen", "fourteen", "fifteen", "sixteen", "seventeen", "eighteen", "nineteen"]

const _small_number_dictionary = Dict{Char, String}(
    '0' => "zero", '1' => "one", '2' => "two", '3' => "three", '4' => "four", '5' => "five", '6' => "six",
    '7' => "seven", '8' => "eight", '9' => "nine")

const _tens = String[
    "twenty", "thirty", "forty", "fifty", "sixty", "seventy", "eighty", "ninety"]

# The scale_numbers are US, Canada and modern British (short scale)
const _scale_modern = String[
    "thousand",             "million",              "billion",                "trillion",               "quadrillion",
    "quintillion",          "sextillion",           "septillion",             "octillion",              "nonillion",
    "decillion",            "undecillion",          "duodecillion",           "tredecillion",           "quattuordecillion",
    "quindecillion",        "sexdecillion",         "septendecillion",        "octodecillion",          "novemdecillion",
    "vigintillion",         "unvigintillion",       "duovigintillion",        "tresvigintillion",       "quattuorvigintillion",
    "quinvigintillion",     "sesvigintillion",      "septemvigintillion",     "octovigintillion",       "novemvigintillion",
    "trigintillion",        "untrigintillion",      "duotrigintillion",       "trestrigintillion",      "quattuortrigintillion",
    "quintrigintillion",    "sestrigintillion",     "septentrigintillion",    "octotrigintillion",      "noventrigintillion",
    "quadragintillion"
    ]

const _scale_traditional_british = String[
    "thousand",              "million",               "thousand million",          "billion",                   "thousand billion",
    "trillion",              "thousand trillion",     "quadrillion",               "thousand quadrillion",      "quintillion",
    "thousand quintillion",  "sextillion",            "thousand sextillion",       "septillion",                "thousand septillion",    "octillion",             "thousand octillion",    "nonillion",                 "thousand nonillion",        "decillion",
    "thousand decillion",    "undecillion",           "thousand undecillion",      "duodecillion",              "thousand duodecillion",
    "tredecillion",          "thousand tredecillion", "quattuordecillion",         "thousand quattuordecillion","quindecillion",
    "thousand quindecillion","sedecillion",           "thousand sedecillion",      "septendecillion",           "thousand septendecillion",
    "octodecillion",         "thousand octodecillion","novendecillion",            "thousand novendecillion",   "vigintillion",
    "thousand vigintillion"
    ]

const _scale_traditional_european = String[
    "thousand",              "million",               "milliard",          "billion",              "billiard",
    "trillion",              "trilliard",             "quadrillion",       "quadrilliard",         "quintillion",
    "quintilliard",          "sextillion",            "sextilliard",       "septillion",           "septilliard",
    "octillion",             "octilliard",            "nonillion",         "nonilliard",           "decillion",
    "decilliard",            "undecillion",           "undecilliard",      "duodecillion",         "duodecilliard",
    "tredecillion",          "tredecilliard",         "quattuordecillion", "quattuordecilliard",   "quindecillion",
    "quindecilliard",        "sedecillion",           "sedecilliard",      "septendecillion",      "septendecilliard",
    "octodecillion",         "octodecilliard",        "novendecillion",    "novendecilliard",      "vigintillion",
    "vigintilliard"
    ]

const limit = BigInt(big(10)^120)
const limit_str = "10^120"


text:
"""
Subsample a dataframe or matrix.

`df`: A dataframe of matrix (rows is number of observations).
`n::Integer`: Size of subset. If `n >= size(df, 1)`, then the original dataset
              is returned. Otherwise, a dataset of size `(n, size(df, 2)` is
              returned.

**NOTE**: requires `DataFrames.jl`
"""
function subsample(df, n::Integer; seed::Integer=nothing)
  seed == nothing || Random.seed!(seed)

  N = size(df, 1)
  if n >= size(df, 1)
    return df
  else
    return df[sample(1:N, n), :]
  end
end

"""
Partition `data::DataFrame` into control and treatment groups according to
`treatmentsym::Symbol` (whether missing(->control) or not(->treatment)),
for a given marker (`markersym::Symbol`).

**NOTE**: requires `DataFrames.jl`
"""
function partition(data, markersym::Symbol; treatmentsym::Symbol=:treatment,
                   log_response=false)
  marker = data[!, markersym]
  yC = marker[ismissing.(data[!, treatmentsym])]
  yT = marker[.!ismissing.(data[!, treatmentsym])]
  if log_response
    yC .= log.(yC)
    yT .= log.(yT)
  end
  return (yC=yC, yT=yT)
end


text:
#=
Character and string classification functions

Copyright 2017-2018 Gandalf Software, Inc., Scott P. Jones,
Licensed under MIT License, see LICENSE.md
=#

# Recommended by deprecate
@static if V6_COMPAT
    text_width(str::AbstractString) = strwidth(str)
    text_width(ch::Char) = charwidth(ch)

    import Base: is_assigned_char, normalize_string, isnumber
    Base.is_assigned_char(ch::Chr) = is_assigned(ch)
    Base.normalize_string(str::Str, opt::Symbol) = normalize(str, opt)
    Base.strwidth(str::Str) = text_width(str)
    Base.charwidth(ch::Chr) = text_width(ch)
    isnumber(val::Chr) = is_numeric(val)
else
    Base.Unicode.normalize(str::Str, opt::Symbol) = normalize(str, opt)
    Base.Unicode.isassigned(ch::Chr) = is_assigned(ch)
    is_graphic(ch::Char) = is_graphic(codepoint(ch))
    is_alphanumeric(ch::Char) = is_alphanumeric(codepoint(ch))
end

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

## character column width function ##

text_width(ch::UInt8)    = Int(ifelse(ch < 0x7f, ch > 0x1f, ch > 0x9f & ch != 0xad))
text_width(ch::UInt16)   = utf8proc_charwidth(ch)
text_width(ch::UInt32)   = utf8proc_charwidth(ch)
text_width(ch::Chr)      = text_width(codepoint(ch))
text_width(ch::ASCIIChr) = Int(32 <= codepoint(ch) <= 126)

text_width(str::Str) = mapreduce(text_width, +, 0, str)
text_width(str::Str{Union{ASCIICSE,Latin_CSEs}}) = length(str)

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

@inline _cat(ch::CodeUnitTypes) = ch <= 0x10ffff ? utf8proc_cat(ch) : Cint(30)
@inline _cat(ch::Chr)     = _cat(codepoint(ch))

# returns code in 0:30 giving Unicode category
@inline category_code(ch::Union{Chr,CodeUnitTypes}) = _cat(ch)

# more human-readable representations of the category code
@inline category_abbrev(ch::CodeUnitTypes) = ch <= 0x10ffff ? utf8proc_cat_abbr(ch) : "In"
@inline category_abbrev(ch::Chr)           = category_abbrev(codepoint(ch))

category_string(ch::CodeUnitTypes) = category_strings[_cat(ch) + 1]
category_string(ch::Chr)           = category_string(codepoint(ch))

is_assigned(ch::CodeUnitTypes) = category_code(ch) != Uni.CN
is_assigned(ch::Chr)           = is_assigned(codepoint(ch))

_cat_mask(a) = a
@inline _cat_mask(a, b) = (1%UInt << a%UInt) | (1%UInt << b%UInt)
@inline _cat_mask(rng::UnitRange) =
    ((2%UInt << rng.stop%UInt) - 1) & ~((1%UInt << rng.start%UInt) - 1)

@inline _check_mask(ch, mask) = ((1%UInt << _cat(ch)%UInt) & mask) != 0

## libc character class predicates ##

# 0xb5, 0xdf, and 0xff cannot be uppercased in LatinStr, although they are lowercase
@inline _can_upper_l(c) = (0xe0 <= c <= 0xfe) & (c != 0xf7)
@inline _can_upper(c) = _islower_a(c) | _can_upper_l(c)

@inline _iscntrl(ch) = (ch <= 0x1f) | (0x7f <= ch <= 0x9f)
@inline _isdigit(ch) = (ch - '0'%UInt8) <= 9
@inline _isxdigit(ch) = _isdigit(ch) | (ch - 'A'%UInt8 < 6) | (ch - 'a'%UInt8 < 6)

const _isupper_mask   = _cat_mask(Uni.LU, Uni.LT)
const _isalpha_mask   = _cat_mask(Uni.LU : Uni.LO)
const _isnumeric_mask = _cat_mask(Uni.ND : Uni.NO)
const _ispunct_mask   = _cat_mask(Uni.PC : Uni.PO)
const _isprint_mask   = _cat_mask(Uni.LU : Uni.ZS)
const _isgraph_mask   = _cat_mask(Uni.LU : Uni.SO)
const _isalnum_mask   = _isnumeric_mask | _isalpha_mask

############################################################################
# Definitions for characters in the ASCII subset of Unicode

const _isnumeric_a = _isdigit
@inline _ispunct_a(ch) = ((1%UInt128 << ch) & 0x2800_0000_b800_0001_8c00_f7ee_0000_0000) != 0
@inline _isspace_a(ch) = (ch == 32) | (9 <= ch <= 13)
@inline _islower_a(ch) = (ch - 'a'%UInt8) < 26
@inline _isupper_a(ch) = (ch - 'A'%UInt8) < 26
@inline _isalpha_a(ch) = _islower_a(ch) | _isupper_a(ch)
@inline _isalnum_a(ch) = _isdigit(ch) | _isalpha_a(ch)
@inline _isprint_a(ch) = 0x20 <= ch < 0x7f
@inline _isgraph_a(ch) = 0x20 < ch < 0x7f

############################################################################
# Definitions for characters in the Latin1 subset of Unicode, but not in the ASCII subset

@inline _isnumeric_l(ch) = (ch <= 0xbe && ((1<<(ch-0xb2)) & 0x1c83) != 0)
@inline _ispunct_l(ch)   = ((1%UInt64 << (ch-0x80)) & 0x88c0_0882_0000_0000) != 0
@inline _isspace_l(ch)   = (ch == 0x85) | (ch == 0xa0)
@inline _islower_l(c)    = ((0xdf <= c <= 0xff) & (c != 0xf7)) | (c == 0xb5)
@inline _isupper_l(c)    = (0xc0 <= c%UInt8 <= 0xde) & (c != 0xd7)
@inline _isalpha_l(c)    = ((0xc0 <= c <= 0xff) & (c != 0xf7) & (c != 0xd7)) | (c == 0xb5)
@inline _isalnum_l(c)    = _isalpha_l(c) || _isnumeric_l(c)
@inline _isprint_l(ch)   = ((0xa0 <= ch <= 0xff) & (ch != 0xad))
@inline _isgraph_l(ch)   = ((0xa0 < ch <= 0xff) & (ch != 0xad))

############################################################################
# Definitions for any Unicode codepoint (requires call to utf8proc) (only used for non-Latin1)

@inline _isnumeric_u(ch) = _check_mask(ch, _isnumeric_mask)
@inline _ispunct_u(ch)   = _check_mask(ch, _ispunct_mask)
@inline _isspace_u(ch)   = _cat(ch) == Uni.ZS
@inline _islower_u(ch)   = _cat(ch) == Uni.LL
@inline _isupper_u(ch)   = _check_mask(ch, _isupper_mask)
@inline _isalpha_u(ch)   = _check_mask(ch, _isalpha_mask)
@inline _isalnum_u(ch)   = _check_mask(ch, _isalnum_mask)
@inline _isprint_u(ch)   = _check_mask(ch, _isprint_mask)
@inline _isgraph_u(ch)   = _check_mask(ch, _isgraph_mask)

############################################################################
# Fallback definitions for all Chr types

@inline is_control(ch::CodeUnitTypes)   = _iscntrl(ch)
@inline is_digit(ch::CodeUnitTypes)     = _isdigit(ch)
@inline is_hex_digit(ch::CodeUnitTypes) = _isxdigit(ch)

@inline is_control(ch::Chr)   = is_control(codepoint(ch))
@inline is_digit(ch::Chr)     = is_digit(codepoint(ch))
@inline is_hex_digit(ch::Chr) = is_hex_digit(codepoint(ch))

@inline is_ascii(ch::Chr)          = is_ascii(codepoint(ch))
@inline is_ascii(ch::Unsigned)     = ch <= 0x7f
@inline is_ascii(ch::ASCIIChr)     = true

@inline is_latin(ch::AbstractChar) = is_latin(codepoint(ch))
@inline is_latin(ch::Unsigned)     = ch <= 0xff

@inline is_bmp(ch::AbstractChar)   = is_bmp(codepoint(ch))
@inline is_bmp(ch::Unsigned)       = ch <= 0xffff && !is_surrogate_codeunit(ch)
@inline is_bmp(ch::UInt8)          = true

@inline is_unicode(ch::AbstractChar) = is_unicode(codepoint(ch))
@inline is_unicode(ch::Unsigned)     = ch <= 0x10ffff && !is_surrogate_codeunit(ch)
@inline is_unicode(ch::UInt8)        = true

const _catfuns =
    ((:numeric,      :numeric),
     (:punctuation,  :punct),
     (:space,        :space),
     (:lowercase,    :lower),
     (:uppercase,    :upper),
     (:alpha,        :alpha),
     (:alphanumeric, :alnum),
     (:printable,    :print),
     (:graphic,      :graph))

for (nnam, fnam) in _catfuns
    isnam  = Symbol(string("is_", nnam))
    namroot  = string("_is", fnam)
    fnam_a  = Symbol(string(namroot, "_a"))
    fnam_al = Symbol(string(namroot, "_al"))
    fnam_ch = Symbol(string(namroot, "_ch"))

    @eval $(fnam_al)(ch) = is_ascii(ch) ? $(fnam_a)(ch) : $(Symbol(string(namroot, "_l")))(ch)
    @eval $(fnam_ch)(ch) = is_latin(ch) ? $(fnam_al)(ch) : $(Symbol(string(namroot, "_u")))(ch)
    @eval $(isnam)(ch::CodeUnitTypes)  = $(fnam_ch)(ch)
    @eval $(isnam)(ch::Chr)            = $(fnam_ch)(codepoint(ch))
    @eval $(isnam)(ch::ASCIIChr)       = $(fnam_a)(codepoint(ch))
    @eval $(isnam)(ch::LatinChars)     = $(fnam_al)(codepoint(ch))
end

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

@static if isdefined(Base, :ismalformed)
    Base.ismalformed(ch::Chr) = false
    Base.isoverlong(ch::Chr) = false
    is_malformed(ch) = ismalformed(ch)
    is_overlong(ch) = isoverlong(ch)
else
    is_malformed(ch) = false
    is_overlong(ch) = false
end

function is_latin(str::MaybeSub{String})
    (siz = sizeof(str)) == 0 && return true
    @preserve str begin
        pnt = pointer(str)
        fin = pnt + siz
        while pnt < fin
            cu = get_codeunit(pnt)
            # cu must be 1) 0-0x7f, or 2) 0xc2 or 0xc3 followed by 0x80-0xbf
            (cu < 0x7f ||
             ((cu - 0xc2) < 0x02 &&
              (pnt += 1) < fin && is_valid_continuation(get_codeunit(pnt)))) ||
              return false
            pnt += 1
        end
        true
    end
end

@inline function check_3byte(cu, pnt)
    b2 = get_codeunit(pnt-1)
    b3 = get_codeunit(pnt)
    is_valid_continuation(b2) && is_valid_continuation(b3) &&
        !is_surrogate_codeunit(((cu & 0x0f)%UInt32 << 12) | ((b2 & 0x3f)%UInt32 << 6) | (b3 & 0x3f))
end

function is_bmp(str::MaybeSub{String})
    (siz = sizeof(str)) == 0 && return true
    @preserve str begin
        pnt = pointer(str)
        fin = pnt + siz
        while pnt < fin
            cu = get_codeunit(pnt)
            # cu must be 1) 0-0x7f, or 2) 0xc2 or 0xc3 followed by 0x80-0xbf
            # c2-df -> de,df
            (cu < 0x7f ||
             ((cu - 0xc2) < 0x1e && (pnt += 1) < fin && checkcont(pnt)) ||
             ((cu - 0xe0) < 0x0f && (pnt += 2) < fin && check_3byte(cu, pnt))) ||
             return false
             pnt += 1
        end
        true
    end
end

function is_latin(str::AbstractString)
    @inbounds for ch in str
        is_latin(ch) || return false
    end
    true
end

function is_bmp(str::AbstractString)
    @inbounds for ch in str
        is_bmp(ch) || return false
    end
    true
end

function is_unicode(str::AbstractString)
    @inbounds for ch in str
        is_unicode(ch) || return false
    end
    true
end


text:
@testset "trace translators" begin

@testset "DeterministicTraceTranslator" begin

    @gen function p1()
        r ~ inv_gamma(1, 1)
        theta ~ uniform(-pi/2, pi/2)
    end

    @gen function p2()
        x ~ normal(0, 1)
        y ~ normal(0, 1)
    end

    @transform f (t1) to (t2) begin
        r = @read(t1[:r], :continuous)
        theta = @read(t1[:theta], :continuous)
        @write(t2[:x], r * cos(theta), :continuous)
        @write(t2[:y], r * sin(theta), :continuous)
    end

    @transform finv (t2) to (t1) begin
        x = @read(t2[:x], :continuous)
        y = @read(t2[:y], :continuous)
        r = sqrt(x^2 + y^2)
        @write(t1[:r], sqrt(x^2 + y^2), :continuous)
        @write(t1[:theta], atan(y, x), :continuous)
    end

    pair_bijections!(f, finv)

    translator = DeterministicTraceTranslator(p2, (), choicemap(), f)
    t1, _ = generate(p1, (), choicemap(:theta => 0))
    t2, log_weight = translator(t1; check=true)
    @test t2[:y] == 0

    translator = DeterministicTraceTranslator(p1, (), choicemap(), finv)
    t2, _ = generate(p2, (), choicemap(:y => 0, :x => 1))
    t1, log_weight = translator(t2; check=true)
    @test t1[:theta] == 0

end

@testset "SimpleExtendingTraceTranslator" begin

    @gen function model(T::Int)
        for t in 1:T
            z = {(:z, t)} ~ normal(0, 1)
            x = {(:x, t)} ~ normal(z, 1)
        end
    end

    @gen function proposal(trace::Trace, x)
        t = get_args(trace)[1] + 1
        {(:z, t)} ~ normal(x, 1)
    end

    translator = SimpleExtendingTraceTranslator(
        p_new_args=(2,), p_argdiffs=(UnknownChange(),),
        new_observations=choicemap((:x, 2) => 5.0),
        q_forward=proposal, q_forward_args=(5.0,))
    t1  = simulate(model, (1,))
    t2, log_weight = translator(t1)

    prop_choices = choicemap((:z, 2) => t2[(:z, 2)])
    prop_weight, _ = assess(proposal, (t1, 5.0), prop_choices)
    constraints = merge(prop_choices, choicemap((:x, 2) => 5.0))
    t3, up_weight, _, _ = update(t1, (2,), (UnknownChange(),), constraints)
    @test log_weight == up_weight - prop_weight

end

@testset "SymmetricTraceTranslator" begin

    @gen function model()
        z ~ bernoulli(0.5)
        if z
            i ~ uniform_discrete(1, 10)
        else
            x ~ uniform(0, 1)
        end
    end

    @gen function proposal(trace)
        if trace[:z]
            dx ~ uniform(0.0, 0.1)
        end
    end

    @transform involution (p1_trace, q1_trace) to (p2_trace, q2_trace) begin
        if @read(p1_trace[:z], :discrete)
            @write(p2_trace[:z], false, :discrete)
            i = @read(p1_trace[:i], :discrete)
            dx = @read(q1_trace[:dx], :continuous)
            x = (i-1)/10 + dx
            @write(p2_trace[:x], x, :continuous)
        else
            @write(p2_trace[:z], true, :discrete)
            x = @read(p1_trace[:x], :continuous)
            i = ceil(x * 10)
            @write(p2_trace[:i], i, :discrete)
            @write(q2_trace[:dx], x - (i-1)/10, :continuous)
        end
    end

    is_involution!(involution)

    translator = SymmetricTraceTranslator(proposal, (), involution)

    t1, _ = generate(model, (), choicemap(:z => false, :x => 0.95))
    t2, log_weight = translator(t1; check=true)
    @test t2[:z] == true && t2[:i] == 10

end

@testset "GeneralTraceTranslator" begin

    @gen function p1()
        x ~ uniform(0, 1)
        y ~ uniform(0, 1)
    end

    @gen function p2()
        i ~ uniform_discrete(1, 10) # interval [(i-1)/10, i/10]
        j ~ uniform_discrete(1, 10) # interval [(j-1)/10, j/10]
    end

    @gen function q1(p1_trace) end

    @gen function q2(p2_trace)
        dx ~ uniform(0.0, 0.1)
        dy ~ uniform(0.0, 0.1)
    end

    @transform f (p1_trace, q1_trace) to (p2_trace, q2_trace) begin
        x = @read(p1_trace[:x], :continuous)
        y = @read(p1_trace[:y], :continuous)
        i = ceil(x * 10)
        j = ceil(y * 10)
        @write(p2_trace[:i], i, :discrete)
        @write(p2_trace[:j], j, :discrete)
        @write(q2_trace[:dx], x - (i-1)/10, :continuous)
        @write(q2_trace[:dy], y - (j-1)/10, :continuous)
    end

    @transform f_inv (p2_trace, q2_trace) to (p1_trace, q1_trace) begin
        i = @read(p2_trace[:i], :discrete)
        j = @read(p2_trace[:j], :discrete)
        dx = @read(q2_trace[:dx], :continuous)
        dy = @read(q2_trace[:dy], :continuous)
        x = (i-1)/10 + dx
        y = (j-1)/10 + dy
        @write(p1_trace[:x], x, :continuous)
        @write(p1_trace[:y], y, :continuous)
    end

    pair_bijections!(f, f_inv)

    translator = GeneralTraceTranslator(
        p_new=p2, p_new_args=(),
        new_observations=choicemap(),
        q_forward=q1, q_forward_args=(),
        q_backward=q2, q_backward_args=(), f=f)

    t1, _ = generate(p1, (), choicemap(:x => 0.05, :y => 0.95))
    t2, log_weight = translator(t1; check=true)
    @test t2[:i] == 1 && t2[:j] == 10

end

end