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

__precompile__()
module JEMSS

# animation
using HTTP
using Sockets
using WebSockets
using JSON

# files
using DelimitedFiles
using LightXML
using Serialization

# misc
using Base64
using CRC32c
using Mmap
using Pkg
using Printf
using Random
import Random: AbstractRNG, GLOBAL_RNG, MersenneTwister
using SparseArrays
using LightGraphs
using ArchGDAL
using Hungarian
using OffsetArrays
using Base.Iterators: Stateful
using Parameters
using Base.Threads

# optimisation
using JuMP
using Cbc
using GLPK
using GLPKMathProgInterface
if haskey(Pkg.installed(), "Gurobi") && haskey(ENV, "GUROBI_HOME") using Gurobi end

# statistics
using Distributions
using HypothesisTests
using StatsBase
using StatsFuns

# simulation functions
export
	initSim, runConfig, # run_config
	simulate!, simulateToTime!, simulateToEnd!, backup!, reset!, simulateNextEvent!, # simulation
	setSimReps!, simulateRep!, simulateReps!, resetRep!, resetReps!, makeRepsRunnable!, # replication
	animate!, animate, # animation
	backupSim!, resetSim! # compat

# file functions
export
	readDlmFileNextLine!, readDlmFile, openNewFile, writeDlmLine!, arrayDict, writeTablesToFile!, writeTablesToFile, readTablesFromFile, readTablesFromData, tableRowsFieldDicts, fileChecksum, serializeToFile, deserializeFile, joinPathIfNotAbs, interpolateString, xmlFileRoot, findElt, eltContent, eltContentVal, eltContentInterpVal, childrenNodeNames, selectXmlFile, # file_io
	readGenConfig, runGenConfig, runGenConfigCalls, makeCalls, # gen_sim_files
	readAmbsFile, readArcsFile, readCallsFile, readCompTableFile, readDemandFile, readDemandCoverageFile, readEventsFile, readGeoFile, readHospitalsFile, readMapFile, readMobilisationDelayFile, readMultiCompTableFile, readNodesFile, readPrioritiesFile, readPriorityListFile, readPriorityListsFile, readRasterFile, readRedispatchFile, readRNetTravelsFile, readStationsFile, readStatsControlFile, readStatsDictFile, readTravelFile, readDeploymentsFile, readZhangIpParamsFile, # read_sim_files
	writeAmbsFile, writeArcsFile, writeCallsFile, writeDemandFile, writeDemandCoverageFile, writeHospitalsFile, writeMapFile, writeMultiCompTableFile, writeNodesFile, writePrioritiesFile, writePriorityListFile, writePriorityListsFile, writeRasterFile, writeRedispatchFile, writeRNetTravelsFile, writeStationsFile, writeTravelFile, writeZhangIpParamsFile, openOutputFiles!, writeOutputFiles, writeMiscOutputFiles, closeOutputFiles!, writeEventToFile!, writeDeploymentsFile, writeBatchMeanResponseDurationsFile, # write_sim_files
	writeStatsFiles, writeAmbsStatsFile, writeCallsStatsFile, writeHospitalsStatsFile, writeStationsStatsFile, writeStatsDictFile # write_sim_files - stats

# move up
export
	initCompTable!, initDmexclp!, initMultiCompTable!, initPriorityList!, initZhangIp!, initTemp0!, initTemp1!, initTemp2!,
	setMoveUpModule!

# misc functions
export
	isBusy, isFree, isWorking, isGoingToStation, isTravelling, # ambulance
	setSimCalls!, # call
	calcCoverBound!, initCoverBound, approxDistr, convolute, makeCdf, binarySearch, calcAmbBusyDurationProbUpperBounds, calcAmbBusyDurationLowerBoundDistrs!, calcNumAmbsMaxCoverageFrac, calcQueuedDurationsMaxCoverageFrac, initCoverBoundSim, simulateCoverBound!, simulateCoverBoundLowerBound!, writeAmbBusyDurationProbUpperBoundsFile, readAmbBusyDurationProbUpperBoundsFile, writeQueuedDurationsMaxCoverageFracFile, readQueuedDurationsMaxCoverageFracFile, # cover bound
	initDemand!, initDemandCoverage!, # demand
	printEvent, # event
	findNearestNode, # grid
	squareDist, normDist, offRoadTravelTime, linearInterpLocation, randLocation, # location
	isFNodeInRGraph, shortestPathNextRNode, shortestPathNextRArc, shortestPathData, shortestPathTravelTime, shortestPathDistance, shortestPath, findRArcFromFNodeToFNode, # network
	rasterRandLocations, printRasterSize, # raster
	shortestRouteTravelTime!, # route
	getCallResponseDurations, getAvgCallResponseDuration, getCallsReachedInTime, countCallsReachedInTime, printSimStats, printAmbsStats, printCallsStats, printHospitalsStats, calcBatchMeans, calcBatchMeanResponseDurations, meanErrorPlot, calcAR0DurbinWatsonTestPValue, tDistrHalfWidth, confInterval, getPeriodStatsList, getRepsPeriodStatsList, getRepPeriodStats, statsDictFromPeriodStatsList, # statistics
	checkCompTable, isCompTableNested, nestCompTable, unnestCompTable, makeRandNestedCompTable, # compliance table
	checkMultiCompTable, # multi-compliance table
	makeRandDeployment, makeRandDeployments, deploymentToStationsNumAmbs, stationsNumAmbsToDeployment, getDeployment, getStationsNumAmbs, setAmbStation!, applyDeployment!, applyStationsNumAmbs!, simulateDeployment!, simulateDeployments!, # deployment
	checkPriorityList, makeRandPriorityList, # priority list
	solveMclp!, solveMclp, # mclp
	solveMexclp!, # mexclp
	solveNestedMclp!, solveNestedMclp, # nested mclp
	solvePMedian, # p-median
	flatten, # dict
	runParallel! # parallel

# types
export
	Location, Point, Node, Arc, Graph, NetTravel, Network, TravelMode, Travel, Route,
	Event, Ambulance, Call, Hospital, Station, Redispatch,
	Map, GridSearchRect, GridRect, Grid, Raster, RasterSampler, DemandMode, Demand, PointsCoverageMode, DemandCoverage,
	CoverBoundSimRep, CoverBoundSim, CoverBoundMode, CoverBound,
	CompTableData, DdsmData, DmexclpData, MultiCompTableData, PriorityListData, ZhangIpData, Temp0Data, Temp1Data, Temp2Data, MoveUpData,
	MeanAndHalfWidth, AmbulanceStats, CallStats, HospitalStats, StationStats, SimPeriodStats, SimStats,
	File, Table, Resimulation, Simulation,
	DistrRng, GenConfig, MobilisationDelay

# defs - consts
export
	Float, FloatSpTime, FloatSpDist, IntRNode, IntFadj, # type alias
	nullIndex, nullX, nullY, nullTime, nullDist, nullHist, NullHist, # nulls
	priorities, numPriorities, # priorities
	Deployment, CompTable, MultiCompTable, NestedCompTable, PriorityList

# defs - enums
export
	Priority, nullPriority, highPriority, medPriority, lowPriority, # priorities
	AmbClass, nullAmbClass, als, bls, # ambulance classes
	EventForm, nullEvent, ambGoesToSleep, ambWakesUp, callArrives, considerDispatch, ambDispatched, ambMobilised, ambReachesCall, ambGoesToHospital, ambReachesHospital, ambBecomesFree, ambReturnsToStation, ambReachesStation, considerMoveUp, ambMoveUpToStation,
	AmbStatus, ambNullStatus, ambSleeping, ambIdleAtStation, ambMobilising, ambGoingToCall, ambAtCall, ambGoingToHospital, ambAtHospital, ambFreeAfterCall, ambReturningToStation, ambMovingUpToStation,
	AmbStatusSet, ambWorking, ambBusy, ambFree, ambTravelling, ambGoingToStation,
	CallStatus, callNullStatus, callScreening, callQueued, callWaitingForAmb, callOnSceneTreatment, callGoingToHospital, callAtHospital, callProcessed,
	RouteStatus, routeNullStatus, routeBeforeStartNode, routeOnPath, routeAfterEndNode,
	MoveUpModule, nullMoveUpModule, compTableModule, dmexclpModule, multiCompTableModule, priorityListModule, zhangIpModule, temp0Module, temp1Module, temp2Module

# deprecated
export
	Depol, makeRandDeploymentPolicy, makeRandDeploymentPolicies, applyDeploymentPolicy!, simulateDeploymentPolicy!, simulateDeploymentPolicies!, writeDeploymentPoliciesFile, readDeploymentPoliciesFile, # renamed "deployment policy" to "deployment"
	isSameLocation, findNearestNodeInGrid

include("defs.jl")

include("misc/dict.jl")
include("misc/histogram.jl")
include("misc/parallel.jl")
include("misc/rand.jl")
include("misc/stream.jl")

include("types/types.jl")
include("types/ambulance.jl")
include("types/call.jl")
include("types/demand.jl")
include("types/event.jl")
include("types/graph.jl")
include("types/grid.jl")
include("types/location.jl")
include("types/network.jl")
include("types/raster.jl")
include("types/route.jl")
include("types/statistics.jl")
include("types/station.jl")
include("types/travel.jl")

include("file/file_io.jl")
include("file/gen_sim_files.jl")
include("file/read_sim_files.jl")
include("file/write_sim_files.jl")

include("decision/call.jl")
include("decision/dispatch.jl")

include("decision/move_up/comp_table.jl")
include("decision/move_up/ddsm.jl")
include("decision/move_up/dmexclp.jl")
include("decision/move_up/move_up_common.jl")
include("decision/move_up/multi_comp_table.jl")
include("decision/move_up/priority_list.jl")
include("decision/move_up/temp0.jl")
include("decision/move_up/temp1.jl")
include("decision/move_up/temp2.jl")
include("decision/move_up/zhang_ip.jl")

include("resimulation.jl")
include("run_config.jl")
include("simulation.jl")
include("replication.jl")
include("statistics.jl")

include("animation/animation.jl")

include("misc/deployment.jl")

include("optim/mclp.jl")
include("optim/mexclp.jl")
include("optim/p_median.jl")
include("optim/cover_bound.jl")
include("optim/nested_mclp.jl")

include("deprecated.jl")

end


text:
using Neuroimaging, Lexicon

save("doc/API.md", Neuroimaging)

ext: jl


text:
function futuresdata(;
  inputpath::String = PARAM[:inputpath],
  futuresfile::String = PARAM[:capacity][:futuresfile],
  futuresfileext::String = PARAM[:capacity][:futuresfileext],
  )

  futures = CSV.read("$inputpath\\$futuresfile.$futuresfileext") |> DataFrame
  select!(futures, PARAM[:capacity][:futurescolumns])

  parsemfloat(f)=something(tryparse(Float64, f), missing)
  futures.test = parsemfloat.(futures.test)

  #display(describe(futures))
  return futures
end

#a function for pulling lots of tickers and automating common parameter choices
function multipull(S;livepull::Bool,
  tickers::Vector{String}=error("tickers kw are required"),
  fields::Vector{String}=error("fields kw are required"),
  startdate::Date=Date(1900,1,1),
  enddate::Date=today(),
  periodicity::String=error("periodicity kw are required"),
  adjustmentAbnormal::Bool = true,
  adjustmentNormal::Bool=	true,
  currency::String=	"USD",
  DateFormat::String = "D",) #is this right?

  #WARNING- what to do about dates?

  #set up a parameter dict
  multipullparams::Dict{Symbol, Any} = Dict(:periodicity=>periodicity,
    :adjustmentAbnormal=>adjustmentAbnormal,
    :adjustmentNormal=>adjustmentNormal,
    :currency=>currency,
    #:DateFormat=>DateFormat,
    )
  if Period=="DAILY"
    multipullparams[:CDR] = "7D" #are tehse right? where is teh correct list of parameters?
    multipullparams[:Quote] = "Close"
  end

  if livepull
    results = BLPData.bdh(S, tickers, fields, startdate, enddate;
      multipullparams...) |> DataFrame
    display(results)
  else
    #serialization- compressed?
  end

  results
end

#main entry point
function bbscript(S, ::Val{:futurespull}; livepull=error("Live pull is required"))
  futures::DataFrame = futuresdata()
  tickers::Vector{String}=futures.fullticker

  pull = multipull(S,
    livepull=livepull,
    tickers=["ES1 Index"],
    fields=["px_last"],
    periodicity="Yearly")

end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
# CodeGenerator
function GetModule(x::CodeGenerator)
    @check_ptrs x
    return LLVM.Module(clang_CodeGenerator_GetModule(x))
end

function ReleaseModule(x::CodeGenerator)
    @check_ptrs x
    return LLVM.Module(clang_CodeGenerator_ReleaseModule(x))
end

function GetDeclForMangledName(x::CodeGenerator, s::String)
    @check_ptrs x
    return Decl(clang_CodeGenerator_GetDeclForMangledName(x, s))
end

function StartModule(x::CodeGenerator, ctx::Context, mod_name::String)
    @check_ptrs x
    @assert ctx.ref != C_NULL "LLVMContextRef has a NULL pointer."
    return LLVM.Module(clang_CodeGenerator_StartModule(x, ctx, mod_name))
end


text:
# Internal facilities for fast random number generation

immutable RandIntSampler  # for generating Int samples in [0, K-1]
	a::Int
	Ku::Uint
	U::Uint

	RandIntSampler(K::Int) = (Ku = uint(K); new(1, Ku, div(typemax(Uint), Ku) * Ku))
	RandIntSampler(a::Int, b::Int) = (Ku = uint(b-a+1); new(a, Ku, div(typemax(Uint), Ku) * Ku))
end

function rand(s::RandIntSampler)
	x = rand(Uint)
	while x >= s.U
		x = rand(Uint)
	end
	s.a + int(rem(x, s.Ku))
end

randi(K::Int) = rand(RandIntSampler(K))
randi(a::Int, b::Int) = rand(RandIntSampler(a, b))

# draw a number from a binomial distribution

rand_binom(n::Real, p::Real) = int(ccall((:rbinom, :libRmath), Float64, (Float64, Float64), n, p))


text:
module DisplayInterface # PushInterface

# https://github.com/Ableton/push-interface/blob/master/doc/AbletonPush2MIDIDisplayInterface.asc#3-display-interface

const DisplayPixelWidth      = 960
const DisplayPixelHeight     = 160

const DataSourceWidth        = 1024                       #  960+64
const LineSize               = 2DataSourceWidth           #  2 * (960+64)

const LineCountPerSendBuffer = Int(DisplayPixelHeight/1)  #  160
const SendBufferSize         = LineCountPerSendBuffer * LineSize

include("DisplayInterface/Canvas.jl")
include("DisplayInterface/framebuffer.jl")

end # module PushInterface.DisplayInterface


text:
include("fundamental_matrix.jl")
include("homography_matrix.jl")


text:
using JavaCall
import Base: map, reduce, count, collect, close
import Base.Iterators

# config
const JSparkConf = @jimport org.apache.spark.SparkConf
# context
const JSparkContext = @jimport org.apache.spark.SparkContext
const JJavaSparkContext = @jimport org.apache.spark.api.java.JavaSparkContext
# SQL
const JSparkSession = @jimport org.apache.spark.sql.SparkSession
const JSparkSessionBuilder = @jimport org.apache.spark.sql.SparkSession$Builder
const JDataFrameReader = @jimport org.apache.spark.sql.DataFrameReader
const JDataFrameWriter = @jimport org.apache.spark.sql.DataFrameWriter
const JDataset = @jimport org.apache.spark.sql.Dataset
const JRelationalGroupedDataset = @jimport org.apache.spark.sql.RelationalGroupedDataset
const JGenericRow = @jimport org.apache.spark.sql.catalyst.expressions.GenericRowWithSchema
const JRowFactory = @jimport org.apache.spark.sql.RowFactory
const JRow = @jimport org.apache.spark.sql.Row
const JColumn = @jimport org.apache.spark.sql.Column
const JSQLFunctions = @jimport org.apache.spark.sql.functions
# RDD
const JRDD = @jimport org.apache.spark.rdd.RDD
const JJavaRDD = @jimport org.apache.spark.api.java.JavaRDD
const JJavaPairRDD = @jimport org.apache.spark.api.java.JavaPairRDD
const JJuliaRDD = @jimport org.apache.spark.api.julia.JuliaRDD
const JJuliaPairRDD = @jimport org.apache.spark.api.julia.JuliaPairRDD
# utils
const JRDDUtils = @jimport org.apache.spark.api.julia.RDDUtils
# Java utils
const JIterator = @jimport java.util.Iterator
const JList = @jimport java.util.List
const JMap = @jimport java.util.Map
const JArrayList = @jimport java.util.ArrayList
const JHashMap = @jimport java.util.HashMap

include("init.jl")
include("serialization.jl")
include("config.jl")
include("context.jl")
include("sql.jl")
include("rdd.jl")
include("attach.jl")
include("worker.jl")


text:
mutable struct TransferEntropy{D,E} <: InfoDist
    states::Vector{Int}
    histories::Vector{Int}
    sources::Vector{Int}
    predicates::Vector{Int}
    bs::NTuple{D,Int}
    bt::NTuple{E,Int}
    k::Int
    Bs::Int
    Bt::Int
    q::Int
    N::Int

    function TransferEntropy(bs::NTuple{D, <:Integer}, bt::NTuple{E, <:Integer}; k::Int=1) where {D,E}
        if D == 0 || E == 0
            throw(MethodError(TransferEntropy, (bs, ts, kwargs...)))
        end
        if any(b -> b < 2, bs) || any(b -> b < 2, bt)
            throw(ArgumentError("the support of each random variable must be at least 2"))
        end
        if k < 1
            throw(ArgumentError("history length must be at least 1"))
        end
        Bs, Bt = prod(bs), prod(bt)
        q = Bt^k
        states = zeros(Int, Bs*Bt*q)
        histories = zeros(Int, q)
        sources = zeros(Int, Bt*q)
        predicates = zeros(Int, Bt*q)
        new{D,E}(states, histories, sources, predicates, bs, bt, k, Bs, Bt, q, 0)
    end
end
TransferEntropy(bs::Integer, bt::Integer; kwargs...) = TransferEntropy((bs,), (bt,); kwargs...)
TransferEntropy(bs::Integer, bt::NTuple; kwargs...) = TransferEntropy((bs,), bt; kwargs...)
TransferEntropy(bs::NTuple, bt::Integer; kwargs...) = TransferEntropy(bs, (bt,); kwargs...)

function TransferEntropy(source::AbstractArray{Int,3}, target::AbstractArray{Int,3}; k::Int=1)
    if isempty(source) || isempty(target)
        throw(ArgumentError("arguments must not be empty"))
    end
    smax = max.(2, maximum(source, dims=(2,3)))
    tmax = max.(2, maximum(target, dims=(2,3)))
    observe!(TransferEntropy(tuple(smax...), tuple(tmax...); k), source, target)
end

function TransferEntropy(source::AbstractArray{Int,2}, target::AbstractArray{Int,2}; k::Int=1)
    if isempty(source) || isempty(target)
        throw(ArgumentError("arguments must not be empty"))
    end
    smax = max.(2, maximum(source, dims=2))
    tmax = max.(2, maximum(target, dims=2))
    observe!(TransferEntropy(tuple(smax...), tuple(tmax...); k), source, target)
end

function TransferEntropy(source::AbstractArray{Int,1}, target::AbstractArray{Int,1}; k::Int=1)
    if isempty(source) || isempty(target)
        throw(ArgumentError("arguments must not be empty"))
    end
    smax = max(2, maximum(source))
    tmax = max(2, maximum(target))
    observe!(TransferEntropy(tuple(smax...), tuple(tmax...); k), source, target)
end

@inline function clear!(dist::TransferEntropy)
    dist.states[:] .= 0
    dist.histories[:] .= 0
    dist.sources[:] .= 0
    dist.predicates[:] .= 0
    dist.N = 0
    dist
end

function observe!(dist::TransferEntropy, source::AbstractArray{Int,3}, target::AbstractArray{Int,3})
    if size(source, 3) != size(target, 3)
        throw(ArgumentError("time series should have the same number of replicates"))
    end

    @views for i in 1:size(source, 3)
        observe!(dist, source[:,:,i], target[:,:,i])
    end

    dist
end

function observe!(dist::TransferEntropy, source::AbstractArray{Int,2}, target::AbstractArray{Int,2})
    if size(source, 2) != size(target, 2)
        throw(ArgumentError("time series should have the same number of timesteps"))
    elseif length(target) <= dist.k
        throw(ArgumentError("target series is too short given k=$(dist.k)"))
    elseif any(b -> b < 1, source) || any(b -> b < 1, target)
        throw(ArgumentError("observations must be positive, nonzero"))
    end

    N = size(target, 2) - 1
    dist.N += N - dist.k + 1
    history = 0
    @views for t in 1:dist.k
        history = dist.Bt*history + index(target[:,t], dist.bt) - 1;
    end
    @views for t in dist.k:N
        x, y = index(source[:,t], dist.bs), index(target[:,t+1], dist.bt)
        future = y - 1
        src = dist.Bs*history + x - 1
        predicate = dist.Bt*history + future
        state = dist.Bs*predicate + x - 1

        dist.states[state + 1] += 1
        dist.histories[history + 1] += 1
        dist.sources[src + 1] += 1
        dist.predicates[predicate + 1] += 1

        history = predicate - dist.q*(index(target[:,t-dist.k+1], dist.bt) - 1)
    end
    dist
end

function observe!(dist::TransferEntropy, source::AbstractArray{Int,1}, target::AbstractArray{Int,1})
    if length(source) != length(target)
        throw(ArgumentError("time series should have the same number of timesteps"))
    elseif length(target) <= dist.k
        throw(ArgumentError("target series is too short given k=$(dist.k)"))
    elseif any(b -> b < 1, source) || any(b -> b < 1, target)
        throw(ArgumentError("observations must be positive, nonzero"))
    end

    N = length(target) - 1
    dist.N += N - dist.k + 1
    history = 0
    @views for t in 1:dist.k
        history = dist.Bt*history + target[t] - 1;
    end
    @views for t in dist.k:N
        future = target[t + 1] - 1
        src = dist.Bs*history + source[t] - 1
        predicate = dist.Bt*history + future
        state = dist.Bs*predicate + source[t] - 1

        dist.states[state + 1] += 1
        dist.histories[history + 1] += 1
        dist.sources[src + 1] += 1
        dist.predicates[predicate + 1] += 1

        history = predicate - dist.q*(target[t - dist.k + 1] - 1)
    end
    dist
end

function estimate(dist::TransferEntropy)
    entropy(dist.sources, dist.N) +
    entropy(dist.predicates, dist.N) -
    entropy(dist.states, dist.N) -
    entropy(dist.histories, dist.N)
end

function transferentropy!(dist::TransferEntropy, source::AbstractArray{Int},
                          target::AbstractArray{Int})
    estimate(observe!(dist, source, target))
end

function transferentropy(source::AbstractArray{Int}, target::AbstractArray{Int}; kwargs...)
    estimate(TransferEntropy(source, target; kwargs...))
end

function transferentropy(::Type{Kraskov1}, xs::AbstractMatrix{Float64}, ys::AbstractMatrix{Float64};
                         k::Int=1, τ::Int=1, delay::Int=1, nn::Int=1, metric::Metric=Chebyshev())
    hs = history(ys, k, τ, delay)

    start = size(ys, 2) - size(hs, 2) + 1
    fs = @view ys[:, start:end]
    ss = @view xs[:, start-delay:end-delay]

    mutualinfo(Kraskov1, fs, ss, hs; nn=nn, metric=metric)
end

function transferentropy(::Type{Kraskov1}, xs::AbstractMatrix{Float64}, ys::AbstractMatrix{Float64},
                         cond::AbstractMatrix{Float64}, conds::AbstractMatrix{Float64}...;
                         k::Int=1, τ::Int=1, delay::Int=1, nn::Int=1, metric::Metric=Chebyshev())
    hs = history(ys, k, τ, delay)

    start = size(ys, 2) - size(hs, 2) + 1
    fs = @view ys[:, start:end]
    ss = @view xs[:, start-delay:end-delay]
    condview = @view cond[:, start-delay:end-delay]
    condviews = [@view c[:, start-delay:end-delay] for c in conds]

    mutualinfo(Kraskov1, fs, ss, hs, condview, condviews...; nn=nn, metric=metric)
end

function transferentropy(::Type{Kraskov}, xs::AbstractMatrix{Float64}, ys::AbstractMatrix{Float64},
                         cond::AbstractMatrix{Float64}...; kwargs...)
    transferentropy(Kraskov1, xs, ys, cond...; kwargs...)
end

function transferentropy(xs::AbstractMatrix{Float64}, ys::AbstractMatrix{Float64},
                         cond::AbstractMatrix{Float64}...; kwargs...)
    transferentropy(Kraskov1, xs, ys, cond...; kwargs...)
end

function transferentropy(xs::AbstractVector{Float64}, ys::AbstractVector{Float64},
                         cond::AbstractVector{Float64}...; kwargs...)
    conds = map(c -> reshape(c, 1, length(c)), cond)
    transferentropy(reshape(xs, 1, length(xs)), reshape(ys, 1, length(ys)), conds...; kwargs...)
end


text:
import REPL

repl_module = Module()
is_running = true
while is_running
    input_fifo_path = ARGS[1]
    output_fifo_path = ARGS[2]

    open(input_fifo_path, "r") do input_fifo
        command = read(input_fifo, String)
        # input command format:
        # eval code
        #  |
        #  --> evals line of code. \n is interpreted as newline
        # eval_scratch code
        # complete code
        # reset
        #  |
        #  --> resets the module state
        # quit

        cmd_end = findfirst(' ', command)
        if cmd_end == nothing
            cmd_end = length(command)+1
        end
        cmd_name = command[1:cmd_end-1]
        if cmd_name == "eval"
           try
               result = Base.eval(repl_module, Meta.parse(command[cmd_end+1:end]))
               open(output_fifo_path, "w") do fifo
                   print(fifo, result)
               end
           catch e
           end
        elseif cmd_name == "complete"
           to_complete = rstrip(command[cmd_end+1:end])
           expr = "REPL.REPLCompletions.completions(\"$(escape_string(to_complete))\", $(length(to_complete)))"
           println(expr)
           completions_res = REPL.REPLCompletions.completions(String(to_complete), length(to_complete), repl_module)
           completions = completions_res[1]
           completion_range = completions_res[2]
           kakoune_completions = String[]
           completion_success = completions_res[3] && !isempty(completion_range) && !isempty(completions)
           if completion_success
               for c in completions
                   completion_cmd = ""
                   completion_unhandled = false
                   if c isa REPL.REPLCompletions.BslashCompletion
	                   completion_text = REPL.REPLCompletions.completion_text(c)
                       completion_cmd = "echo \'a\'"
	                   completion_menu = completion_text
                   elseif c isa REPL.REPLCompletions.ModuleCompletion
	                   completion_text = REPL.REPLCompletions.completion_text(c)
                       completion_cmd = ""
                       completion_menu = "$(c.mod)   {MenuInfo}$(c.mod)"
                   else
                       completion_unhandled = true
                   end

                   if !completion_unhandled
	                   push!(kakoune_completions, "\"$completion_text|$completion_cmd|$completion_menu\"")
                   end
               end
           end
           completion_success &= !isempty(kakoune_completions)
           open(output_fifo_path, "w") do fifo
               if completion_success
                   output = "$(completion_range.start):$(length(completion_range))"
                   output *= " "
                   output *= join(kakoune_completions, " ")
                   print(fifo, output)
                   println(output)
               end
           end
        elseif cmd_name == "reset"
           global repl_module = create_repl_module()
        elseif cmd_name == "quit"
           global is_running = false
        end
       flush(stdout)
    end
end

max_stars_repo_name: doppioandante/kakoune-julia-mode

max_issues_repo_name: doppioandante/kakoune-julia-mode

max_forks_repo_name: doppioandante/kakoune-julia-mode


text:
# Internal node (INode). For a large dictionary, the MuttsDict consists of a
# tree of depth 7, where the first 6 nodes along any path are INodes.
# INodes are trie nodes with a branching factor 2^k for some k.  The branching
# factor k increases with N, according to the NodeConfig — see _config_for_size(N).
#
# The 64-bit hash value of a key is partitioned among the nodes as follows:
# 0000000000000000000000000000000000000000000000000000000000000000
# +-------+-------+-------+-------+-------+-------+--------------+
# inode1  inode2  inode3  inode4  inode5  inode6  leaf
#
# For example, to navigate from inode1 to the appropriate second-level inode
# for a key, we would choose child (hashcode >> 56) & mask,
# where mask=length(INode.entries)-1.
mutable struct INode{K,V}
    entries::Union{Vector{INode{K,V}},Vector{Leaf{K,V}}}
    is_mutable::Bool

    function INode{K,V}(capacity::Int, child) where {K,V}
        @dassert1 ispow2(capacity)
        entries = Vector{typeof(child)}(undef, capacity)
        for i=1:capacity
            @inbounds entries[i] = child
        end
        new(entries, true)
    end

    # Enlarge the capacity of an INode.  For example, when going from
    # children[1..4] to capacity 8, the resulting INode would have entries
    # [children[1], ..., children[4], children[1], ..., children[4]]
    function INode{K,V}(capacity::Int, children::Vector{T}) where {K,V,T}
        @dassert1 ispow2(capacity)
        @dassert1 ispow2(length(children))
        @dassert1 capacity > length(children)
        entries = typeof(children)(undef, capacity)
        pos = 1
        for i=1:capacity
            @dassert1 !is_mutable(children[pos])
            entries[i] = children[pos]
            pos += 1
            if pos == length(children)+1
                pos = 1
            end
        end
        new(entries,true)
    end

    function INode{K,V}(entries::Vector{T}, is_mutable::Bool) where {K,V,T}
        new(entries, is_mutable)
    end
end

is_mutable(inode::INode{K,V}) where {K,V} = inode.is_mutable

function _get_mutable_version(inode::INode{K,V}, config::NodeConfig, height::Int) where {K,V}
    depth = 1 + length(config.inode_capacity) - height
    should_have_capacity = config.inode_capacity[depth]   # TODO: @inbounds
    if length(inode.entries) != should_have_capacity
        # Need to enlarge the branching factor of this inode. Mark all descendents as
        # immutable, since they will now be reachable by multiple paths
        mark_immutable!(inode)

        # Construct a new inode with increased capacity.  When an inode with
        # entries=[x0,x1,x2,x3] gets increased in size by one bit, the new
        # entries is [x0,x1,x2,x3,x0,x1,x2,x3].
        return INode{K,V}(should_have_capacity, inode.entries)
    elseif inode.is_mutable
        return inode
    else
        new_inode = INode{K,V}(copy(inode.entries), true)
        return new_inode
    end
end

function mark_immutable!(inode::INode{K,V}) where {K,V}
    (!inode.is_mutable) && return
    entries = inode.entries
    if isa(entries, Vector{Leaf{K,V}})
        # Branch specialized for Vector{Leaf{K,V}}
        for child in entries
            mark_immutable!(child)
        end
    else
        # Branch specialized for Vector{INode{K,V}}
        for child in entries
            mark_immutable!(child)
        end
    end
    inode.is_mutable = false
end


text:
import Reproduce

abstract type LoggerKeyData end

function lg_step!(self::LoggerKeyData, env, agent, s, a, s_next, r, is_terminal, cur_step_in_episode, cur_step_total)
    error("lg_step! for this logger is unimplemented")
end

function lg_episode_end!(self::LoggerKeyData, cur_step_in_episode, cur_step_total)
    error("lg_episode_end! for this logger is unimplemented")
end

function save_log(self::LoggerKeyData, save_dict::Dict)
    error("save_log for this logger is unimplemented")
end

abstract type ErrorRecorder <: LoggerKeyData end

include("error.jl")
include("goal_visitation.jl")
include("episode_length.jl")
include("mountain_car_error.jl")
include("cumulants.jl")
include("temp_print.jl")
include("value_map.jl")
include("autostep_stepsize.jl")
include("tabular_tmaze_old_error.jl")
include("tabular_tmaze_error_map.jl")
include("intrinsic_reward.jl")
include("behaviour_action_values.jl")
include("one_d_state_visitation.jl")
include("ttmaze_state_visitation.jl")
include("emphasis.jl")
include("ISRatio.jl")

# Module for scoping key names
module LoggerKey
    const GOAL_VISITATION = "GOAL_VISITATION"
    const ONED_GOAL_VISITATION = "ONED_GOAL_VISITATION"
    const EPISODE_LENGTH = "EPISODE_LENGTH"
    const MC_UNIFORM_ERROR = "MC_UNIFORM_ERROR"
    const MC_START_ERROR = "MC_START_ERROR"
    const TEMP_PRINT = "TEMP_PRINT"
    const TTMAZE_ERROR = "TTMAZE_ERROR"
    const VALUE_MAP = "VALUE_MAP"
    const TWODGRIDWORLDERROR = "TWODGRIDWORLDERROR"
    const TWODGRIDWORLDERRORDPI = "TWODGRIDWORLDERRORDPI"
    const AUTOSTEP_STEPSIZE = "AUTOSTEP_STEPSIZE"
    const TTMAZE_UNIFORM_ERROR = "TTMAZE_UNIFORM_ERROR"
    const TTMAZE_OLD_ERROR = "TTMAZE_OLD_ERROR"
    const TTMAZE_ERROR_MAP = "TTMAZE_ERROR_MAP"
    const INTRINSIC_REWARD = "INTRINSIC_REWARD"
    const TTMAZE_DIRECT_ERROR = "TTMAZE_DIRECT_ERROR"
    const TTMAZE_ROUNDROBIN_ERROR = "TTMAZE_ROUNDROBIN_ERROR"
    const BEHAVIOUR_ACTION_VALUES = "BEHAVIOUR_ACTION_VALUES"
    const ONED_STATE_VISITATION = "ONED_STATE_VISITATION"
    const CUMULANTS = "CUMULANTS"
    const WC_PER_DEMON = "WC_PER_DEMON"
    const ONEDTMAZEERROR = "ONEDTMAZEERROR" # this old error metric gives true uniform distribution, but this is not reachable with no action noise
    const ONEDTMAZEERROR_DPI = "ONEDTMAZEERROR_DPI" # error based on following d_pi for each GVF in the onedtmaze
    const ONEDTMAZEERROR_UNIFORM = "ONEDTMAZEERROR_UNIFORM" # error with states following RR but using uniform aciton weightings
    const ONEDTMAZEERROR_DMU = "ONEDTMAZEERROR_DMU" # State weighting according to round robin with exploring starts
    const TTMAZE_STATE_VISITATION = "TTMAZE_STATE_VISITATION"
    const EMPHASIS = "EMPHASIS"
    const IS_RATIO = "IS_RATIO"
end

module LoggerInitKey
    const TOTAL_STEPS = "TOTAL_STEPS"
    const INTERVAL = "INTERVAL"
    const ENV = "ENV"
end

const LOGGER_KEY_MAP = Dict(
    LoggerKey.GOAL_VISITATION => GoalVisitation,
    LoggerKey.EPISODE_LENGTH => EpisodeLength,
    LoggerKey.MC_UNIFORM_ERROR => MCErrorUniform,
    LoggerKey.MC_START_ERROR => MCErrorStartStates,
    LoggerKey.TTMAZE_ERROR => TTMazeError,
    LoggerKey.TEMP_PRINT => TempPrint,
    LoggerKey.VALUE_MAP => ValueMap,
    LoggerKey.ONEDTMAZEERROR => OneDTMazeError,
    LoggerKey.TWODGRIDWORLDERROR => TwoDGridWorldError,
    LoggerKey.TWODGRIDWORLDERRORDPI => TwoDGridWorldErrorCenterDPI,
    LoggerKey.ONED_GOAL_VISITATION => OneDGoalVisitation,
    LoggerKey.AUTOSTEP_STEPSIZE => AutostepStepSize,
    LoggerKey.TTMAZE_UNIFORM_ERROR => TTMazeUniformError,
    LoggerKey.TTMAZE_OLD_ERROR => TTMazeOldError,
    LoggerKey.TTMAZE_ERROR_MAP => TTMazeErrorMap,
    LoggerKey.TTMAZE_ROUNDROBIN_ERROR => TTMazeRoundRobinError,
    LoggerKey.INTRINSIC_REWARD => IntrinsicRewardLogger,
    LoggerKey.TTMAZE_DIRECT_ERROR => TTMazeDirectError,
    LoggerKey.BEHAVIOUR_ACTION_VALUES => BehaviourActionValues,
    LoggerKey.ONED_STATE_VISITATION => OneDStateVisitation,
    LoggerKey.CUMULANTS => Cumulants,
    LoggerKey.WC_PER_DEMON => WC_Demon_Logger,
    LoggerKey.ONEDTMAZEERROR_DPI => OneDTMazeError_DPI,
    LoggerKey.ONEDTMAZEERROR_UNIFORM => OneDTMazeError_Uniform,
    LoggerKey.TTMAZE_STATE_VISITATION => TTMazeStateVisitation,
    LoggerKey.EMPHASIS => Emphasis,
    LoggerKey.IS_RATIO => ISRatio,
    LoggerKey.ONEDTMAZEERROR_DMU => OneDTMazeError_dmu,
)

# Common logger for all experiments. It has multiple functionalities so pass in what you need to get started
mutable struct Logger
    save_file:: String
    # Tracks what needs to be logged
    logger_keys::Array{String}
    # Array of structs for logger logic
    logger_key_data::Array{LoggerKeyData}

    # This should only track info that makes sense to track across ALL loggers.
    cur_step_in_episode::Int64
    cur_step_total::Int64

    function Logger(parsed, logger_init_info, save_file)
        logger_keys = parsed["logger_keys"]
        logger_key_data = LoggerKeyData[]
        save_results = Dict()

        cur_step_in_episode = 0;
        cur_step_total = 0;

        # Set up for each key
        for key in logger_keys
            if haskey(LOGGER_KEY_MAP, key)
                push!(logger_key_data, LOGGER_KEY_MAP[key](logger_init_info))
            else
                throw(ArgumentError("Logger Key $(key) is not in the LOGGER_KEY_MAP"))
            end
        end
        new(save_file, logger_keys, logger_key_data, cur_step_in_episode, cur_step_total)
    end
end

function logger_start!(self::Logger, env, agent)
    for data in self.logger_key_data
        if applicable(lg_start!, data, env, agent)
            lg_start!(data, env, agent)
        # else
        #     @warn "No lg_start! defined for $(typeof(data))"
        end
    end
end

function logger_step!(self::Logger, env, agent, s, a, s_next, r, is_terminal)
    self.cur_step_in_episode += 1;
    self.cur_step_total += 1;

    for data in self.logger_key_data
        lg_step!(data, env, agent, s, a, s_next, r, is_terminal, self.cur_step_in_episode, self.cur_step_total)
    end
end

function logger_episode_end!(self::Logger)
    for data in self.logger_key_data
        lg_episode_end!(data, self.cur_step_in_episode, self.cur_step_total)
    end

    self.cur_step_in_episode = 0
end

function logger_save(self::Logger)
    save_dict = Dict()
    for data in self.logger_key_data
        save_log(data, save_dict)
    end
    save_results(self.save_file, save_dict)
end

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
# Question
# What is the largest prime factor of the number 600851475143 ?

# Time
# O(n) any better solution??
function sieve(a)
    # find all the prime numbers less than or equal to a
    sieve = collect(1:a)
    # now from this iterate from 2 and remove all their multiples
    index = ones(Bool, a)
    while true

    end

end

function main()
    val = 600851475143
    sq = convert(Int, floor(sqrt(val)))
    # now search for all the factors which
    # we could use sieve of erosthenes to findout the prime numbers
    sieve(sq)
    println(sq)
end

@time main()


text:
function AbstractPlotting.plot!(
        lscene::LScene, P::AbstractPlotting.PlotFunc,
        attributes::AbstractPlotting.Attributes, args...;
        kw_attributes...)

    plot = AbstractPlotting.plot!(lscene.scene, P, attributes, args...; kw_attributes...)

    plot
end

function AbstractPlotting.plot!(P::AbstractPlotting.PlotFunc, ls::LScene, args...; kw_attributes...)
    attributes = AbstractPlotting.Attributes(kw_attributes)
    AbstractPlotting.plot!(ls, P, attributes, args...)
end

function LScene(fig_or_scene; bbox = nothing, scenekw = NamedTuple(), kwargs...)

    topscene = get_topscene(fig_or_scene)

    default_attrs = default_attributes(LScene, topscene).attributes
    theme_attrs = subtheme(topscene, :LScene)
    attrs = merge!(merge!(Attributes(kwargs), theme_attrs), default_attrs)

    layoutobservables = LayoutObservables{LScene}(attrs.width, attrs.height, attrs.tellwidth, attrs.tellheight,
        attrs.halign, attrs.valign, attrs.alignmode; suggestedbbox = bbox)

    scene = Scene(topscene, lift(round_to_IRect2D, layoutobservables.computedbbox); scenekw...)

    ls = LScene(fig_or_scene, layoutobservables, attrs, Dict{Symbol, Any}(), scene)

    # register as current axis
    # TODO: is this a good place for that? probably not
    if fig_or_scene isa Figure
        AbstractPlotting.current_axis!(fig_or_scene, ls)
    end

    ls
end


text:
@testset "atom" begin
    sdfa = SDFileAtom(:Fe, 2, 1, nothing, [1.0, 2.0, 0.0])
    @test sdfa.symbol == :Fe
    @test sdfa.charge == 2
    @test sdfa.coords == [1.0, 2.0, 0.0]
    @test atomnumber(sdfa) == 26
    @test atomweight(sdfa) == 55.845

    smia = SmilesAtom(:C, 0, 1, 13.0, false, 2)
    @test !smia.isaromatic
    @test smia.stereo == 2
    @test atomname(smia) == "Carbon"
    @test atomweight(smia) == 13.0

    @test atomnumber(:Ra) == 88
    @test atomname(:Ti) == "Titanium"
end


text:
function truncate_increments!(ΔW, A)
    if A > 0
        for i in eachindex(ΔW)
            if ΔW[i] < -A
                ΔW[i] = -A
            elseif ΔW[i] > A
                ΔW[i] = A
            end
        end
    end
end

function update_params!(int::IntegratorSIRK, sol::AtomicSolutionSDE)
    # set time for nonlinear solver and copy previous solution
    int.params.t  = sol.t
    int.params.q .= sol.q
    int.params.ΔW .= sol.ΔW
    int.params.ΔZ .= sol.ΔZ

    # truncate the increments ΔW with A
    truncate_increments!(int.params.ΔW, int.params.A)
end

function update_params!(int::Union{IntegratorSIPRK,IntegratorSISPRK}, sol::AtomicSolutionPSDE)
    # set time for nonlinear solver and copy previous solution
    int.params.t  = sol.t
    int.params.q .= sol.q
    int.params.p .= sol.p
    int.params.ΔW .= sol.ΔW
    int.params.ΔZ .= sol.ΔZ

    # truncate the increments ΔW with A
    truncate_increments!(int.params.ΔW, int.params.A)
end

function update_params!(int::IntegratorWIRK, sol::AtomicSolutionSDE)
    # set time for nonlinear solver and copy previous solution
    int.params.t  = sol.t
    int.params.q .= sol.q
    int.params.ΔW .= sol.ΔW
    int.params.ΔZ .= sol.ΔZ
end

@doc raw"""
Update solution for stochastic Runge-Kutta methods (SIRK and WIRK)

- `x`: the solution vector to be updated
- `V`: the matrix containing the drift vector evaluated at the internal stages v(Q_i) (SIRK) or v(Q0_i) (WIRK)
- `B`: the array containing the diffusion matrix evaluated at the internal stages B(Q_i) (SIRK) or B(Q1^(l)_i) (WIRK)
- `bdrift`: the Runge-Kutta coefficients for the drift part
- `bdiff`:  the Runge-Kutta coefficients for the diffusion part
- `Δt`: the time step
- `ΔW`: the increments of the Brownian motion (SFIRK) or the increments represented by the random variables Î^(k) (WFIRK)
"""
function update_solution!(sol::AtomicSolutionSDE{T}, V::Vector{Vector{T}}, B::Vector{Matrix{T}},
                            bdrift::Vector{T}, bdiff::Vector{T}, Δt::T, ΔW::Vector{T}, Δy::Vector{T}=zero(ΔW)) where {T}

    @assert length(bdrift) == length(bdiff) == length(V) == length(B)

    for i in eachindex(V, B)
        @assert length(sol.q) == length(V[i]) == size(B[i], 1)
        @assert length(ΔW)== size(B[i], 2)
    end

    local Δx::T

    # Contribution from the drift part
    for k in eachindex(sol.q)
        Δx = 0
        for i in eachindex(bdrift, V)
            Δx += bdrift[i] * V[i][k]
        end
        update!(sol, Δt*Δx, k)
    end

    # Contribution from the diffusion part
    for k in eachindex(sol.q)
        Δy .= 0
        for i in eachindex(bdiff, B)
            for l in eachindex(Δy)
                Δy[l] += bdiff[i] * B[i][k,l]
            end
        end
        update!(sol, dot(Δy,ΔW), k)
    end
end

@doc raw"""
Update solution for stochastic Runge-Kutta methods (SERK)

- `x`: the solution vector to be updated
- `V`: the matrix containing the drift vector evaluated at the internal stages v(Q_i)
- `B`: the array containing the diffusion matrix evaluated at the internal stages B(Q_i)
- `bdrift`: the Runge-Kutta coefficients for the drift part
- `bdiff`:  the Runge-Kutta coefficients for the ΔW terms of the diffusion part
- `bdiff2`: the Runge-Kutta coefficients for the ΔZ terms of the diffusion part
- `Δt`: the time step
- `ΔW`: the increments of the Brownian motion
- `ΔZ`: the integrals of the increments of the Brownian motion
"""
function update_solution!(sol::AtomicSolutionSDE{T}, V::Vector{Vector{T}}, B::Vector{Matrix{T}},
                            bdrift::Vector{T}, bdiff::Vector{T}, bdiff2::Vector{T},
                            Δt::T, ΔW::Vector{T}, ΔZ::Vector{T}, Δy::Vector{T}=zero(ΔW)) where {T}

    @assert length(bdiff2) == length(B)

    for i in eachindex(B)
        @assert length(ΔZ)== size(B[i], 2)
    end

    # Contributions from the drift and diffusion parts
    update_solution!(sol, V, B, bdrift, bdiff, Δt, ΔW, Δy)

    # Contribution from the diffusion part (ΔZ terms)
    for k in eachindex(sol.q)
        Δy .= 0
        for i in eachindex(bdiff2, B)
            for l in eachindex(Δy)
                Δy[l] += bdiff2[i] * B[i][k,l]
            end
        end
        update!(sol, dot(Δy,ΔZ)/Δt, k)
    end
end

@doc raw"""
Update solution for stochastic partitioned Runge-Kutta methods

- `q`, `p`: the solution vector to be updated
- `V`, `F`: the matrix containing the drift vectors evaluated at the internal stages v(Q_i), f(Q_i)
- `B`, `G`: the array containing the diffusion matrices evaluated at the internal stages B(Q_i), G(Q_i)
- `bqdrift`, `bpdrift`: the Runge-Kutta coefficients for the drift parts of the q and p equations
- `bqdiff`, `bpdiff`:   the Runge-Kutta coefficients for the diffusion parts of the q and p equations
- `Δt`: the time step
- `ΔW`: the increments of the Brownian motion
"""
function update_solution!(sol::AtomicSolutionPSDE{T},
                          V::Vector{Vector{T}}, F::Vector{Vector{T}},
                          B::Vector{Matrix{T}}, G::Vector{Matrix{T}},
                          bqdrift::Vector{T}, bqdiff::Vector{T},
                          bpdrift::Vector{T}, bpdiff::Vector{T},
                          Δt::T, ΔW::Vector{T}, Δy::Vector{T}=zero(ΔW), Δz::Vector{T}=zero(ΔW)) where {T}

    @assert length(bqdrift) == length(bqdiff) == length(bpdrift) == length(bpdiff) == length(V) == length(F) == length(B) == length(G)

    for i in eachindex(V, F, B, G)
        @assert length(sol.q) == length(sol.p) == length(V[i]) == length(F[i]) == size(B[i], 1) == size(G[i], 1)
        @assert length(ΔW) == size(B[i], 2) == size(G[i], 2)
    end

    local Δq::T
    local Δp::T

    # Contribution from the drift part
    for k in eachindex(sol.q, sol.p)
        Δq = 0
        Δp = 0
        for i in eachindex(bqdrift, bpdrift, V, F)
            Δq += bqdrift[i] * V[i][k]
            Δp += bpdrift[i] * F[i][k]
        end
        update!(sol, Δt*Δq, Δt*Δp, k)
    end

    # Contribution from the diffusion part
    for k in eachindex(sol.q, sol.p)
        Δy .= 0
        Δz .= 0
        for i in eachindex(bqdiff, bpdiff, B, G)
            for l in eachindex(Δy, Δz)
                Δy[l] += bqdiff[i] * B[i][k,l]
                Δz[l] += bpdiff[i] * G[i][k,l]
            end
        end
        update!(sol, dot(Δy,ΔW), dot(Δz,ΔW), k)
    end
end

@doc raw"""
Update solution for stochastic split partitioned Runge-Kutta methods

- `q`, `p`: the solution vector to be updated
- `V`, `F1`, `F2`: the matrix containing the drift vectors evaluated at the internal stages v(Q_i), fi(Q_i)
- `B`, `G1`, `G2`: the array containing the diffusion matrices evaluated at the internal stages B(Q_i), Gi(Q_i)
- `bqdrift`, `bpdrift1`, `bpdrift2`: the Runge-Kutta coefficients for the drift parts of the q and p equations
- `bqdiff`, `bpdiff1`, `bpdiff2`:    the Runge-Kutta coefficients for the diffusion parts of the q and p equations
- `Δt`: the time step
- `ΔW`: the increments of the Brownian motion
"""
function update_solution!(sol::AtomicSolutionPSDE{T},
                            V::Vector{Vector{T}}, F1::Vector{Vector{T}}, F2::Vector{Vector{T}},
                            B::Vector{Matrix{T}}, G1::Vector{Matrix{T}}, G2::Vector{Matrix{T}},
                            bqdrift::Vector{T}, bqdiff::Vector{T},
                            bpdrift1::Vector{T}, bpdrift2::Vector{T},
                            bpdiff1::Vector{T}, bpdiff2::Vector{T},
                            Δt::T, ΔW::Vector{T}, Δy::Vector{T}=zero(ΔW), Δz::Vector{T}=zero(ΔW)) where {T}

   @assert length(bqdrift) == length(bqdiff) == length(bpdrift1) == length(bpdrift2) == length(bpdiff1) == length(bpdiff2) == length(V) == length(F1) == length(F2) == length(B) == length(G1) == length(G2)

   for i in eachindex(V, F1, F2, B, G1, G2)
       @assert length(sol.q) == length(sol.p) == length(V[i]) == length(F1[i]) == length(F2[i]) == size(B[i], 1) == size(G1[i], 1) == size(G2[i], 1)
       @assert length(ΔW) == size(B[i], 2) == size(G1[i], 2) == size(G2[i], 2)
   end

   local Δq::T
   local Δp::T

   # Contribution from the drift part
   for k in eachindex(sol.q, sol.p)
       Δq = 0
       Δp = 0
       for i in eachindex(bqdrift, bpdrift1, bpdrift2, V, F1, F2)
           Δq += bqdrift[i] * V[i][k]
           Δp += bpdrift1[i] * F1[i][k] + bpdrift2[i] * F2[i][k]
       end
       update!(sol, Δt*Δq, Δt*Δp, k)
   end

   # Contribution from the diffusion part
   for k in eachindex(sol.q, sol.p)
       Δy .= 0
       Δz .= 0
       for i in eachindex(bqdiff, bpdiff1, bpdiff2, B, G1, G2)
           for l in eachindex(Δy,Δz)
               Δy[l] += bqdiff[i]  * B[i][k,l]
               Δz[l] += bpdiff1[i] * G1[i][k,l] + bpdiff2[i] * G2[i][k,l]
           end
       end
       update!(sol, dot(Δy,ΔW), dot(Δz,ΔW), k)
   end
end

@doc raw"""
Update solution for weak Runge-Kutta methods WERK

- `x`:  the solution vector to be updated
- `V`:  the matrix containing the drift vector evaluated at the internal stages v(Q_i)
- `B1`: the array containing the diffusion matrix evaluated at the internal stages H^(l)_i, such that B1[:,l,i] is evaluated at H^(l)_i
- `B2`: the array containing the diffusion matrix evaluated at the internal stages Ĥ^(l)_i, such that B2[:,l,i] is evaluated at Ĥ^(l)_i
- `α`:  the Runge-Kutta coefficients for the drift part
- `β1`: the Runge-Kutta coefficients for the diffusion term with the random increments
- `β2`: the Runge-Kutta coefficients for the second diffusion term
- `Δt`: the time step
- `ΔW`: the increments of the Brownian motion represented by the random variables Î^(k)
"""
function update_solution!(sol::AtomicSolutionSDE{T}, V::Vector{Vector{T}},
                            B1::Vector{Matrix{T}}, B2::Vector{Matrix{T}}, α::Vector{T},
                            β1::Vector{T}, β2::Vector{T}, Δt::T, ΔW::Vector{T}, Δy::Vector{T}=zero(ΔW)) where {T}
    @assert length(α) == length(β1) == length(β2) == length(V) == length(B1) == length(B2)

    for i in eachindex(V, B1, B2)
        @assert length(sol.q) == length(V[i]) == size(B1[i], 1) == size(B2[i], 1)
        @assert length(ΔW)== size(B1[i], 2) == size(B2[i], 2)
    end

    local Δx::T

    # Contribution from the drift part
    for k in eachindex(sol.q)
        Δx = 0
        for i in eachindex(α, V)
            Δx += α[i] * V[i][k]
        end
        update!(sol, Δt*Δx, k)
    end

    # Contribution from the diffusion term with the random variables I^(k)_i
    for k in eachindex(sol.q)
        Δy .= 0
        for i in eachindex(β1, B1)
            for l in eachindex(Δy)
                Δy[l] += β1[i] * B1[i][k,l]
            end
        end
        update!(sol, dot(Δy,ΔW), k)
    end

    # Contribution from the second diffusion term
    for k in eachindex(sol.q)
        Δx = 0
        for i in eachindex(β2, B2)
            for l in axes(B2[i], 2)
                Δx += β2[i] * B2[i][k,l]
            end
        end
        update!(sol, sqrt(Δt)*Δx, k)
    end
end


text:
struct AugmentedLagrangian
    λ               # duals - general constraints
    λl              # duals - primal lower bounds
    λu              # duals - primal upper bounds

    ρ               # penalties - general constraints
    ρl              # penalties - primal lower bounds
    ρu              # penalties - primal upper bounds
    ρ0              # initial penalty
    s               # penalty scaling

    as              # active set for general constraints
    asl             # active set for primal lower bounds
    asu             # active set for primal upper bounds

    c               # memory for constraints
    ∇c              # memory for constraints Jacobian

    xl              # primal lower bounds
    xu              # primal upper bounds
    cl              # memory for xl - x
    cu              # memory for x - xu

    idx_ineq        # inequality constraint indices
    idx_l           # primal lower bound indices
    idx_u           # primal upper bound indices
    idx_soc         # second-order cone constraint indices
end

function augmented_lagrangian(n, m;
        ρ0 = 1.0, s = 10.0,
        idx_ineq = (1:0),
        xl = -Inf * ones(n), xu = Inf * ones(n),
        idx_soc = (1:0))

    bool_l = isfinite.(xl)
    bool_u = isfinite.(xu)
    idx_l = (1:length(xl))[bool_l]
    idx_u = (1:length(xu))[bool_u]
    ml = sum(bool_l)
    mu = sum(bool_u)

    λ = zeros(m)
    λl = zeros(ml)
    λu = zeros(mu)

    ρ = ρ0 * ones(m)
    ρl = ρ0 * ones(ml)
    ρu = ρ0 * ones(mu)

    as = ones(Bool, m)
    asl = ones(Bool, ml)
    asu = ones(Bool, mu)

    c = zeros(m)
    ∇c = spzeros(m, n)

    cl = zeros(ml)
    cu = zeros(mu)

    AugmentedLagrangian(λ, λl, λu,
        ρ, ρl, ρu, ρ0, s,
        as, asl, asu,
        c, ∇c,
        xl, xu, cl, cu,
        idx_ineq, idx_l, idx_u, idx_soc)
end

function update!(al::AugmentedLagrangian)
    dual_update!(al)
    penalty_update!(al)
    active_set_update!(al)
end

function dual_update!(al::AugmentedLagrangian)
    al.λ .+= al.ρ .* al.c
    al.λl .+= al.ρl .* al.cl
    al.λu .+= al.ρu .* al.cu

    inequality_projection!(al)
    #TODO: soc project
end

function inequality_projection!(al::AugmentedLagrangian)
    al.λ[al.idx_ineq] = max.(0.0, view(al.λ, al.idx_ineq))
    al.λl .= max.(0.0, al.λl)
    al.λu .= max.(0.0, al.λu)
end

function soc_projection!(al::AugmentedLagrangian)
    return nothing
end

function penalty_update!(al::AugmentedLagrangian)
    al.ρ .*= al.s
    al.ρl .*= al.s
    al.ρu .*= al.s
end

function active_set_update!(al::AugmentedLagrangian)
    fill!(al.as, 1)
    fill!(al.asl, 1)
    fill!(al.asu, 1)

    al.as[al.idx_ineq] = ((view(al.λ, al.idx_ineq) .> 0.0) .| (view(al.c, al.idx_ineq) .>= 0.0)) # (view(al.c, al.idx_ineq) .> 0.0) #
    al.asl .= ((al.λl .> 0.0) .| (al.cl .>= 0.0)) # (al.cl .> 0.0) #
    al.asu .= ((al.λu .> 0.0) .| (al.cu .>= 0.0)) # (al.cu .> 0.0) #
    #TODO: soc as
end

function bounds!(al::AugmentedLagrangian, x)
    al.cl .= view(al.xl, al.idx_l) - view(x, al.idx_l)
    al.cu .= view(x, al.idx_u) - view(al.xu, al.idx_u)
end

function constraint_violation(al::AugmentedLagrangian, x)
    c!(al.c, x)
    bounds!(al, x)
    al.c[al.idx_ineq] = max.(0.0, al.c[al.idx_ineq])
    return max(norm(al.c, Inf),
               norm(max.(0.0, al.cl), Inf),
               norm(max.(0.0, al.cu), Inf))
end

function reset!(al::AugmentedLagrangian)
    fill!(al.λ, 0.0)
    fill!(al.λl, 0.0)
    fill!(al.λu, 0.0)
    fill!(al.ρ, al.ρ0)
    fill!(al.ρl, al.ρ0)
    fill!(al.ρu, al.ρ0)
    fill!(al.as, 1)
    fill!(al.asl, 1)
    fill!(al.asu, 1)
    fill!(al.c, 0.0)
    fill!(al.cl, 0.0)
    fill!(al.cu, 0.0)
end

# n = 20
# m = 10
# al = augmented_lagrangian(n, m, xl = zeros(n))
# update!(al)
# reset!(al)
#
# al.asl
# al.λl
# al.cl
# ((al.λl .== 0.0) .& (al.cl .< 0.0))


text:
using Documenter, ContinuousTransformations

makedocs(modules = [ContinuousTransformations],
         format = :html,
         clean = true,
         sitename = "ContinuousTransformations.jl",
         authors = "Tamás K. Papp",
         checkdocs = :all,
         linkcheck = true,
         assets = ["assets/custom.css"],
         html_prettyurls = haskey(ENV, "TRAVIS"), # clean URLs building on Travis
         pages = [
             "Overview" => "index.md",
             "General API" => "general.md",
             "Intervals and univariate transformations" => "univariate.md",
             "Grouped transformations" => "grouped.md",
             "Wrapped transformations" => "wrapped.md",
             "Internals" => "internals.md"
         ])

deploydocs(repo = "github.com/tpapp/ContinuousTransformations.jl.git",
           target = "build",
           deps = nothing,
           make = nothing,
           julia = "0.6")


text:
using RadiativeTransfer
using RadiativeTransfer.Architectures
using RadiativeTransfer.Absorption
using RadiativeTransfer.Scattering
using Test
using DelimitedFiles
using Statistics
using ProgressMeter
using WignerSymbols
using Distributions
using JLD2
using CUDA

# Test the Cross Section module
@testset "RadiativeTransfer.Absorption" begin include("test_Absorption.jl") end

# Test the Phase Function module
@testset "RadiativeTransfer.Scattering" begin include("test_Scattering.jl") end


text:
#=
# 2. Basic flow with a stationary body
In this notebook we will simulate the flow past a stationary body.
=#

using ViscousFlow
#-
using Plots

#=
### The basic steps
From the previous notebook, we add one additional step:
* **Specify the problem**: Set the Reynolds number, the free stream, and any other problem parameters
* **Discretize**: Set up a solution domain, choose the grid Reynolds number and the critical time step limits
* **Set up bodies**: *Create the body or bodies and specify their motions, if any*
* **Construct the system structure**: Create the operators that will be used to perform the simulation
* **Initialize**: Set the initial flow field and initialize the integrator
* **Solve**: Solve the flow field
* **Examine**: Examine the results
The rest of the steps are nearly the same as in the previous example.

As before, we initialize the parameters dictionary:
=#
my_params = Dict()

#=
### Problem specification
Set the Reynolds number and free stream. We will set the free stream to be in the $x$ direction, with speed equal to 1.
=#
my_params["Re"] = 200
my_params["freestream speed"] = 1.0
my_params["freestream angle"] = 0.0

#=
### Discretize
We set the grid Re to 4.0 here to get a quicker solution, though it is generally
better to make this smaller (it defaults to 2.0).
=#
xlim = (-1.0,5.0)
ylim = (-2.0,2.0)
my_params["grid Re"] = 4.0
g = setup_grid(xlim,ylim,my_params)

#=
### Set up bodies
Here, we will set up a rectangle of half-height 0.5 and half-width 0.25
at 45 degrees angle of attack. We also need to supply the spacing between points. For this, we use
the function [`surface_point_spacing`](@ref)
=#
Δs = surface_point_spacing(g,my_params)
body = Rectangle(0.5,0.25,Δs)

#=
We place the body at a desired location and orientation with the `RigidTransform`
function. This function creates an operator `T` that acts in-place on the body:
after the operation is applied, `body` is transformed to the correct location/orientation.
=#
cent = (0.0,0.0) # center of body
α = 45π/180 # angle
T = RigidTransform(cent,α)
T(body) # transform the body to the current configuration

# Let's plot it just to make sure
plot(body,xlim=xlim,ylim=ylim)

#=
### Construct the system structure
This step is like the previous notebook, but now we also provide the body as an argument. It is important
to note that we have not provided any explicit information about the boundary conditions on our shape.
It therefore assumes that we want to enforce zero velocity on the shape. We will show another
example later in which we change this.
=#
sys = viscousflow_system(g,body,phys_params=my_params);

#=
### Initialize
Now, we initialize with zero vorticity. Note that we do this by calling
`init_sol` with no argument except for `sys` itself.
=#
u0 = init_sol(sys)

#=
and now create the integrator, with a long enough time span to hold the whole
solution history:
=#
tspan = (0.0,20.0)
integrator = init(u0,tspan,sys)

#=
### Solve
Now we are ready to solve the problem. Let's advance the solution to $t = 1$.
=#
@time step!(integrator,1.0)

#=
### Examine
Let's look at the flow field at the end of this interval
=#
plot(
plot(vorticity(integrator),sys,title="Vorticity",clim=(-15,15),levels=range(-15,15,length=30), color = :RdBu,ylim=ylim),
plot(streamfunction(integrator),sys,title="Streamlines",ylim=ylim,color = :Black),
    size=(700,300)
    )

#=
#### Compute the force history
To do this, we supply the solution history `sol`, the system `sys`, and the index
of the body (1).
=#
sol = integrator.sol;
fx, fy = force(sol,sys,1);

#=
Plot the histories. Note that we are actually plotting the drag and lift
coefficient histories here:
$$ C_D = \dfrac{F_x}{\frac{1}{2}\rho U_\infty^2 L}, \quad C_L = \dfrac{F_y}{\frac{1}{2}\rho U_\infty^2 L} $$
Since the quantities in this simulation are already scaled by $\rho$, $U_\infty$, and $L$
(because $\rho$ has been scaled out of the equations, and the free stream speed is
set to 1 and the height of the shape to 1), then we obtain these coefficients by
simply dividing by 1/2, or equivalently, by multiplying by 2:
=#
plot(
plot(sol.t,2*fx,xlim=(0,Inf),ylim=(0,4),xlabel="Convective time",ylabel="\$C_D\$",legend=:false),
plot(sol.t,2*fy,xlim=(0,Inf),ylim=(-4,4),xlabel="Convective time",ylabel="\$C_L\$",legend=:false),
    size=(800,350)
)

# The mean drag and lift coefficients (omitting the first two steps) are
meanCD = GridUtilities.mean(2*fx[3:end])
#-
meanCL = GridUtilities.mean(2*fy[3:end])


text:
module testUniverseModels

using Base.Test

include(string(Pkg.dir("AssetMgmt"), "/src/AssetMgmt.jl"))

println("\n Running UniverseModel tests\n")

end


text:
using Tissue
import Tissue as T
import Base.Threads

const MAX_GEN = 3

accumulator = Threads.Atomic{Int}(0)

mutable struct SourceCalculator
    last::Int
    SourceCalculator() = new(1)
end

function T.process(c::SourceCalculator)
    if c.last > MAX_GEN
        return nothing
    end

    ret_val = c.last
    c.last += 1

    return ret_val
end

struct AddConstantCalculator
    constant::Int
end

function T.process(c::AddConstantCalculator, num_in::Int; graph::Graph)::Int
    return num_in + c.constant
end

struct SubtractConstantCalculator
    constant::Int
end

function T.process(c::SubtractConstantCalculator, num_in::Int)::Int
    return num_in - c.constant
end

struct MultiplyCalculator end

function T.process(c::MultiplyCalculator, first_num::Int, second_num::Int)::Int
    return first_num * second_num
end

mutable struct PrinterCalculator end

function T.process(c::PrinterCalculator, num_to_print)
    accumulator[] += num_to_print

    nothing
end

@graph NumberGraph begin
    # 1. Declare calculators.
    @calculator source = SourceCalculator()
    @calculator add1 = AddConstantCalculator(1)
    @calculator sub2 = SubtractConstantCalculator(2)
    @calculator mult = MultiplyCalculator()
    @calculator printer = PrinterCalculator()

    # 2. Declare the streams which connect the calculators together
    @bindstreams add1 num_in = source
    @bindstreams sub2 num_in = source
    @bindstreams mult (first_num = add1) (second_num = sub2)
    @bindstreams printer num_to_print=mult
end

graph = NumberGraph()

###############################
# Before the graph is started
###############################
@test typeof(T.get_calculator(T.get_source_calculator_wrapper(graph))) ==
      SourceCalculator

@test length(T.get_calculator_wrappers(graph)) == 5

expected_values = Dict(
    SourceCalculator => Dict(
        :len_input_ch => 0,
        :len_output_ch => 2,
        :is_sink => false,
        :has_graph_kw => false,
    ),
    AddConstantCalculator => Dict(
        :len_input_ch => 1,
        :len_output_ch => 1,
        :is_sink => false,
        :has_graph_kw => true,
    ),
    SubtractConstantCalculator => Dict(
        :len_input_ch => 1,
        :len_output_ch => 1,
        :is_sink => false,
        :has_graph_kw => false,
    ),
    MultiplyCalculator => Dict(
        :len_input_ch => 2,
        :len_output_ch => 1,
        :is_sink => false,
        :has_graph_kw => false,
    ),
    PrinterCalculator => Dict(
        :len_input_ch => 1,
        :len_output_ch => 0,
        :is_sink => true,
        :has_graph_kw => false,
    ),
)

for cw in T.get_calculator_wrappers(graph)
    vals_for_calc = expected_values[typeof(T.get_calculator(cw))]

    @test length(T.get_input_channels(cw)) == vals_for_calc[:len_input_ch]
    @test length(T.get_output_channels(cw)) == vals_for_calc[:len_output_ch]
    @test T.is_sink(cw) == vals_for_calc[:is_sink]
    @test T.process_has_graph_kw(cw) == vals_for_calc[:has_graph_kw]
end

###############################
# After the graph is started
###############################

@test T.is_done(graph) == false
@test T.get_num_sinks_not_init(graph) == 1
@test length(T.get_cw_tasks(graph)) == 0

T.start(graph)

# The SourceCalculator does not run in a cw_task
@test length(T.get_cw_tasks(graph)) == 4

T.wait_until_done(graph)

@test length(T.get_cw_tasks(graph)) == 0
@test T.is_done(graph) == true

@test accumulator[] == 2


text:
# Alternatively, we could use `InteractiveUtils.subtypes` to get a list of it.
threshold_methods = (
    Balanced,
    Entropy,
    Intermodes,
    MinimumError,
    MinimumIntermodes,
    Moments,
    Otsu,
    UnimodalRosin,
    Yen
)

"""
    SingleHistogramThreshold <: AbstractImageBinarizationAlgorithm
    SingleHistogramThreshold(alg::ThresholdAlgorithm; nbins=256)

    binarize([T,] img, f::ThresholdAlgorithm; nbins=256)
    binarize!([out,] img, f::ThresholdAlgorithm; nbins=256)

Binarizes the image `img` using the threshold found by given threshold finding algorithm `alg`.

# Output

Return the binarized image as an `Array{Gray{T}}` of size `size(img)`. If `T` is not specified, it
is inferred from `out` and `img`.

# Arguments

The function argument is described in more detail below.

##  `img::AbstractArray`

The image that needs to be binarized.  The image is automatically converted to `Gray` in order to
construct the requisite graylevel histogram.

## `alg::ThresholdAlgorithm`

`ThresholdAlgorithm` is an Abstract type defined in `ThresholdAlgorithm.jl`, it provides various
threshold finding algorithms:

$(mapreduce(x->"- [`"*string(x)*"`](@ref)", (x,y)->x*"\n"*y, threshold_methods))

For the more detailed explaination and the construction, please refer to each concrete algorithm.
For example, type `?Otsu` in REPL will give you more details on how to use `Otsu` methods.

## `nbins::Integer`

The number of discrete bins that used to build the histogram. A smaller `nbins` could possibly gives
a less noisy, or in other words, a smoother output. The default value is `256`.

# Examples

All the usage follows the same pattern, take `Otsu` as an example:

```julia
using TestImages, ImageBinarization

img = testimage("cameraman")
img_binary = binarize(img, Otsu())
```

It is less convenient, but still, you could also construct a `SingleHistogramThreshold` by yourself:

```julia
using TestImages, ImageBinarization

img = testimage("cameraman")
f = SingleHistogramThreshold(Otsu(), nbins=256)
img_binary = binarize(img, f)
```

"""
struct SingleHistogramThreshold{T} <: AbstractImageBinarizationAlgorithm
    alg::T
    nbins::Int
    function SingleHistogramThreshold(alg::T; nbins::Integer) where T <: ThresholdAlgorithm
        new{T}(alg, nbins)
    end
end

function binarize!(out::GenericGrayImage, img, f::ThresholdAlgorithm, args...; nbins::Integer=256, kwargs...)
    binarize!(out, img, SingleHistogramThreshold(f, nbins=nbins), args...; kwargs...)
end
function binarize!(img::GenericGrayImage, f::ThresholdAlgorithm, args...; nbins::Integer=256, kwargs...)
    binarize!(img, SingleHistogramThreshold(f, nbins=nbins), args...; kwargs...)
end
function binarize(::Type{T}, img, f::ThresholdAlgorithm, args...; nbins::Integer=256, kwargs...) where T
    binarize(T, img, SingleHistogramThreshold(f, nbins=nbins), args...; kwargs...)
end
function binarize(img, f::ThresholdAlgorithm, args...; nbins::Integer=256, kwargs...)
    binarize(ccolor(Gray, eltype(img)), img, SingleHistogramThreshold(f, nbins=nbins), args...; kwargs...)
end
function binarize(img::AbstractArray{T}, f::ThresholdAlgorithm, args...; kwargs...) where T <: Number
    # issue #46: Do not promote Number to Gray{<:Number}
    binarize(T, img, f, args...; kwargs...)
end

function (f::SingleHistogramThreshold)(out::GenericGrayImage, img::GenericGrayImage)
    edges, counts = build_histogram(img, f.nbins)
    t = find_threshold(f.alg, counts[1:end], edges)
    @. out = img > t # here we rely on implicit type conversion to `eltype(out)`
end

function (f::SingleHistogramThreshold)(out::GenericGrayImage, img::AbstractArray{<:Color3})
    # map `img` to grayspace while keeping the storage type
    # TODO: this mapping operation can be done lazily to reduce memory allocation
    GT = base_color_type(eltype(out)){eltype(eltype(img))}
    f(out, GT.(img))
end


text:
## gradient methods

"""
    gradient(inFLD::gcmfaces)

Compute spatial derivatives.
"""
function gradient(inFLD::gcmfaces)
(dFLDdx, dFLDdy)=gradient(inFLD,true)
return dFLDdx, dFLDdy
end

"""
    gradient(inFLD::gcmfaces,doDIV::Bool)

Compute spatial derivatives with or without dividing by grid scale.
"""
function gradient(inFLD::gcmfaces,doDIV::Bool)

exFLD=exchange(inFLD,1)
dFLDdx=similar(inFLD)
dFLDdy=similar(inFLD)

for a=1:inFLD.nFaces;
  (s1,s2)=fsize(exFLD,a)
  tmpA=view(exFLD.f[a],2:s1-1,2:s2-1)
  tmpB=tmpA-view(exFLD.f[a],1:s1-2,2:s2-1)
  tmpC=tmpA-view(exFLD.f[a],2:s1-1,1:s2-2)
  if doDIV
    dFLDdx.f[a]=tmpB./MeshArrays.DXC.f[a]
    dFLDdy.f[a]=tmpC./MeshArrays.DYC.f[a]
  else
    dFLDdx.f[a]=tmpB
    dFLDdy.f[a]=tmpC
  end
end

return dFLDdx, dFLDdy
end

"""
    gradient(inFLD::gcmfaces,iDXC::gcmfaces,iDYC::gcmfaces)

Compute spatial derivatives and multiply by inverse grid scale.
"""
function gradient(inFLD::gcmfaces,iDXC::gcmfaces,iDYC::gcmfaces)

exFLD=exchange(inFLD,1)
dFLDdx=similar(inFLD)
dFLDdy=similar(inFLD)

for a=1:inFLD.nFaces;
  (s1,s2)=fsize(exFLD,a)
  tmpA=view(exFLD.f[a],2:s1-1,2:s2-1)
  tmpB=tmpA-view(exFLD.f[a],1:s1-2,2:s2-1)
  tmpC=tmpA-view(exFLD.f[a],2:s1-1,1:s2-2)
  dFLDdx.f[a]=tmpB.*iDXC.f[a]
  dFLDdy.f[a]=tmpC.*iDYC.f[a]
end

return dFLDdx, dFLDdy
end

## mask methods

"""
    mask(fld::gcmfaces)

Call mask(fld,NaN)
"""
function mask(fld::gcmfaces)
fldmsk=mask(fld,NaN)
return fldmsk
end

"""
    mask(fld::gcmfaces, val::Number)

Replace non finite values with val
"""
function mask(fld::gcmfaces, val::Number)
  fldmsk=similar(fld)
  for a=1:fld.nFaces
    tmp1=copy(fld.f[a])
    replace!(x -> !isfinite(x) ? val : x, tmp1 )
    fldmsk.f[a]=tmp1
  end
  return fldmsk
end

"""
    mask(fld::gcmfaces, val::Number)

Replace noval instances with val
"""
function mask(fld::gcmfaces, val::Number, noval::Number)
  fldmsk=similar(fld)
  for a=1:fld.nFaces
    tmp1=copy(fld.f[a])
    replace!(x -> x==noval ? val : x, tmp1  )
    fldmsk.f[a]=tmp1
  end
  return fldmsk
end

## convergence methods

"""
    convergence(uFLD::gcmfaces,vFLD::gcmfaces)

Compute convergence of a vector field
"""
function convergence(uFLD::gcmfaces,vFLD::gcmfaces)

#important note:
#  Normally uFLD, vFLD should not contain any NaN;
#  if otherwise then something this may be needed:
#  uFLD=mask(uFLD,0.0); vFLD=mask(vFLD,0.0);

CONV=similar(uFLD)

(tmpU,tmpV)=exch_UV(uFLD,vFLD)
for a=1:tmpU.nFaces
  (s1,s2)=fsize(uFLD,a)
  tmpU1=view(tmpU.f[a],1:s1,1:s2)
  tmpU2=view(tmpU.f[a],2:s1+1,1:s2)
  tmpV1=view(tmpV.f[a],1:s1,1:s2)
  tmpV2=view(tmpV.f[a],1:s1,2:s2+1)
  CONV.f[a]=tmpU1-tmpU2+tmpV1-tmpV2
end

return CONV
end

## smooth function

"""
    smooth(FLD::gcmfaces,DXCsm::gcmfaces,DYCsm::gcmfaces)

Smooth out scales below threshold
"""
function smooth(FLD::gcmfaces,DXCsm::gcmfaces,DYCsm::gcmfaces)

#important note:
#input FLD should be land masked (NaN/1) by caller if needed

#get land masks (NaN/1):
mskC=fill(1.0,FLD) + 0.0 * mask(FLD)
(mskW,mskS)=gradient(FLD,false)
mskW=fill(1.0,FLD) + 0.0 * mask(mskW)
mskS=fill(1.0,FLD) + 0.0 * mask(mskS)

#replace NaN with 0. in FLD and land masks:
FLD=mask(FLD,0.0)
mskC=mask(mskC,0.0)
mskW=mask(mskW,0.0)
mskS=mask(mskS,0.0)

#get inverse grid spacing:
iDXC=similar(FLD)
iDYC=similar(FLD)
for a=1:FLD.nFaces;
  iDXC.f[a]=1.0./MeshArrays.DXC.f[a]
  iDYC.f[a]=1.0./MeshArrays.DYC.f[a]
end

#Before scaling the diffusive operator ...
tmp0=DXCsm*iDXC*mskW;
tmp00=maximum(tmp0);
tmp0=DYCsm*iDYC*mskS;
tmp00=max(tmp00,maximum(tmp0));

#... determine a suitable time period:
nbt=ceil(1.1*2*tmp00^2);
dt=1.;
T=nbt*dt;
#println("nbt="*"$nbt")

#diffusion operator times DYG / DXG:
KuxFac=mskW*DXCsm*DXCsm/T/2.0*MeshArrays.DYG;
KvyFac=mskS*DYCsm*DYCsm/T/2.0*MeshArrays.DXG;

#time steping factor:
dtFac=dt*mskC/MeshArrays.RAC;

#loop:
for it=1:nbt
  (dTdxAtU,dTdyAtV)=gradient(FLD,iDXC,iDYC);
  tmpU=similar(FLD)
  tmpV=similar(FLD)
  for a=1:FLD.nFaces
      tmpU.f[a]=dTdxAtU.f[a].*KuxFac.f[a];
      tmpV.f[a]=dTdyAtV.f[a].*KvyFac.f[a];
  end
  tmpC=convergence(tmpU,tmpV);
  for a=1:FLD.nFaces
      FLD.f[a]=FLD.f[a]-dtFac.f[a].*tmpC.f[a];
  end
end

#Apply land mask (NaN/1) to end result:
mskC=mask(mskC,NaN,0.0)
FLD=mskC*FLD

return FLD

end

##


text:
using Test
using PSID, DataDeps, JSON3

@show pwd()
x = dirname(pathof(PSID))
fx = "$x/allfiles_hash.json"
skipdata = try
    PSID.verifyfiles(fx)
    println("Found all files, running full tests")
    false
catch
    println("Did not find data files, running partial tests")
    true
end

if skipdata
    Base.download("https://raw.githubusercontent.com/aaowens/PSID.jl/master/examples/user_input.json", "user_input.json")
    Base.download("https://drive.google.com/uc?authuser=0&id=1nz1UaVGcj0ur2Bp3ev7a8agJbj0A5JTF&export=download", "J265684_codebook.zip")
    run(DataDeps.unpack_cmd("J265684_codebook.zip", "$(pwd())", ".zip", ""))
    Base.download("https://psidonline.isr.umich.edu/help/xyr/psid.xlsx", "psid.xlsx")
    userinput_json = "user_input.json"
    isfile(userinput_json) || error("$userinput_json not found in current directory")
    isdir("output") || mkdir("output")
    isdir("datafiles") || mkdir("datafiles")
    PSID.process_codebook()
    PSID.process_input("user_input.json")
    JSON3.read(read("output/user_output.json", String), Vector{PSID.VarInfo5})
    #famdatas, inddata = PSID.unzip_data()
    #PSID.construct_alldata(famdatas, inddata)
else
    makePSID("user_input.json")
end


text:
# ---
# title: 35. Search Insert Position
# id: problem35
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Array, Binary Search
# link: <https://leetcode.com/problems/search-insert-position/description/>
# hidden: true
# ---
#
# Given a sorted array of distinct integers and a target value, return the index
# if the target is found. If not, return the index where it would be if it were
# inserted in order.
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [1,3,5,6], target = 5
#     Output: 2
#
#
# **Example 2:**
#
#
#
#     Input: nums = [1,3,5,6], target = 2
#     Output: 1
#
#
# **Example 3:**
#
#
#
#     Input: nums = [1,3,5,6], target = 7
#     Output: 4
#
#
# **Example 4:**
#
#
#
#     Input: nums = [1,3,5,6], target = 0
#     Output: 0
#
#
# **Example 5:**
#
#
#
#     Input: nums = [1], target = 0
#     Output: 0
#
#
#
#
# **Constraints:**
#
#   * `1 <= nums.length <= 104`
#   * `-104 <= nums[i] <= 104`
#   * `nums` contains **distinct** values sorted in **ascending** order.
#   * `-104 <= target <= 104`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
# Set model up
m = Model990(testing = true)
m <= Setting(:alternative_policy, DSGE.taylor99())
Γ0, Γ1, C, Ψ, Π = DSGE.taylor99_eqcond(m)

# Check policy rule
eq = m.equilibrium_conditions
endo = m.endogenous_states
@testset "Check Taylor99 rule has the right coefficients." begin
    @test Γ0[eq[:eq_mp], endo[:R_t]] == 1
    @test Γ0[eq[:eq_mp], endo[:π_a_t]] == -1.5/4
    @test Γ0[eq[:eq_mp], endo[:y_t]]  == -1/4
    @test Γ0[eq[:eq_mp], endo[:y_f_t]] == 1/4
end


text:
function DormandPrince853Integrator(arg0::jdouble, arg1::jdouble, arg2::Vector{jdouble}, arg3::Vector{jdouble})
    return DormandPrince853Integrator((jdouble, jdouble, Vector{jdouble}, Vector{jdouble}), arg0, arg1, arg2, arg3)
end

function DormandPrince853Integrator(arg0::jdouble, arg1::jdouble, arg2::jdouble, arg3::jdouble)
    return DormandPrince853Integrator((jdouble, jdouble, jdouble, jdouble), arg0, arg1, arg2, arg3)
end

function get_a(obj::DormandPrince853Integrator)
    return jcall(obj, "getA", Vector{Vector{jdouble}}, ())
end

function get_b(obj::DormandPrince853Integrator)
    return jcall(obj, "getB", Vector{jdouble}, ())
end

function get_c(obj::DormandPrince853Integrator)
    return jcall(obj, "getC", Vector{jdouble}, ())
end

function get_order(obj::DormandPrince853Integrator)
    return jcall(obj, "getOrder", jint, ())
end


text:
function KernelMatrix(args...;UniFrac=theUniFrac)
	N = size(UniFrac)[1]
	D2 = UniFrac .^2

    I0 = Array{Float64}(N,N)
    fill!(I0, 1.0)
    I = Diagonal(I0)

    one = Array{Float64}(N,1)
    fill!(one, 1.0)
    one_t = transpose(one)

    K = -0.5 * (I - one .* one_t / N) * D2 * (I - one .* one_t / N)
    #U = svdfact(K)[:U]
    #S = svdfact(K)[:S]
    #Vt = svdfact(K)[:Vt]
    # psdK = U .* abs(S) * Vt
    evecK = eigvecs(K,permute=false,scale=true)
    #writedlm("dd",evecK)
    #pwd()
    #convert(Array{Float64, 2}, evecK)
    evalK = abs.(eigvals(K))
		evalK[evalK.<1e-4]=1e-4
    #orderK= sort!(evalK)[1:2]
    #evalK[evalK .==orderK[1]] = min(1e-4,orderK[2]/2)
    diagK = Diagonal(evalK)
    psdK = evecK * diagK * transpose(evecK)

    #mineig=minimum(abs(eigvals(K)))
    #psdK=psdK-mineig*I
    #psdK=convert(Array(Float64,2), psdK)
    #cholfact(psdK,[pivot=Val{true}])
    return psdK
    # return K
end


text:
using Printf                    #

"""
    mutable struct LinkedInts
        next :: Vector{Int}
        prev :: Vector{Int}
        back :: Vector{Int}

        free_ptr :: Int
        free_cap :: Int
        root     :: Int

        block :: Vector{Int}
        size  :: Vector{Int}
    end

Block linked list, there are `length(next) = length(prev)` indices, `free_cap`
of which are free and the other are used. `root` gives the last one ones, i.e.
starting from `root` and following `prev` until it gives zero should give the
`length(prev) - free_cap` used indices. Similarly, `free_ptr` gives the last
index free, i.e. starting from `free_ptr` and following `prev` until it gives
zero should give the `free_cap` free indices.

`length(next) == length(prev) == length(back)`.
* `next[i]` is the minimum `j > i` such that `j` is used, or 0 if `i` is the
  last one used, i.e. `i = root`.
* `prev[i]` is the maximum `j < i` such that `j` is used, or 0 if `i` is the
  first one used.
* `back[i]` is such that `block[back[i]]` is the first index of the block
  in which `i` is.

* `free_cap` is the number of free slots.
* `free_ptr` last index of the `free_cap` free slots.
* `root` last used index.

`length(block) == length(size)`.
* `block`: mapping from block index to first index of the block
* `size`: mapping from block index to length of the block
"""
mutable struct LinkedInts
    next :: Vector{Int}
    prev :: Vector{Int}
    back :: Vector{Int}

    free_ptr :: Int
    free_cap :: Int
    root     :: Int

    block :: Vector{Int}
    size  :: Vector{Int}
end

function LinkedInts(capacity=128)
    return LinkedInts(Int[], Int[], Int[],
                      0, 0, 0,
                      Int[], Int[])
end

num_allocated(s::LinkedInts) = count(b -> b > 0, s.block)
allocatedlist(s::LinkedInts) = findall(s.block .> 0)
allocated(s::LinkedInts, id :: Int) = id > 0 && id <= length(s.block) && s.block[id] > 0
blocksize(s::LinkedInts, id :: Int) = s.size[id]
Base.length(s::LinkedInts) = length(s.next)
numblocks(s::LinkedInts) = length(s.block)

function Base.show(f::IO, s::LinkedInts)
    print(f,"LinkedInts(\n")
    @printf(f,"  Number of blocks: %d\n", length(s.block))
    @printf(f,"  Number of elements: %d\n", length(s.next))
    print(f,"  Blocks:\n")
    for i in 1:length(s.block)
        if s.block[i] > 0
            idxs = getindexes(s,i)
            print(f,"    #$i: $idxs\n")
        end
    end
    p = s.free_ptr
    freelst = Int[]
    while p > 0
        push!(freelst,p)
        p = s.prev[p]
    end
    print(f,"  Free: $freelst\n")

    println(f,"  free_ptr = $(s.free_ptr)")
    println(f,"  root     = $(s.root)")
    println(f,"  next     = $(s.next)")
    println(f,"  prev     = $(s.prev)")

    print(f,")")
end

"""
    ensurefree(s::LinkedInts, N :: Int)

Ensure that there are at least `N` elements free, and allocate as necessary.
"""
function ensurefree(s::LinkedInts, N :: Int)
    if s.free_cap < N
        num = N - s.free_cap

        cap = length(s.next)
        first = cap+1
        last  = cap+num

        append!(s.next, Int[i+1 for i in first:last])
        append!(s.prev, Int[i-1 for i in first:last])
        append!(s.back, zeros(Int, num))

        s.next[last] = 0
        s.prev[first] = s.free_ptr
        if s.prev[first] > 0
            s.next[s.prev[first]] = first
        end
        s.free_ptr = last
        s.free_cap += num

        return num
    else
        return 0
    end
end

function allocate_block(s::LinkedInts, N::Int, id::Integer)
    @assert(N > 0)
    ensurefree(s, N)
    # remove from free list
    ptre = s.free_ptr
    # ptre is the last index
    ptrb = ptre
    for i = 1:N-1
        s.back[ptrb] = id
        ptrb = s.prev[ptrb]
    end
    s.back[ptrb] = id
    s.block[id] = ptrb
    # ptrb is the first index

    prev = s.prev[ptrb]

    if prev > 0
        s.next[prev] = 0
    end

    s.free_ptr = s.prev[ptrb]
    s.free_cap -= N

    # insert into list `idx`
    s.prev[ptrb] = s.root
    if s.root > 0
        s.next[s.root] = ptrb
    end
    s.root = ptre

    #if ! checkconsistency(s)
    #    println("List = ",s)
    #    assert(false)
    #end

    return id
end

function create_block(s::LinkedInts, N::Int)
    @assert(N > 0)
    push!(s.size, N)
    push!(s.block, 0)
    @assert length(s.size) == length(s.block)
    return length(s.block)
end

"""
    newblock(s::LinkedInts, N :: Int)

Add a new block in list `idx`
"""
function newblock(s::LinkedInts, N::Int)::Int
    id = create_block(s, N)
    allocate_block(s, N, id)
    return id
end

"""
Move a block to the free list.
"""
function deleteblock(s::LinkedInts, id :: Int)
    if s.size[id] > 0
        ptrb = s.block[id]
        N = s.size[id]
        ptre = ptrb
        s.back[ptre] = 0
        for i in 2:N
            ptre = s.next[ptre]
            s.back[ptre] = 0
        end
        prev = s.prev[ptrb]
        next = s.next[ptre]

        # remove from list and clear the block id
        if s.root == ptre s.root = prev end
        if prev > 0 s.next[prev] = next end
        if next > 0 s.prev[next] = prev end

        s.size[id]  = 0
        s.block[id] = 0

        # add to free list
        if s.free_ptr > 0
            s.next[s.free_ptr] = ptrb
        end
        s.prev[ptrb] = s.free_ptr
        s.free_ptr = ptre
        s.next[ptre] = 0

        s.free_cap += N
    end
end

# TODO merge getoneindex and getindex
function getoneindex(s::LinkedInts, id :: Int)
    N = s.size[id]
    if N < 1
        error("No values at id")
    end

    s.block[i]
end

"""
    getindex(s::LinkedInts, id::Int)

Shortcut for `getindexes(s, id)[1]` when `s.size[id]` is 1.
"""
function getindex(s::LinkedInts, id::Int)
    @assert s.size[id] == 1
    @assert s.back[s.block[id]] == id
    return s.block[id]
end

"""
    getindexes(s::LinkedInts, id :: Int)

Return the vector of indices for the block `id`.
"""
function getindexes(s::LinkedInts, id :: Int)
    @assert length(s.next) == length(s.prev) == length(s.back)
    N = s.size[id]
    r = Vector{Int}(undef, N)
    p = s.block[id]
    for i in 1:N
        @assert s.back[p] == id
        r[i] = p
        p = s.next[p]
    end
    return r
end

function getindexes(s::LinkedInts, id::Int, target::Vector{Int}, offset::Int)
    N = s.size[id]
    p = s.block[id]
    for i in 1:N
        @assert s.back[p] == id
        target[i+offset-1] = p
        p = s.next[p]
    end
    return N
end

function getindexes(s::LinkedInts, ids::Vector{Int})
    N = sum(map(id -> s.size[id], ids))
    r = Vector{Int}(undef,N)
    offset = 1
    for id in ids
        offset += getindexes(s, id, r, offset)
    end
    return r
end

"""
Get a list if the currently free elements.
"""
function getfreeindexes(s::LinkedInts)
    N = s.free_cap
    r = Array{Int}(undef,N)
    ptr = s.free_ptr
    for i in 1:N
        @assert iszero(s.back[ptr])
        r[N-i+1] = ptr
        ptr  = s.prev[ptr]
    end
    r
end

"""
Get a list if the currently used elements.
"""
function getusedindexes(s::LinkedInts)
    N = length(s.next) - s.free_cap
    r = Array{Int}(undef,N)
    ptr = s.root
    for i in 1:N
        @assert !iszero(s.back[ptr])
        r[N-i+1] = ptr
        ptr  = s.prev[ptr]
    end
    r
end

"""
Check consistency of the internal structures.
"""
function checkconsistency(s::LinkedInts) :: Bool
    if length(s.prev) != length(s.next)
        return false
    end

    N = length(s.prev)

    if ! (all(i -> s.prev[i] == 0 || s.next[s.prev[i]] == i, 1:N) &&
          all(i -> s.next[i] == 0 || s.prev[s.next[i]] == i, 1:N))
        @assert(false)
    end

    mark = fill(false,length(s.prev))

    p = s.free_ptr
    while p != 0
        @assert iszero(s.back[ptr])
        mark[p] = true
        p = s.prev[p]
    end

    p = s.root
    while p != 0
        @assert(!mark[p])
        @assert !iszero(s.back[ptr])
        mark[p] = true
        p = s.prev[p]
    end

    if !all(mark)
        println(s)
        println(mark)
        @assert(all(mark))
    end

    return true
end


text:
#. PREAMBLE
using StaticArrays
using Random
using Singular

import AbstractAlgebra
import Nemo

const AA = AbstractAlgebra

#. MONOMIALS

#.. Abstract context

""" Abstract type for monomial contexts
Must provide:
    - termorder(::MonomialContext{T}) :: TermOrder
    - exponents(::MonomialContext{T}, ::T) :: SVector
    - (::MonomialContext{T})(::SVector)

It is probably interesting to redefine:
    - mul(::MonomialContext{T}, ::T, ::T) :: T
    - div(::MonomialContext{T}, ::T, ::T) :: T
    - lcm(::MonomialContext{T}, ::T, ::T) :: T
    - gcd(::MonomialContext{T}, ::T, ::T) :: T
"""
abstract type MonomialContext{T} <: Context{T} end

ordervector(ctx::Γ, x::T) where {T, Γ<:MonomialContext{T}} =
    ordervector(termorder(ctx), exponents(ctx, x))

divides(ctx::Γ, x::T, y::T) where {T, Γ<:MonomialContext{T}} =
    divides(termorder(ctx), exponents(ctx, x), exponents(ctx, y))

iscompatible(ctx::Γ, x::T, y::T) where {T, Γ<:MonomialContext{T}} =
    iscompatible(termorder(ctx), exponents(ctx, x), exponents(ctx, y))

lt(ctx::Γ, x::T, y::T) where {T, Γ<:MonomialContext{T}} =
    lt(termorder(ctx), exponents(ctx, x), exponents(ctx, y))

function mul(ctx::MonomialContext{T}, x::T, y::T) where {T}
    ctx(exponents(ctx, x)+exponents(ctx, y))
end

function div(ctx::MonomialContext{T}, x::T, y::T) where {T}
    ctx(exponents(ctx, x)-exponents(ctx, y))
end

function lcm(ctx::MonomialContext{T}, x::T, y::T) where {T}
    ctx(max.(exponents(ctx, x), exponents(ctx, y)))
end

function gcd(ctx::MonomialContext{T}, x::T, y::T) where {T}
    ctx(min.(exponents(ctx, x), exponents(ctx, y)))
end

#.. Monomial type

@generated function makehash(::Val{N}, m) where {N}
    rng = MersenneTwister(18121987)
    hash = :( $(UInt64(0)) )
    for i in 1:N
        hash = :($hash + $(rand(rng, UInt))*UInt(m[$i]))
    end
    return hash
end

makehash(m::SVector{N}) where N = makehash(Val(N), m)

# Monomial in N variables with exponent type E
# We require E to be signed because it makes grevlex easier.
# [TODO is that a good reason?]
struct Monomial{N, E<:Signed}
    exponents::SVector{N, E}
    hash :: UInt

    Monomial{N, E}(e, h :: UInt) where {N, E} = new(e, h)
    Monomial{N, E}(e) where {N, E} = new(e, makehash(Val{N}(), e))
    Monomial{N, E}(m :: Monomial) where {N, E} = new(m.exponents, m.hash)
end

Monomial(e) = Monomial{length(e), eltype(e)}(e, makehash(e))

Base.hash(m::Monomial) = m.hash
Base.show(io::IO, m::Monomial) = Base.show(io, convert(Vector{Int}, m.exponents))
Base.isequal(x::Monomial, y::Monomial) = x.exponents == y.exponents

# Compute a bitmask for the monomial M, TODO: make this generated
function bitmask(m::Monomial{N, E}, pwrs::Dict{Int, Vector{E}}; masktype = UInt32) where {N, E<:Signed}
    @inbounds arr = BitArray([p <= m.exponents[i] for i in 1:N for p in pwrs[i]])
    reverse(arr).chunks[1] % masktype
end

Base.one(::Type{Monomial{N, E}}) where {N, E} = Monomial{N, E}(SVector{N, E}(zeros(E, N)))
Base.isone(m::Monomial{N, E}) where {N, E} = m == one(typeof(m))

#.. Monomial context

struct MonomialΓ{N, E, O<:TermOrder{N}}<:MonomialContext{Monomial{N, E}}
    order::O
    vars::Vector{String}
end

# Int16 exponents allows for degree 32727, this leaves room for many applications
function monomialctx(;exponents=Int16, order::TermOrder, vars=[], varprefix="x", indexed=false, kwargs...)
    if isempty(vars)
        vars = [varprefix*string(i) for i in 1:nvars(order)]
    end
    mo = MonomialΓ{nvars(order), exponents, typeof(order)}(order, vars)
    if indexed
        return ixmonomialctx(mo; kwargs...)
    else
        return mo
    end
end

params(::MonomialΓ{N, E}) where {N, E} = N, E
exponents(::MonomialΓ{N, E}, m :: Monomial{N, E}) where {N, E}= m.exponents
exponenttype(::MonomialΓ{N, E}) where {N, E} = E
nvars(::MonomialΓ{N}) where N = N
variables(ctx::MonomialΓ) = ctx.vars
termorder(ctx::MonomialΓ) = ctx.order
Base.one(ctx::MonomialΓ) = Base.one(eltype(ctx))

function mul(::MonomialContext{T}, x::T, y::T) where {T<:Monomial}
    T(x.exponents + y.exponents, x.hash + y.hash)
end

function div(::MonomialContext{T}, x::T, y::T) where {T<:Monomial}
    T(x.exponents - y.exponents, x.hash - y.hash)
end

#... AA interoperability

function (ctx :: MonomialContext)(p :: AA.MPolyElem)
    AA.isterm(p) || error("Not a monomial")
    return ctx(AA.exponent_vector(p, 1))
end

function (R :: AA.MPolyRing)(ctx::MonomialContext{M}, m::M) where M
    br = AA.base_ring(R)

    # there may be a shortcut in AA, I haven't looked
    return R([one(AA.base_ring(R))],
             [convert(Vector{Int}, exponents(ctx, m))])
end

#. POLYNOMIALS

#.. Polynomial type

# The AbstractVector interface is useful, for example, to sort the terms
mutable struct Polynomial{M, T}<:AbstractVector{Tuple{M, T}}
    mo::Vector{M}
    co::Vector{T}
end

Base.size(p::Polynomial) = size(p.co)
Base.@propagate_inbounds Base.getindex(p::Polynomial, i...) = (getindex(p.mo, i...), getindex(p.co, i...))
# may not be ideal
Base.@propagate_inbounds Base.getindex(p::Polynomial, r::UnitRange) = [(p.mo[i], p.co[i]) for i in r]
Base.@propagate_inbounds function Base.setindex!(p::Polynomial{M, T}, t::Tuple{M, T}, i...) where {M, T}
    setindex!(p.mo, t[1], i...)
    setindex!(p.co, t[2], i...)
end
Base.Sort.defalg(p::Polynomial) =  Base.DEFAULT_UNSTABLE

Base.isempty(p::Polynomial) = isempty(p.co)
Base.@propagate_inbounds coefficient(p::Polynomial, i) = p.co[i]
Base.@propagate_inbounds monomial(p::Polynomial, i) = p.mo[i]
Base.@propagate_inbounds leadingmonomial(p::Polynomial) = monomial(p, 1)
Base.@propagate_inbounds function tail(p::Polynomial)
    if length(p.mo) > 1
        typeof(p)(p.mo[2:end], p.co[2:end])
    else
        zero(p)
    end
end
monomials(p::Polynomial) = p.mo
Base.zero(::Type{Polynomial{M, T}}) where {M, T} = Polynomial(M[], T[])
Base.zero(p::Polynomial) = zero(typeof(p))
Base.iszero(p::Polynomial) = p == zero(p)

ismonic(p::Polynomial) = !isempty(p) && @inbounds isone(coefficient(p, 1))

Base.copy(p::Polynomial) = typeof(p)(copy(p.mo), copy(p.co))

#.. Context

struct PolynomialΓ{M, T, MΓ<:Context{M}, TΓ<:Context{T}}<:Context{Polynomial{M, T}}
    mo::MΓ
    co::TΓ
end

function polynomialctx(coefficients; monomials=nothing, kwargs...)
    if isnothing(monomials)
        monomials = monomialctx(;kwargs...)
    end
    return PolynomialΓ{eltype(monomials), eltype(coefficients), typeof(monomials),
                       typeof(coefficients)}(monomials, coefficients)
end

nvars(ctx::PolynomialΓ) = nvars(ctx.mo)
variables(ctx::PolynomialΓ) = variables(ctx.mo)
termorder(ctx::PolynomialΓ) = termorder(ctx.mo)

function degree(ctx::PolynomialΓ{M}, p::Polynomial{M}) where M
    maximum(m -> degree(ctx.mo, m), p.mo)
end

function lt(ctx::PolynomialΓ, a, b)
    lt(ctx.mo, leadingmonomial(a), leadingmonomial(b))
end

function lt(ctx::Γ, a::Tuple{M, T}, b :: Tuple{M, T}) where {M, T, Γ<:PolynomialΓ{M, T}}
    lt(ctx.mo, a[1], b[1])
end

function normalize!(ctx::Γ, p::T; onlysort::Bool=false) where {T<:Polynomial, Γ<:Context{T}}
    sort!(p, by = mon_coeff -> first(mon_coeff), rev=true, order=order(ctx.mo))
    onlysort && return

    # collapse like terms, remove zeros
    w = 0                   # read and write indices
    r = 1
    lenp = length(p)
    while r <= lenp
        m = monomial(p, r)
        c = coefficient(p, r)

        while true
            r += 1
            if r > lenp || monomial(p, r) != m
                break
            end
            c = add(ctx.co, c, coefficient(p, r))
        end

        if !iszero(c)
            w += 1
            p[w] = (m, c)
        end
    end

    resize!(p.co, w)
    resize!(p.mo, w)
    sizehint!(p.co, w)
    sizehint!(p.mo, w)

    return
end

function monic!(ctx::TΓ,
                p::Polynomial{M, T}) where {M, T, TΓ <: Context{T}}
    isempty(p) && return
    @inbounds begin
        isone(coefficient(p, 1)) && return

        mult = inv(ctx, coefficient(p, 1))
        for (i, c) in enumerate(p.co)
            p.co[i] = mul(ctx, mult, c)
        end
    end
    return
end

#.. Basic operations

function selectcoefficient(ctx::Γ, p::P, var::Int, exp::Int) where {P <: Polynomial, Γ<:Context{P}}
    newp = ctx([], [])
    for (m, c) in p
        if exponents(ctx.mo, m)[var] == exp
            if exp == 0
                push!(newp.mo, m)
            else
                newm = ctx.mo(setindex(exponents(ctx.mo, m), 0, var))
                push!(newp.mo, newm)
            end
            push!(newp.co, c)
        end
    end
    return newp
end

function add!(ctx::Γ, p::P, q::P) where {P<:Polynomial, Γ<:Context{P}}
    append!(p.co, q.co)
    append!(p.mo, q.mo)
    normalize!(ctx, p)
    return
end

function add(ctx::Γ, p::P, q::P) where {P<:Polynomial, Γ<:Context{P}}
    cp = copy(p)
    add!(ctx, cp, q)
    return cp
end

function random_lin_comb(ctx::PolynomialΓ{M, T},
                         ps::Vector{Polynomial{M, T}}) where {M, T}

    reduce((p1, p2) -> add(ctx, mul_scalar(ctx, p1, rand_elem(ctx.co)),
                           mul_scalar(ctx, p2, rand_elem(ctx.co))), ps)
end

# Share the coefficient array
function mul(ctx::PolynomialΓ{M, T}, p::Polynomial{M, T}, m::M) where {M, T}
    newmo = [mul(ctx.mo, m, mp) for mp in p.mo]
    eltype(ctx)(newmo, p.co)
end

function mul_scalar(ctx::PolynomialΓ{M, T}, p::Polynomial{M, T}, c::T) where {M, T}
    newco = [mul(ctx.co, c, coeff) for coeff in p.co]
    eltype(ctx)(p.mo, newco)
end

#.. AA/Singular interoperability

abstractalgebra(ctx :: PolynomialΓ) = AA.PolynomialRing(abstractalgebra(ctx.co), variables(ctx))

function (ctx :: MonomialΓ)(m ::AA.MPolyElem)
    exp = first(exponent_vectors(m))
    ctx(Monomial{nvars(ctx), exponenttype(ctx)}(exp))
end

function (ctx :: PolynomialΓ)(p :: AA.MPolyElem)
    mo = [ctx.mo(m) for m in AA.monomials(p)]
    if typeof(p) <: Singular.spoly
        co = [ctx.co(Int(c)) for c in Singular.coefficients(p)]
    else
        co = [ctx.co(AA.coeff(p, i)) for i in 1:length(p)]
    end
    p = ctx(mo, co)
    sort!(p, rev=true, order=order(ctx))
    return p
end

function (R :: AA.MPolyRing)(ctx, p :: Polynomial)
    br = AA.base_ring(R)
    if typeof(R) <: Singular.PolyRing
        C = MPolyBuildCtx(R)
        for (m, c) in p
            push_term!(C, br(Int(c)), Vector{Int}(exponents(ctx.mo, m)))
        end
        finish(C)
    else
        return R([br(coefficient(p, i)) for i in 1:length(p)],
                 [convert(Vector{Int}, exponents(ctx.mo, monomial(p, i))) for i in 1:length(p)])
    end
end


text:
# Like serialize, but will also make a directory if needed.
function save(file, x)
    dir = dirname(file)
    isdir(dir) || mkpath(dir)
    serialize(file, x)
    return nothing
end

function transfer(pipe, @nospecialize(x))
    serialize(TRANSFERPATH, x)
    println(pipe, "measurements")
    return nothing
end

function paramtransfer(pipe, @nospecialize(x))
    serialize(PARAMPATH, x)
    println(pipe, "params")
    return nothing
end

# Processing pipeline for name formatting
modify(x) = x
modify(s::Union{String,Symbol}) = replace(String(s), "-"=>"_")
modify(::Val{T}) where {T} = T
modify(x::LFSR) = ceil(Int, log2(length(x)))

# If passed a function, get the name of the function
modify(f::Function) = string(last(split(string(f), ".")))

function make_params(f, nt::NamedTuple{names}) where {names}
    return NamedTuple{(:benchmark, names...)}(modify.((f, nt...)))
end

getnames(::Type{<:NamedTuple{names}}) where {names} = names

#####
##### Custom Threading
#####

# Break up an array into a bunch of views and distribute them across threads.
#
# This makes sure that things are getting split up how we expect them to be.
function threadme(f, A, args...; prepare = false, iterations = 1)
    nthreads = Threads.nthreads()
    @assert iszero(mod(length(A), nthreads))

    step = div(length(A), nthreads)
    Threads.@threads for i in 1:Threads.nthreads()
        threadid = Threads.threadid()
        start = step * (i-1) + 1
        stop = step * i
        x = view(A, start:stop)

        # Run the inner loop
        for j in 1:iterations
            f(x, args...)
        end
    end
    return nothing
end

max_stars_repo_name: darchr/AutoTM

max_issues_repo_name: darchr/AutoTM

max_forks_repo_name: darchr/AutoTM


text:
using ReinforcementLearning
using StableRNGs
using Flux
using Flux.Losses

function RL.Experiment(
    ::Val{:JuliaRL},
    ::Val{:BCQ},
    ::Val{:Pendulum},
    type::AbstractString;
    save_dir = nothing,
    seed = 123,
)
    rng = StableRNG(seed)
    inner_env = PendulumEnv(T = Float32, rng = rng)
    A = action_space(inner_env)
    low = A.left
    high = A.right
    ns = length(state(inner_env))
    na = 1
    latent_dims = 10

    trajectory_num = 10000
    dataset_size = 10000
    batch_size = 64

    env = ActionTransformedEnv(
        inner_env;
        action_mapping = x -> low + (x[1] + 1) * 0.5 * (high - low),
    )
    init = glorot_uniform(rng)

    create_policy_net() = NeuralNetworkApproximator(
        model = PerturbationNetwork(
            base = Chain(
                Dense(ns + na, 64, relu; init = glorot_uniform(rng)),
                Dense(64, 64, relu; init = glorot_uniform(rng)),
                Dense(64, na; init = glorot_uniform(rng))
            ),
            ϕ = 0.05f0,
        ),
        optimizer = ADAM(0.003),
    )

    create_q_net() = NeuralNetworkApproximator(
        model = Chain(
            Dense(ns + na, 64, relu; init = init),
            Dense(64, 64, relu; init = init),
            Dense(64, 1; init = init),
        ),
        optimizer = ADAM(0.003),
    )

    create_vae_net() = NeuralNetworkApproximator(
        model = VAE(
            encoder = GaussianNetwork(
                pre = Chain(
                    Dense(ns + na, 64, relu),
                    Dense(64, 64, relu),
                ),
                μ = Chain(Dense(64, latent_dims, init = init)),
                logσ = Chain(Dense(64, latent_dims, init = init)),
            ),
            decoder = Chain(
                Dense(ns + latent_dims, 64, relu; init = init),
                Dense(64, 64, relu; init = init),
                Dense(64, na; init = init),
            ),
            latent_dims = latent_dims,
        ),
        optimizer = ADAM(0.003),
    )

    agent = Agent(
        policy = OfflinePolicy(
            learner = BCQLearner(
                policy = create_policy_net() |> cpu,
                target_policy = create_policy_net() |> cpu,
                qnetwork1 = create_q_net() |> cpu,
                qnetwork2 = create_q_net() |> cpu,
                target_qnetwork1 = create_q_net() |> cpu,
                target_qnetwork2 = create_q_net() |> cpu,
                vae = create_vae_net() |> cpu,
                γ = 0.99f0,
                τ = 0.005f0,
                λ = 0.75f0,
                p = 10,
                batch_size = batch_size,
                start_step = 1000,
                update_freq = 1,
            ),
            dataset = gen_JuliaRL_dataset(:SAC, :Pendulum, type; dataset_size = dataset_size),
            continuous = true,
            batch_size = batch_size,
        ),
        trajectory = CircularArraySARTTrajectory(
            capacity = trajectory_num,
            state = Vector{Float32} => (ns,),
            action = Vector{Float32} => (na,),
        ),
    )

    stop_condition = StopAfterStep(trajectory_num, is_show_progress=!haskey(ENV, "CI"))
    hook = TotalRewardPerEpisode()
    Experiment(agent, env, stop_condition, hook, "BCQ <-> Pendulum ($type dataset)")
end

using Plots
pyplot() #hide
ex = E`JuliaRL_BCQ_Pendulum(medium)`
run(ex)
plot(ex.hook.rewards)
savefig("assets/JuliaRL_BCQ_Pendulum_medium.png") #hide

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl