text:
struct Track
    x::Vector{Float64}
    y::Vector{Float64}
    lane_width::Vector{Float64}
    x′::Vector{Float64}
    y′::Vector{Float64}
    lane_width′::Vector{Float64}
    boundary_left::Matrix{Float64}
    boundary_right::Matrix{Float64}
    sample_factor::Int
end

function Track(infile::String, width::Vector{Float64}; sample_factor=20)
    csv_data = CSV.File(infile, header=false)
    csv_data.cols == 2 || error("Can only have 2 columns for a track file")
    length(csv_data.Column1) == length(width) || error("Supplied width vector does not match length of track file")
    x = csv_data.Column1
    y = csv_data.Column2
    lane_width = width
    x′ = x[1:sample_factor:end]
    y′ = y[1:sample_factor:end]
    lane_width′ = lane_width[1:sample_factor:end]

    l_boundary, r_boundary = calculate_boundary(x, y, lane_width)

    return Track(x, y, lane_width, x′, y′, lane_width′, l_boundary, r_boundary, sample_factor)
end

function Track(infile::String; width::Float64=15.0, sample_factor=20)
    csv_data = CSV.File(infile, header=false)
    lane_width = ones(Float64, length(csv_data)) * width
    return Track(infile, lane_width, sample_factor=sample_factor)
end

function calculate_boundary(x::Vector{Float64}, y::Vector{Float64}, w::Vector{Float64})

    r_bound = zeros(Float64, length(x), 2)
    l_bound = zeros(Float64, length(x), 2)

    Δx = x[2] - x[1]
    Δy = y[2] - y[1]
    p = [-Δy, Δx] ./ norm([-Δy, Δx])
    l_bound[1,:] = [x[1], y[1]] + w[1]*p
    r_bound[1,:] = [x[1], y[1]] - w[1]*p

    for ii in 2:(length(x)-1)
        Δx = x[ii+1] - x[ii-1]
        Δy = y[ii+1] - y[ii-1]
        p = [-Δy, Δx] ./ norm([-Δy, Δx])
        l_bound[ii,:] = [x[ii], y[ii]] + w[ii]*p
        r_bound[ii,:] = [x[ii], y[ii]] - w[ii]*p
    end

    Δx = x[end] - x[end-1]
    Δy = y[end] - y[end-1]
    p = [-Δy, Δx] ./ norm([-Δy, Δx])
    l_bound[end,:] = [x[end], y[end]] + w[end]*p
    r_bound[end,:] = [x[end], y[end]] - w[end]*p

    return l_bound, r_bound
end

function within_track(track::Track, pos::Vector{Int})
    return within_track(track::Track, Float64.(pos))
end

function within_track(track::Track, pos::Vector{Float64})
    length(pos) == 2 || error("Position is only 2D")

    dists = (track.x′ .- pos[1]).^2 .+ (track.y′ .- pos[2]).^2

    _ , min_idx  = findmin(dists)
    min_idx = Int(min_idx)
    min_idx_m1 = mod1(min_idx - 1, length(track.x′))
    min_idx_p1 = mod1(min_idx + 1, length(track.x′))
    dist_m1 = norm([track.x′[min_idx_m1], track.y′[min_idx_m1]] - pos)
    dist_p1 = norm([track.x′[min_idx_p1], track.y′[min_idx_p1]] - pos)
    min_idx_2 = dist_m1 <= dist_p1 ? min_idx_m1 : min_idx_p1

    # Projecting the point (p3) onto the line defined by the closest two points (p1, p2)
    # The line is parametrized as p1 + t(p2 - p1)
    # The projection falls where t = [(p3-p1) ⋅ (p2-p1)] / |p2-p1|^2
    p1 = [track.x′[min_idx], track.y′[min_idx]]
    p2 = [track.x′[min_idx_2], track.y′[min_idx_2]]
    p3 = pos
    t = ((p3-p1) ⋅ (p2-p1)) /  ((p2-p1) ⋅ (p2-p1))
    projected_point = p1 + t * (p2 - p1)
    dist_to_pt = norm(projected_point - pos)
    return (within=(dist_to_pt < track.lane_width′[min_idx]), dist=dist_to_pt)

end

max_stars_repo_name: sisl/MPOPIS

max_issues_repo_name: sisl/MPOPIS

max_forks_repo_name: sisl/MPOPIS


text:
# AbstractString constants for NLEQ options
const OPT_RTOL                  = "relativeTolerance"
const OPT_QSUCC                 = "successiveCall"
const OPT_MODE                  = "mode"
const OPT_JACGEN                = "jacobianMethod"
const OPT_JACFCN                = "jacobianFunction"
const OPT_MSTOR                 = "jacobiStorage"
const OPT_ML                    = "lowerBandwidth"
const OPT_MU                    = "upperBandwidth"
const OPT_ISCAL                 = "scaling"
const OPT_PRINTWARNING          = "printWarning"
const OPT_PRINTITERATION        = "printIterationMonitor"
const OPT_PRINTIOWARN           = "printIOwarning"
const OPT_PRINTIOMON            = "printIOmonitor"
const OPT_PRINTIOSOL            = "printIOsolution"
const OPT_PRINTSOLUTION         = "printSolution"
const OPT_NONLIN                = "nonlinearType"
const OPT_QRANK1                = "rank1BroydenUpdates"
const OPT_QORDI                 = "ordinaryNewton"
const OPT_QSIMPL                = "simpleNewton"
const OPT_NOROWSCAL             = "automaticRowScaling"
const OPT_BOUNDEDDAMP           = "boundedDampingStrategy"
const OPT_IORMON                = "convergenceOrderMonitor"
const OPT_NITMAX                = "maxIter"
const OPT_FCBAND                = "boundedDampingFactor"
const OPT_SIGMA                 = "broydenDecisionParameter"
const OPT_SIGMA2                = "dampingFactorParameter"
const OPT_AJDEL                 = "relativePerturbation"
const OPT_AJMIN                 = "thresholdValue"
const OPT_ETADIF                = "targetRelativePerturbation"
const OPT_ETAINI                = "initialDenominatorDifference"
const OPT_NBROY                 = "maxBroydenSteps"
const OPT_FCSTART               = "initialDampingFactor"
const OPT_FCMIN                 = "minimumDampingFactor"
const OPT_IRANK                 = "initialRank"
const OPT_COND                  = "subConditionNumber"
const OPT_STORE                 = "storeIntermediateValues"

# Statistics
const STATS_NITER               = "numberOfNewtonIterations"
const STATS_NCORR               = "numberOfCorrectorSteps"
const STATS_NFCN                = "numberOfFunctionEvals"
const STATS_NFCNJ               = "numberOfFunctionEvalsForJac"
const STATS_NJAC                = "numberOfJacobianCalls"
const STATS_NREJR1              = "numRejectedRank1NewtonSteps"
const STATS_XSCAL               = "lastScalingVector"
const STATS_RTOL                = "finalRelativeAccuracy"
const STATS_XITER               = "allNewtonIterates"
const STATS_NATLEVEL            = "seqNaturalLevels"
const STATS_SIMLEVEL            = "seqNaturalLevelsSimpl"
const STATS_STDLEVEL            = "seqStandardLevels"
const STATS_PRECISION           = "achievedPrecision"
const STATS_DAMPINGFC           = "seqDampingFactors"
const STATS_NEW                 = "consecutiveRank1Updates"
const STATS_ICONV               = "convergenceMonitorState"
const STATS_CONV                = "conv"
const STATS_SUMX                = "sumx"
const STATS_DLEVF               = "dlevf"
const STATS_SUBCOND             = "subcondition"
const STATS_SENS                = "sensitivity"
const STATS_IFAIL               = "failureIndicator"

# Temporary Workspace variables
const WK_A                      = "a"
const WK_DXSAVE                 = "dxsave"
const WK_QA_DXSAVE              = "qa_dxsave"
const WK_DX                     = "dx"
const WK_DXQ                    = "dxq"
const WK_QA                     = "qa"
const WK_XA                     = "xa"
const WK_XWA                    = "xwa"
const WK_F                      = "f"
const WK_FA                     = "fa"
const WK_ETA                    = "eta"
const WK_XW                     = "xw"
const WK_FW                     = "fw"
const WK_DXQA                   = "dxqa"
const WK_SUMXA0                 = "sumxa0"
const WK_SUMXA1                 = "sumxa1"
const WK_QU                     = "qu"
const WK_T1                     = "t1"
const WK_T2                     = "t2"
const WK_T3                     = "t3"
const WK_FCMON                  = "fcmon"
const WK_FCA                    = "fca"
const WK_FCKEEP                 = "fckeep"
const WK_FCPRI                  = "fcpri"
const WK_DMYCOR                 = "dmycor"
const WK_SUMXS                  = "sumxs"
const WK_SENS1                  = "sens1"

"""
function getMachineConstants(T::DataType)

Get the machine constants depending on the DataType T.

Supported DataTypes are Float64 and BigFloat.
"""
function getMachineConstants(T::DataType)
    if T != Float64 && T != DataType
        println("ERROR: Wrong DataType provided. You provided $T, expected Float64 or BigFloat.")
        error("getMachineConstants - wrong DataType provided");
    end

    epMach = eps(T)
    small  = sqrt(realmin(T)/eps(T))
    great  = sqrt(realmax(T)/10.0)

    return (epMach, small, great)
end

ext: jl


text:
using TreeViews, IterTools
# Dicts
TreeViews.hastreeview(::AbstractDict) = true
TreeViews.numberofnodes(d::AbstractDict) = length(d)
TreeViews.nodelabel(io::IO, d::AbstractDict, i::Integer, ::MIME"text/plain") = print(io, nth(keys(d), i))
TreeViews.treenode(d::AbstractDict, i::Integer) = d[nth(keys(d), i)]

# Vectors
TreeViews.hastreeview(::AbstractVector) = true
TreeViews.numberofnodes(d::AbstractVector) = length(d)
TreeViews.nodelabel(io::IO, d::AbstractVector, i::Integer, ::MIME"text/plain") = print(io, "")
TreeViews.treenode(d::AbstractVector, i::Integer) = d[i]


text:
export iauFk5hz
"""
Transform an FK5 (J2000.0) star position into the system of the
Hipparcos catalogue, assuming zero Hipparcos proper motion.

This function is part of the International Astronomical Union's
SOFA (Standards Of Fundamental Astronomy) software collection.

Status:  support function.

Given:
   r5           double   FK5 RA (radians), equinox J2000.0, at date
   d5           double   FK5 Dec (radians), equinox J2000.0, at date
   date1,date2  double   TDB date (Notes 1,2)

Returned:
   rh           double   Hipparcos RA (radians)
   dh           double   Hipparcos Dec (radians)

Notes:

   1. This function converts a star position from the FK5 system to
   the Hipparcos system, in such a way that the Hipparcos proper
   motion is zero.  Because such a star has, in general, a non-zero
   proper motion in the FK5 system, the function requires the date
   at which the position in the FK5 system was determined.

   2. The TT date date1+date2 is a Julian Date, apportioned in any
   convenient way between the two arguments.  For example,
   JD(TT)=2450123.7 could be expressed in any of these ways,
   among others:

         date1          date2

      2450123.7           0.0       (JD method)
      2451545.0       -1421.3       (J2000 method)
      2400000.5       50123.2       (MJD method)
      2450123.5           0.2       (date & time method)

   The JD method is the most natural and convenient to use in
   cases where the loss of several decimal digits of resolution
   is acceptable.  The J2000 method is best matched to the way
   the argument is handled internally and will deliver the
   optimum resolution.  The MJD method and the date & time methods
   are both good compromises between resolution and convenience.

   3. The FK5 to Hipparcos transformation is modeled as a pure
   rotation and spin;  zonal errors in the FK5 catalogue are not
   taken into account.

   4. The position returned by this function is in the Hipparcos
   reference system but at date date1+date2.

   5. See also iauFk52h, iauH2fk5, iauHfk5z.

Called:
   iauS2c       spherical coordinates to unit vector
   iauFk5hip    FK5 to Hipparcos rotation and spin
   iauSxp       multiply p-vector by scalar
   iauRv2m      r-vector to r-matrix
   iauTrxp      product of transpose of r-matrix and p-vector
   iauPxp       vector product of two p-vectors
   iauC2s       p-vector to spherical
   iauAnp       normalize angle into range 0 to 2pi

Reference:

   F.Mignard & M.Froeschle, 2000, Astron.Astrophys. 354, 732-739.

This revision:  2013 June 18

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# void iauFk5hz(double r5, double d5, double date1, double date2,
#               double *rh, double *dh)

function iauFk5hz(r5::Real, d5::Real, date1::Real, date2::Real)
   ref_rh = Ref{Float64}(0.0)
   ref_dh = Ref{Float64}(0.0)

   ccall((:iauFk5hz, libsofa_c), Cvoid,
         (Cdouble, Cdouble, Cdouble, Cdouble,
         Ref{Cdouble}, Ref{Cdouble}),
          convert(Float64, r5), convert(Float64, d5),
          convert(Float64, date1), convert(Float64, date2),
          ref_rh, ref_dh)

   return ref_rh[], ref_dh[]
end


text:
# ---- generated SymmetricTensors ----

"""
  randsymarray(T, dim ::Int, N::Int = 4)

Returns N-dimmensional random super-symmetric array with elements of type T drawn from uniform distribution on [0,1),
dim denotes data size.

"""
function randsymarray(::Type{T}, dim::Int, N::Int = 4) where T<:Number
  t = zeros(fill(dim, N)...,)
  for i in pyramidindices(N,dim)
    n = rand(T)
    for j in collect(permutations(i))
        @inbounds t[j...] = n
    end
  end
  t
end

"""
  randsymarray(T, dim::Int, N::Int = 4)

Returns N-dimmensional random super-symmetric array with Float64 elements drawn from uniform distribution on [0,1),
dim denotes data size.

"""
randsymarray(dim::Int, N::Int = 4) = randsymarray(Float64, dim, N)

"""

  fixpointperms(j::NTuple{N, Int}) where N

Returns Vactor{Vector}, a fix point permutation of given multiindex

TODO: this is a naive implementation
"""
function fixpointperms(j::NTuple{N, Int}) where N
  r = []
  for p in permutations(1:N)
    if j == j[p]
      push!(r, p)
    end
  end
  r
end

"""
  randblock(::Type{T}, dims::NTuple{N, Int}, j::NTuple{N, Int})

Returns a block of size dims and position j by a uniformly distributed random number
of type T

"""
function randblock(::Type{T}, dims::NTuple{N, Int}, j::NTuple{N, Int}) where {T<:Number, N}
  t = zeros(dims)
  ofset = vcat([0], cumsum(counts([j...])))
  fp = fixpointperms(j)
  for i in 1:(prod(dims))
    i = Tuple(CartesianIndices(dims)[i])
    if mapreduce(k -> issorted(i[ofset[k]+1:ofset[k+1]]), *, 1:length(ofset)-1)
      n = rand(T)
      for k in fp
          @inbounds t[i[k]...] = n
      end
    end
  end
  t
end

"""
rand(SymmetricTensor{T, N}, n::Int, b::Int = 2)

Returns N-dimensional random SymmetricTensor with elements of type T drawn from uniform distribution on [0,1),
n denotes data size and b denotes block size.

"""
function rand(::Type{SymmetricTensor{T, N}}, n::Int, b::Int = 2) where {T<:Number, N}
  sizetest(n, b)
  nbar = mod(n,b)==0 ? n÷b : n÷b + 1
  ret = arraynarrays(Float64, fill(nbar, N)...,)
  for j in pyramidindices(N, nbar)
    dims = (mod(n,b) == 0 || !(nbar in j)) ? (fill(b,N)...,) : map(i -> (i == nbar) ? n - b*(nbar-1) : b, j)
    if j == (unique(j)...,)
      @inbounds ret[j...] = rand(T, dims...)
    else
      @inbounds ret[j...] = randblock(T, dims, j)
    end
  end
  SymmetricTensor(ret; testdatstruct = false)
end


text:
module ShaderManager

using CoreExtended
using FileManager
using JLGEngine

using ..GraphicsManager

API = nothing

""" TODO """
function init()
	global API = GRAPHICSDRIVER().ShaderManager
end

""" TODO """
type ShaderProgram <: JLGEngine.IComponent
	id::Symbol
  source::FileSource
	obj::AbstractGraphicsShaderProgram
	ShaderProgram(id=:_) = new(id,FileSource(),API.createShaderProgram())
	ShaderProgram(obj::AbstractGraphicsShaderProgram) = new(:_,FileSource(),obj)
end

""" TODO """
function init(this::ShaderProgram)
end

""" TODO """
function resets()
	API.reset()
end

""" TODO """
function unlinks()
	API.unlink()
end

""" TODO """
function link(this::ShaderProgram)
	API.link(this.obj)
end

""" TODO """
getSource(this::ShaderProgram) = this.source

""" TODO """
setPath(this::ShaderProgram, path::String) = (this.source.path = path)

# preset for manager: includes code here
presetManager(ShaderProgram)

listWatchShaderProgram = []

""" TODO """
listenReload(this::ShaderProgram, f::Function) =
	registerFileUpdate(getSource(this), (source::FileSource,p) -> push!(listWatchShaderProgram, (p,f)), this)

""" TODO """
function watchShaderReload(OnReload=()->nothing)
	global listWatchShaderProgram
	count=length(listWatchShaderProgram)
	if count != 0
		println("Reload Shaders ($count)...")
		for t in listWatchShaderProgram
			p = t[1]
			f = t[2]
			reload(p)
			if isa(f, Function) execute(f,p) end
		end
		listWatchShaderProgram = []
		OnReload()
	end
end

""" TODO """
function create(name::Symbol, path::String, listen::Union{Void,Function}=nothing)
	exists = haskey(list,name)
	this = create(name)
	if !exists
		setPath(this, path)
		if listen != nothing listenReload(this, listen) end
	end
	this
end

""" TODO """
setListenerOnShaderPropertyUpdate(name::Symbol, listen::Tuple) = API.setListenerOnShaderPropertyUpdate(name, listen)

""" TODO """
getShaderProperties(program::ShaderProgram, buffer::AbstractGraphicsData) = API.getShaderProperties(program.obj, buffer)

""" TODO """
render(program::ShaderProgram, buffer::AbstractGraphicsData, data::AbstractGraphicsData) = API.render(program.obj, buffer, data)

#watchShaderReload(f::Function) = API.watchShaderReload(f)
""" TODO """
reload(program::ShaderProgram) = API.reload(program.obj, program.source)

end #ShaderManager

max_stars_repo_name: Gilga/JGE

max_issues_repo_name: Gilga/JGE

max_forks_repo_name: Gilga/JGE


text:
const FRAGMENT_SHADER = 0x00008B30
const VERTEX_SHADER   = 0x00008B31

function Shader(typ::Integer, src::String)
	s = Shader(typ)
	ShaderSource(s, src)
	CompileShader(s)
	return s
end

function Shader(typ::Integer)
	ret = ccall( (:glCreateShader, lib), Shader, (GLenum,), typ)
	if ret == 0
		GetError()
		error("shader creation failed")
	end
	return ret
end

function ShaderSource(s::Shader, src::String)
	ccall( (:glShaderSource, lib), Void,
		(Shader, GLsizei, Ptr{Ptr{GLchar}}, Ptr{GLint}),
		s, 1, [bytestring(src)], [length(src)])
end

function CompileShader(s::Shader)
	ccall( (:glCompileShader, lib), Void, (Shader,), s)
	if !GetShader(s, COMPILE_STATUS)
		error(GetShaderInfoLog(s))
	end
end

const COMPILE_STATUS = 0x8B81

function GetShader(s::Shader, param::Integer)
	ret = GLint[-1]
	ccall( (:glGetShaderiv, lib), Void,
		(Shader, GLenum, Ptr{GLint}), s, param, ret)
	ret = ret[1]
	if ret == -1
		GetError()
		error("no output")
	elseif contains((COMPILE_STATUS, DELETE_STATUS), param)
		return ret == 1
	else
		return ret
	end
end

function GetShaderInfoLog(s::Shader)
	size = GetShader(s, INFO_LOG_LENGTH)
	buf = Array(GLchar, size)
	ccall( (:glGetShaderInfoLog, lib), Void,
		(Shader, GLsizei, Ptr{GLsizei}, Ptr{GLchar}),
		s, size, &size, buf)
	return bytestring(buf)
end


text:
# # Approximate Inference on Sparse GPs using VI

# Loading the necessary packages and setting seed.

using AbstractGPs, Plots, Random
Random.seed!(1234);

# Loading toy regression
# [dataset](https://github.com/GPflow/docs/blob/master/doc/source/notebooks/basics/data/regression_1D.csv)
# taken from GPFlow examples.

x = [0.8658165855998895, 0.6661700880180962, 0.8049218148148531, 0.7714303440386239,
    0.14790478354654835, 0.8666105548197428, 0.007044577166530286, 0.026331737288148638,
    0.17188596617099916, 0.8897812990554013, 0.24323574561119998, 0.028590102134105955];
y = [1.5255314337144372, 3.6434202968230003, 3.010885733911661, 3.774442382979625,
    3.3687639483798324, 1.5506452040608503, 3.790447985799683, 3.8689707574953,
    3.4933565751758713, 1.4284538820635841, 3.8715350915692364, 3.7045949061144983];
scatter(x, y, xlabel="x", ylabel="y")

# Split the observations into train and test set.

(x_train, y_train) = (x[begin:8], y[begin:8]);
(x_test, y_test) = (x[9:end], y[9:end]);

# Instantiate the kernel.

k = Matern52Kernel()

# Instantiate a Gaussian Process with the given kernel `k`.

f = GP(k);

# Instantiate a `FiniteGP`, a finite dimentional projection at the inputs of the dataset
# observed under Gaussian Noise with $\sigma = 0.001$ .

fx = f(x_train, 0.001);

# Calculating the exact posterior over `f` given `y`. The GP's kernel currently has some
# arbitrary fixed parameters.

p_fx = posterior(fx, y_train)

logpdf(p_fx(x_test), y_test)

# ## Variational Inference

# Sanity check for the Evidence Lower BOund (ELBO) implemented according to
# M. K. Titsias's _Variational learning of inducing variables in sparse Gaussian processes_

elbo(fx, y_train, f(rand(7)))

# We will be using `Optim.jl` package's `LBFGS` algorithm to maximize the given ELBO.

using Optim

# Create a helper function for optimization. It takes in the parameters (both variational
# and model parameters) and returns the negative of the ELBO.

function optim_function(params; x=x_train, y=y_train)
    kernel = ScaledKernel(
        transform(
            Matern52Kernel(),
            ScaleTransform(exp.(params[1]))
            ),
        exp.(params[2])
    )
    f = GP(kernel)
    fx = f(x, 0.1)
    return -elbo(fx, y, f(params[3:end]))
end

# Initialize the parameters (Varitational and Model parameters)
x0 = rand(7)

# Sanity check for the helper function. We intend to minimize the result of this function.

-optim_function(x0)

# Optimize using `Optim.jl` package's `LBFGS` algorithm.

opt = optimize(optim_function, x0, LBFGS())

# Optimal parameters:

opt.minimizer

# ELBO with optimal parameters. We see that there is significant improvement when compared
# to the initial parameters.

-optim_function(opt.minimizer; x=x_test, y=y_test)

# Visualize the posterior.

opt_kernel = ScaledKernel(
    transform(
        Matern52Kernel(),
        ScaleTransform(exp.(opt.minimizer[1]))
        ),
    exp.(opt.minimizer[2])
)
opt_f = GP(opt_kernel)
opt_fx = opt_f(x_train, 0.1)
ap = approx_posterior(VFE(), opt_fx, y_train, opt_f(opt.minimizer[3:end]));

# Average log-marginal-probability of data with posterior kernel parameter samples sampled
# using ESS. We can observe that there is significant improvement over exact posterior with
# default kernel parameters.

logpdf(ap(x_test), y_test)

plt = plot(ap, 0:0.001:1, label="Approx Posterior")
plot!(plt, p_fx, 0:0.001:1, label="Exact Posterior")
scatter!(
    plt,
    opt.minimizer[3:end],
    mean(rand(ap(opt.minimizer[3:end], 0.1), 100), dims=2),
    label="Pseudo-points"
)
scatter!(plt, x_train, y_train, label="Train Data")
scatter!(plt, x_test, y_test, label="Test Data")
plt


text:
using Reel
using Test
using Plots
using Compose
import Cairo

@testset "render plots" begin
    function render(t, dt)
        # t is the time into the sequence
        # dt is the time to advance for the next frame

        # any expression that results in an object which can be
        # rendered as png or jpg
        plot([x -> sin(x+t*π), x -> cos(x+t*π)], 0, 6)
    end

    film = roll(render, fps=30, duration=2)

    write("output.gif", film) # Write to a gif file
    write("file.webm", film) # Write to a webm video
    write("file.mp4", film)  # An mp4 formatted video

    @show typeof((render(t,0.1) for t in 0:0.1:2))
    roll((render(t,0.1) for t in 0:0.1:2), fps=30)
    write("output2.gif", film) # Write to a gif file

    ts = []
    film = roll(fps=1, duration=5) do t, dt
        push!(ts, t)
        plot([x -> sin(x+t*π), x -> cos(x+t*π)], 0, 6)
    end
    @test all(diff(ts).==1.0)
end

@testset "do block" begin
    film = roll(fps=30, duration=2) do t, dt
        plot([x -> sin(x+t*π), x -> cos(x+t*π)], 0, 6)
    end

    write("output.gif", film)
end

@testset "Compose" begin
    Compose.set_default_graphic_size(3inch, 3inch) # Square

    # draw a regular n-gon
    ngon(n) = compose(context(units=UnitBox(-1, -1, 2, 2)), fill("lightblue"),
                       polygon([(cos(x), sin(x)) for x in π/2:2π/n:3.5π]))

    film = roll(map(ngon, collect(3:10)), fps=5)

    write("output.gif", film)
end

@testset "Frames" begin
    frames = Frames(MIME("image/png"), fps=2)

    for i=1:31
        t = i * 0.1
        push!(frames, plot([x -> sin(x+t*π), x -> cos(x+t*π)], 0, 6))
    end

    write("output.gif", frames)
    write("output.webm", frames)
end


text:
def print_str(str)
    for i=0, i < |str|, i=i+1
        print(str[i])

def println(var) {
    if var is Object
        print_str(var.to_s())
    else print(var)
	print_str('\n')
}

def print_strln(str) {
	print_str(str)
	print_str('\n')
}

def print_strs(a) print_str(flatten(a))

def print_strsln(a) print_strln(flatten(a))

max_stars_repo_name: jeffreymanzione/jeff-vm-lang

max_issues_repo_name: jeffreymanzione/jeff-vm-lang

max_forks_repo_name: jeffreymanzione/jeff-vm-lang


text:
### A Pluto.jl notebook ###
# v0.17.2

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
        el
    end
end

# ╔═╡ aab79d69-3435-48be-a8c2-752b55867468
using HypertextLiteral

# ╔═╡ d1db2fc1-f465-49e3-b6a7-848137c8948d
md"""
## Basic examples
"""

# ╔═╡ ce83a3d3-ed13-4b0a-80a4-a27f45c49b67

# ╔═╡ 0ebaa1fe-9eea-4c5e-9f3e-9498fefa2869
1:20 ⊆ 0:.2:30

# ╔═╡ 63998497-5fcf-440e-9e2b-d138f0bffbfe

# ╔═╡ 5e6f075e-9719-4ee3-9841-99a153f1f83b
md"""
The range is exact: you can even use rational numbers!
"""

# ╔═╡ 2576d239-2252-48ee-8999-8a3ee2b70a98
md"""
## Alignment
"""

# ╔═╡ eeed3ebd-0431-4fbe-8f1e-09f14df2f87e
import AbstractPlutoDingetjes

# ╔═╡ 2a7d7469-270d-4165-b167-e1dbeb0a9598
const show_compat_warning_ref = Ref{Bool}(false)

# ╔═╡ 8ac624df-306f-4a87-b25a-eb8cd9d1f30b
function closest(range::AbstractRange, x::Real)
	rmin = minimum(range)
	rmax = maximum(range)

	if x <= rmin
		rmin
	elseif x >= rmax
		rmax
	else
		rstep = step(range)

		int_val = (x - rmin) / rstep
		range[round(Int, int_val) + 1]
	end
end

# ╔═╡ 34553132-0392-45a2-8d6d-3c8c2252d01f
begin
	struct RangeSlider
		range::AbstractRange{<:Real}
		show_value::Bool
		default::AbstractRange

		RangeSlider(range::AbstractRange;
				default::AbstractRange=range,
				show_value::Bool=true,
				# for compat
				left=nothing,
				right=nothing,
			) = let

			if show_compat_warning_ref[] && (left !== nothing || right !== nothing)
				@warn "Compat: The keyword arguments `left` and `right` will be removed in the future, use the `default` keyword argument instead."
			end
			start = something(left, minimum(default))
			stop = something(right, maximum(default))

			default = start:step(range):stop

			@assert default ⊆ range "The default range is not a subset of the slidable range. Verify that the beginning and end of the default range are inside the slidable range."
			new(range, show_value, default)
		end
	end

	function Base.show(io, m::MIME"text/html", rs::RangeSlider)
		if !AbstractPlutoDingetjes.is_supported_by_display(io, AbstractPlutoDingetjes.Bonds.transform_value)
			write(io, "❌ Update Pluto to use this element.")
			return
		end
		show(io, m, @htl("""
		<link
			href="https://cdn.jsdelivr.net/npm/nouislider@15.5.0/dist/nouislider.css" rel="stylesheet"/>"""))
	end

	AbstractPlutoDingetjes.Bonds.initial_value(rs::RangeSlider) = rs.default

	function AbstractPlutoDingetjes.Bonds.transform_value(rs::RangeSlider, js_val::Any)
		if js_val isa Vector && length(js_val) == 2

			rounded_range = closest.([rs.range], js_val)
			rounded_range[1] : step(rs.range) : rounded_range[2]
		else
			rs.default
		end
	end

end

# ╔═╡ 2150eb64-4f79-11ec-1eec-a9db61b6d8f3
export RangeSlider

# ╔═╡ 8c2c225b-d142-451e-9a04-296230b73bf3
zbond = @bind z RangeSlider(0.0:π:100)

# ╔═╡ 891e953e-7f6c-4def-8ba0-91f351c5b283
z

# ╔═╡ 3764fff6-3f77-41dd-91fc-228d00474ab5
zbond

# ╔═╡ 48531527-c8ae-46f2-a6b9-d9056ef583bf
@htl("""

<input type=range>
asdf
$(zbond)
<input type=range>
""")

# ╔═╡ ecce7dfb-9885-4fda-9ec0-fd4189b7bc06
z

# ╔═╡ ea0bf76f-e51b-471d-b958-d331b3a64bc8
@bind x RangeSlider(1:100; default=6:20, show_value=false)

# ╔═╡ 8da460e7-c1b6-45b0-aa12-5eab97cdcf24
x

# ╔═╡ 8ec7e04b-3c32-4b38-9aa5-38c14b742393
@bind rational RangeSlider(1:1//3:1000)

# ╔═╡ 1e86ff39-8cc2-4dcf-80a9-b0d45a1e7a0e
rational

# ╔═╡ b9066464-d4df-4a2c-82f4-34eaa50bb487
Float64.(rational)

# ╔═╡ c06b7708-f508-40c1-a34c-c64f17f04c74
md"""
Hello world!

Hello world! $(RangeSlider(1:10; show_value=true)) How are you?

Hello world! $(RangeSlider(1:10; show_value=false)) How are you?

Hello world!

Hello world!
"""

# ╔═╡ 128642a9-5f52-4712-b60b-cc7814a2401b
let
	show_compat_warning_ref[] = false
	result = RangeSlider(1:10; right=6), RangeSlider(1:10; left=6), RangeSlider(1:10; left=3, right=6)
	show_compat_warning_ref[] = true
	result
end

# ╔═╡ ab48cd4a-e4c3-4eca-84a1-689e874ce0e4
@bind closest_test html"<input type=range step=0.000001 max=5>"

# ╔═╡ 4f2f6baa-1deb-4132-bfe5-4a6ef1e8ec86
closest_test

# ╔═╡ 0a4bc67c-67e9-4001-a1c1-c079a984215d
closest(2:.3:4, closest_test isa Missing ? 0 : closest_test)

# ╔═╡ ec579632-c508-4ada-8561-1ed26f1a60cc
RangeSlider(0.0:π:100; show_value=false)

# ╔═╡ 18e7f591-1d3e-42da-84d6-bc5350f85535
RangeSlider(0:5:100; default=50:70)

# ╔═╡ f6a97b41-3513-4dc7-b767-d966c682956a
try
	RangeSlider(0:5:100; default=50:200)
catch e
	sprint(showerror, e) |> Text
end

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
AbstractPlutoDingetjes = "6e696c72-6542-2067-7265-42206c756150"
HypertextLiteral = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"

[compat]
AbstractPlutoDingetjes = "~1.1.1"
HypertextLiteral = "~0.9.3"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

[[AbstractPlutoDingetjes]]
deps = ["Pkg"]
git-tree-sha1 = "0bc60e3006ad95b4bb7497698dd7c6d649b9bc06"
uuid = "6e696c72-6542-2067-7265-42206c756150"
version = "1.1.1"

[[ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[HypertextLiteral]]
git-tree-sha1 = "2b078b5a615c6c0396c77810d92ee8c6f470d238"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.3"

[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[Random]]
deps = ["Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╠═2150eb64-4f79-11ec-1eec-a9db61b6d8f3
# ╟─d1db2fc1-f465-49e3-b6a7-848137c8948d
# ╠═8c2c225b-d142-451e-9a04-296230b73bf3
# ╠═891e953e-7f6c-4def-8ba0-91f351c5b283
# ╠═3764fff6-3f77-41dd-91fc-228d00474ab5
# ╟─ce83a3d3-ed13-4b0a-80a4-a27f45c49b67
# ╠═ea0bf76f-e51b-471d-b958-d331b3a64bc8
# ╠═8da460e7-c1b6-45b0-aa12-5eab97cdcf24
# ╠═0ebaa1fe-9eea-4c5e-9f3e-9498fefa2869
# ╟─63998497-5fcf-440e-9e2b-d138f0bffbfe
# ╟─5e6f075e-9719-4ee3-9841-99a153f1f83b
# ╠═8ec7e04b-3c32-4b38-9aa5-38c14b742393
# ╠═1e86ff39-8cc2-4dcf-80a9-b0d45a1e7a0e
# ╠═b9066464-d4df-4a2c-82f4-34eaa50bb487
# ╟─2576d239-2252-48ee-8999-8a3ee2b70a98
# ╠═48531527-c8ae-46f2-a6b9-d9056ef583bf
# ╠═c06b7708-f508-40c1-a34c-c64f17f04c74
# ╠═ecce7dfb-9885-4fda-9ec0-fd4189b7bc06
# ╠═eeed3ebd-0431-4fbe-8f1e-09f14df2f87e
# ╠═aab79d69-3435-48be-a8c2-752b55867468
# ╠═128642a9-5f52-4712-b60b-cc7814a2401b
# ╠═2a7d7469-270d-4165-b167-e1dbeb0a9598
# ╠═34553132-0392-45a2-8d6d-3c8c2252d01f
# ╠═8ac624df-306f-4a87-b25a-eb8cd9d1f30b
# ╠═ab48cd4a-e4c3-4eca-84a1-689e874ce0e4
# ╠═4f2f6baa-1deb-4132-bfe5-4a6ef1e8ec86
# ╠═0a4bc67c-67e9-4001-a1c1-c079a984215d
# ╠═ec579632-c508-4ada-8561-1ed26f1a60cc
# ╠═18e7f591-1d3e-42da-84d6-bc5350f85535
# ╠═f6a97b41-3513-4dc7-b767-d966c682956a
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
#=##############################################################################################
Filename: rotor_sweep.jl
Author: Ryan Anderson
Contact: rymanderson@gmail.com
README: this is a template file. Convenience methods are provided to prepare a single rotor and
        sweep over advance ratios.
=###############################################################################################

# initialize parameters
"""
    RotorSweepParameters{V1,V2,V3,V4,V5} <: Parameters

# Fields:

* `omegas::V1`
* `Js::V2`
* `Ts::V2`
* `Qs::V2`
* `Ps::V2`
* `CTs::V2`
* `CQs::V2`
* `ηs::V2`
* `us::V3`
* `vs::V3`
* `rotor_names::V4`
* `plot_directory::V5`
* `plot_base_name::V5`
* `plot_extension::V5`

"""
struct RotorSweepParameters{TF} <: Parameters
    omegas::Array{TF,2}
    Js::Array{TF,2}
    Ts::Array{TF,2}
    Qs::Array{TF,2}
    Ps::Array{TF,2}
    CTs::Array{TF,2}
    CQs::Array{TF,2}
    ηs::Array{TF,2}
    us::Vector{Vector{Vector{TF}}}
    vs::Vector{Vector{Vector{TF}}}
    rotor_names::Union{Vector{String}, Vector{LaTeXString}}
    plot_directory::String
    plot_base_name::String
    plot_extension::String
end

# function (rs::RotorSweepParameters{V1, V2, V3, V4, V5})(Js, omega...)

# end

# rs  = RotorSweepParameter(args...)

# rs(Js, omega, nblades...)

# build freestream_function
# function get_J(parameters::RotorSweepParameters{V1,V2,V3,V4,V5,V6}, ri, stepi) where {V1, V2 <:AbstractArray{<:Any,1}, V3, V4, V5, V6}
#     parameters.Js[stepi]
# end

# function get_J(parameters::RotorSweepParameters{V1,V2,V3,V4}, ri, stepi) where {V1, V2 <:AbstractArray{<:Any,2}, V3, V4, V5}
function get_J(parameters, ri, stepi)
    parameters.Js[ri,stepi]
end

"""
    rotor_sweep_template(Js, omegas, rotors;
        rotor_names = ["rotor 1"],
        plot_directory = joinpath(topdirectory, "data", "airfoil", "polars", TODAY),
        plot_base_name = "default",
        plot_extension = ".pdf",
        step_symbol = L"J",
        kwargs...
    )

# Arguments:

* `Js::Array{Float64,2}` : each [i,j]th element contains the desired advance ratio of the ith rotor at the jth step
* `omegas::Array{Float64,2}` : each [i,j]th element contains the commanded rotational velocity in rad/s of the ith rotor at the jth step
* `rotors::CCBladeSystem` : `CCBladeSystem` object defines the rotors.

# Keyword Arguments:

* `Res_list`
* `rotor_names`
* `plot_directory`
* `plot_base_name`
* `plot_extension`
* `step_symbol`

"""
function rotor_sweep_template(Js, omegas, rotors;
    rotor_names = ["rotor 1"],
    plot_directory = joinpath(topdirectory, "data", "airfoil", "polars", TODAY),
    plot_base_name = "default",
    plot_extension = ".pdf",
    step_symbol = L"J",
    kwargs...
)

    # prepare subsystems
    wings = nothing
    inertia = nothing
    nonliftingbodies = nothing
    structures = nothing
    motors = nothing
    batteries = nothing

    # build system struct
    aircraft = Aircraft(wings, rotors, inertia, nonliftingbodies, structures, motors, batteries)

    # compile actions
    actions = [solve_rotor_nondimensional]

    # initialize parameters
    step_range = 1:size(Js)[2] # use time to define each part of the sweep
    # omegas_zero, Js_zero, CTs, CQs, ηs = solve_rotor(aircraft, step_range)
    omegas_zero, Js_zero, Ts, Qs, Ps, CTs, CQs, ηs, us, vs = solve_rotor_nondimensional(aircraft, step_range)

    # check data
    @assert typeof(Js) <: AbstractArray{<:Any,2} "Js must be a 2-dimensional array; got $(typeof(Js))"
    @assert size(omegas_zero) == size(omegas) "`omegas` has improper dimensions. \n\tExpected: $(size(omegas_zero))\n\tGot: $(size(omegas))"
    @assert size(Js_zero) == size(Js) "`Js` has improper dimensions. \n\tExpected: $(size(Js_zero))\n\tGot: $(size(Js))"
    # @warn isdir(plot_directory) "plot_directory does not exist at $plot_directory\n\ttry: `mkdir $plot_directory`"

    if !isdir(plot_directory); mkpath(plot_directory); @warn "plot directory $plot_directory does not exist; creating..."; end

    # build parameters
    parameters = RotorSweepParameters(omegas, Js, Ts, Qs, Ps, CTs, CQs, ηs, us, vs, rotor_names, plot_directory, plot_base_name, plot_extension)

    function freestream_function(aircraft, parameters, environment, step_range, stepi)

        J = get_J(parameters, 1, stepi)
        omega = parameters.omegas[1] # get omega for the first rotor
        n = omega / 2 / pi
        D = aircraft.rotor_system.rotors[1].Rtip * 2

        # calculate freestream
        Vinf = J * n * D
        alpha = 0.0
        beta = 0.0
        Omega = StaticArrays.@SVector zeros(3)
        freestream = Freestream(Vinf, alpha, beta, Omega)

        return freestream
    end

    # build environment_function
    function environment_function(aircraft, parameters, step_range, stepi)
        Environment()
    end

    # assemble postactions
    postactions = [post_plot_rotor_sweep]

    # build objective_function
    objective_function(aircraft, parameters, freestream, environment, step_range) = 0.0

    return aircraft, parameters, actions, freestream_function, environment_function, postactions, objective_function, step_range, step_symbol
end

"""
    rotor_sweep_template(Js, omegas,
        nblades_list, rhub_list, rtip_list, mach_correction,
        radii_list, chords_list, twists_list, ccblade_af_objects_list,
        index, positions, orientations, spin_directions;
            kwargs...
    )

Convenience function to create rotors with pre-created polars.

# Arguments

* `nblades_list::Vector{Int}` : number of blades of the ith rotor
* `rhub_list::Vector{Float64}` : hub radius of the ith rotor
* `rtip_list::Vector{Float64}` : tip radius of the ith rotor
* `mach_correction::Union{Nothing, CCBlade.Mach_correction}` : mach correction method used
* `radii_list::Vector{Vector{Float64}}` : ith element is a vector of radial stations defining the ith rotor
* `chords_list::Vector{Vector{Float64}}` : ith element is a vector of chords defining the ith rotor
* `twists_list::Vector{Vector{Float64}}` : ith element is a vector of twist angles in radians defining the ith rotor
* `ccblade_af_objects_list::Vector{Vector{String}}` : ith element is a vector of airfoil polar objects corresponding to the ith rotor
* `index::Vector{Int64}` : ith element is the index of `ccbladesystem.rotors` (the returned `CCBladeSystem` object) corresponding to the ith rotor of the system
* `positions::Vector{Vector{Float64}}` : ith element is a vector describing the hub position of the ith rotor in the body frame
* `orientations::Vector{Vector{Float64}}` : ith element is a vector describing the hub orientation (positive thrust direction) of the ith rotor in the body frame (back-right-up)
* `spin_directions::Vector{Bool}` : ith element is `true` if the ith rotor rotates in the same direction as the positive thrust direction by the right hand rule

"""
function rotor_sweep_template(Js, omegas, nblades_list, rhub_list, rtip_list, mach_correction, radii_list, chords_list, twists_list, ccblade_af_objects_list, index, positions, orientations, spin_directions; kwargs...)

    # prepare subsystems
    rotors = CCBladeSystem(nblades_list, rhub_list, rtip_list, mach_correction, radii_list, chords_list, twists_list, ccblade_af_objects_list, index, positions, orientations, spin_directions; kwargs...)

    return rotor_sweep_template(Js, omegas, rotors; kwargs...)
end

"""
    rotor_sweep_template(Js, omegas,
        nblades_list, rhub_list, rtip_list, mach_correction,
        radii_list, chords_list, twists_list,

        index, positions, orientations, spin_directions;
            kwargs...
    )

Convenience function to create rotors and polars.

# Arguments

* `Res_lists::Vector{Vector{Vector{Int64}}}` : [i][j]th element is a vector of Reynold's numbers at which the ith rotor, jth section polar is evaluated
* `Ms_lists::Vector{Vector{Vector{Int64}}}` : [i][j]th element is a vector of Mach numbers at which the ith rotor, jth radial section polar is evaluated
* `contour_paths_list::Vector{Vector{String}}` : [i][j]th element is the path to the .dat airfoil contour file corresponding to the jth radial section of the ith rotor

"""
function rotor_sweep_template(Js, omegas, nblades_list, rhub_list, rtip_list, mach_correction, radii_list, chords_list, twists_list, Res_lists, Ms_lists, contour_paths_list, index, positions, orientations, spin_directions;
    polar_directory = joinpath(topdirectory, "data", "airfoil", "polars", TODAY),
    kwargs...
)

    # prepare subsystems
    rotors = CCBladeSystem(nblades_list, rhub_list, rtip_list, mach_correction,
        radii_list, chords_list, twists_list,
        Res_lists, Ms_lists, contour_paths_list,
        index, positions, orientations, spin_directions;
            polar_directory,
            kwargs...
    )

    return rotor_sweep_template(Js, omegas, rotors; kwargs...)
end

max_stars_repo_name: byuflowlab/AircraftSystems

max_issues_repo_name: byuflowlab/AircraftSystems

max_forks_repo_name: byuflowlab/AircraftSystems


text:
function e007(N)
    return get_nth_prime(N)
end

function solution007()
    N = 10001
    answer = e007(N)
    println("The answer to problem 007 is $answer")
end

max_stars_repo_name: mithodin/ttpe

max_issues_repo_name: mithodin/ttpe

max_forks_repo_name: mithodin/ttpe


text:
################################################################################
#
#  Conversion to Nemo/Hecke for storage
#
################################################################################

function toNemo(f::IOStream, a::Dict{T, S}; name::String="R") where {T, S}
  print(f, name, " = Dict{$T,$S}(\n")
  k = collect(keys(a))
  for i=1:length(k)-1
    print(f, k[i], " => ", a[k[i]], ", \n")
  end
  print(f, k[end], " => ", a[k[end]], ");\n")
end

function toNemo(s::String, a::Dict; name::String="R", mode::String ="w")
  f = open(s, mode)
  toNemo(f, a, name = name)
  close(f)
end

function toNemo(f::IOStream, a::Array{T, 1}; name::String="R") where T
  print(f, name, " = [\n")
  for i=1:(length(a)-1)
    print(f, a[i], ",\n")
  end
  print(f, a[end], "];\n")
end

function toNemo(s::String, a::Array{T, 1}; name::String="R", mode::String ="w") where T
  f = open(s, mode)
  toNemo(f, a, name = name)
  close(f)
end


text:
function realtimeop!(Hs::AdiaComputer)
    Hs.state = trotter(-im*(1-Hs.location)*Hs.HB,-im*Hs.location*Hs.HP,3)*Hs.state
    Hs.location += Hs.dt/Hs.maxtime
end

function next_timestep!(Hs::AdiaComputer;evopercentage::Real=1/3,nev=6)
    @assert (0<=Hs.location+evopercentage)&&( Hs.location+evopercentage-1<0.1) "evolutoin percentage out of bounds(should be in [0,1])"

    const evotime = Hs.maxtime
    const dt      = Hs.dt

    # calculate eigen values
    for i=Hs.location*evotime:dt:(Hs.location+evopercentage)*evotime
        realtimeop!(Hs)
        Hs.eigens = [Hs.eigens;eigs(Hamiltonian(Hs),nev=nev,which=:SM)[1].']
    end

    return Hs
end

export next_timestep!

max_stars_repo_name: Roger-luo/QuCmp

max_issues_repo_name: Roger-luo/QuCmp

max_forks_repo_name: Roger-luo/QuCmp


text:
import Base: iterate

export qmr, qmr!

mutable struct LanczosDecomp{T, rT, matT, mat_tT, vecT}
    A::matT
    At::mat_tT

    v_prev::vecT # Orthonormal basis vectors for A
    v_curr::vecT # Orthonormal basis vectors for A
    v_next::vecT # Orthonormal basis vectors for A

    w_prev::vecT # Orthonormal basis vectors for A'
    w_curr::vecT # Orthonormal basis vectors for A'
    w_next::vecT # Orthonormal basis vectors for A'

    α::T
    β_prev::T
    β_curr::T
    δ::T
    resnorm::rT
end
function LanczosDecomp(x, A::matT, b; initially_zero = false) where {matT}
    T = eltype(A)

    # we choose:
    # v_0 = 0, v_1 = b - Ax
    # w_0 = 0, w_1 = v_1
    # v_2 and w_2 are uninitialized and used as workspace

    v_prev = zero(x)
    v_curr = copy(b)
    v_next = similar(x)

    if !initially_zero
        # r0 = A*x - b
        mul!(v_next, A, x)
        axpy!(-one(T), v_next, v_curr)
    end
    resnorm = norm(v_curr)
    rmul!(v_curr, inv(resnorm))

    w_prev = zero(x)
    w_curr = copy(v_curr)
    w_next = similar(x)

    α = zero(T)
    β_prev = zero(T)
    β_curr = zero(T)
    δ = zero(T)

    LanczosDecomp(
        A, adjoint(A),
        v_prev, v_curr, v_next,
        w_prev, w_curr, w_next,
        α, β_prev, β_curr, δ,
        resnorm)
end

start(::LanczosDecomp) = 1
done(l::LanczosDecomp, iteration::Int) = false
function iterate(l::LanczosDecomp, iteration::Int=start(l))
    # Following Algorithm 7.1 of Saad
    if done(l, iteration)
        return nothing
    end

    mul!(l.v_next, l.A, l.v_curr)

    l.α = dot(l.v_next, l.w_curr) # v_next currently contains A*v_curr
    axpy!(-conj(l.α), l.v_curr, l.v_next)
    if iteration > 1 # β_1 = 0
        axpy!(-conj(l.β_curr), l.v_prev, l.v_next)
    end

    mul!(l.w_next, l.At, l.w_curr)
    axpy!(-l.α, l.w_curr, l.w_next)
    if iteration > 1 # δ_1 = 0
        axpy!(-l.δ, l.w_prev, l.w_next)
    end

    vw = dot(l.v_next, l.w_next)
    l.δ = sqrt(abs(vw))
    if iszero(l.δ)
        return nothing
    end

    l.β_prev = l.β_curr
    l.β_curr = vw / l.δ

    rmul!(l.v_next, inv(l.δ))
    rmul!(l.w_next, inv(l.β_curr))

    l.w_next, l.w_curr, l.w_prev = l.w_prev, l.w_next, l.w_curr
    l.v_next, l.v_curr, l.v_prev = l.v_prev, l.v_next, l.v_curr

    return nothing, iteration + 1
end

mutable struct QMRIterable{T, xT, rT, lanczosT}
    x::xT

    lanczos::lanczosT
    resnorm::rT
    tol::rT
    maxiter::Int

    g::Vector{T} # right-hand size following Hessenberg multiplication
    H::Vector{T} # Hessenberg Matrix, computed on CPU

    c_prev::T
    c_curr::T
    s_prev::T
    s_curr::T

    p_prev::xT
    p_curr::xT
end

function qmr_iterable!(x, A, b;
                       abstol::Real = zero(real(eltype(b))),
                       reltol::Real = sqrt(eps(real(eltype(b)))),
                       maxiter::Int = size(A, 2),
                       initially_zero::Bool = false,
                       lookahead::Bool = false)
    T = eltype(x)

    lanczos = LanczosDecomp(x, A, b, initially_zero = initially_zero)

    resnorm = lanczos.resnorm
    g = [resnorm; zero(T)]
    H = zeros(T, 4)

    # Givens rotations
    c_prev, s_prev = one(T), zero(T)
    c_curr, s_curr = one(T), zero(T)

    # Projection operators
    p_prev = zero(x)
    p_curr = zero(x)

    tolerance = max(reltol * lanczos.resnorm, abstol)

    QMRIterable(x,
        lanczos, resnorm, tolerance,
        maxiter,
        g, H,
        c_prev, c_curr, s_prev, s_curr,
        p_prev, p_curr
    )
end

converged(q::QMRIterable) = q.resnorm ≤ q.tol
start(::QMRIterable) = 1
done(q::QMRIterable, iteration::Int) = iteration > q.maxiter || converged(q)

function iterate(q::QMRIterable, iteration::Int=start(q))
    # Check for termination first
    if done(q, iteration)
        return nothing
    end

    iterate(q.lanczos, iteration)

    # for classical Lanczos algorithm, only need 4 elements of last column
    # of Hessenberg matrix
    # H[1] is h_m,(m-2), H[2] is h_m,(m-1), H[3] is h_m,m, and H[4] is h_m,(m+1)
    q.H[2] = conj(q.lanczos.β_prev) # β_m
    q.H[3] = conj(q.lanczos.α)      # α_m
    q.H[4] = q.lanczos.δ            # δ_(m+1)

    # Rotation on H[1] and H[2]. Note that H[1] = 0 initially
    if iteration > 2
        q.H[1] = q.s_prev * q.H[2]
        q.H[2] = q.c_prev * q.H[2]
    end

    # Rotation on H[2] and H[3]
    if iteration > 1
        tmp    = -conj(q.s_curr) * q.H[2] + q.c_curr * q.H[3]
        q.H[2] =       q.c_curr  * q.H[2] + q.s_curr * q.H[3]
        q.H[3] = tmp
    end

    a, b = q.H[3], q.H[4]
    # Compute coefficients for Ω_m, and apply it
    c, s, q.H[3] = givensAlgorithm(q.H[3], q.H[4])

    # Apply Ω_m to rhs
    q.g[2] = -conj(s) * q.g[1]
    q.g[1] =       c  * q.g[1]

    # H is now properly factorized
    # compute p_m
    # we re-use q.lanczos.v_next as workspace for p_m
    q.lanczos.v_next .= q.lanczos.v_prev # we need v_m, not v_m+1
    iteration > 1 && axpy!(-q.H[2], q.p_curr, q.lanczos.v_next)
    iteration > 2 && axpy!(-q.H[1], q.p_prev, q.lanczos.v_next)
    rmul!(q.lanczos.v_next, inv(q.H[3]))

    # iterate x_n = x_n-1 + γ_m p_m
    axpy!(q.g[1], q.lanczos.v_next, q.x)

    # transfer coefficients of Givens rotations for next iteration
    q.c_prev, q.s_prev, q.c_curr, q.s_curr = q.c_curr, q.s_curr, c, s
    q.p_prev .= q.p_curr
    q.p_curr .= q.lanczos.v_next
    q.g[1] = q.g[2]

    # approximate residual
    # See Proposition 7.3 of Saad
    q.resnorm = abs(q.g[2])

    q.resnorm, iteration + 1
end

"""
    qmr(A, b; kwargs...) -> x, [history]

Same as [`qmr!`](@ref), but allocates a solution vector `x` initialized with zeros.
"""
qmr(A, b; kwargs...) = qmr!(zerox(A, b), A, b; initially_zero = true, kwargs...)

"""
    qmr!(x, A, b; kwargs...) -> x, [history]

Solves the problem ``Ax = b`` with the Quasi-Minimal Residual (QMR) method.

# Arguments
- `x`: Initial guess, will be updated in-place;
- `A`: linear operator;
- `b`: right-hand side.

## Keywords

- `initally_zero::Bool`: If `true` assumes that `iszero(x)` so that one
  matrix-vector product can be saved when computing the initial residual
  vector;
- `maxiter::Int = size(A, 2)`: maximum number of iterations;
- `abstol::Real = zero(real(eltype(b)))`,
  `reltol::Real = sqrt(eps(real(eltype(b))))`: absolute and relative
  tolerance for the stopping condition
  `|r_k| ≤ max(reltol * |r_0|, abstol)`, where `r_k = A * x_k - b`
- `log::Bool`: keep track of the residual norm in each iteration;
- `verbose::Bool`: print convergence information during the iteration.

# Return values

**if `log` is `false`**

- `x`: approximate solution.

**if `log` is `true`**

- `x`: approximate solution;

- `history`: convergence history.

[^Saad2003]:
    Saad, Y. (2003). Interactive method for sparse linear system.
[^Freund1990]:
    Freund, W. R., & Nachtigal, N. M. (1990). QMR : for a Quasi-Minimal
    Residual Linear Method Systems. (December).
"""
function qmr!(x, A, b;
              abstol::Real = zero(real(eltype(b))),
              reltol::Real = sqrt(eps(real(eltype(b)))),
              maxiter::Int = size(A, 2),
              lookahead::Bool = false,
              log::Bool = false,
              initially_zero::Bool = false,
              verbose::Bool = false)
    history = ConvergenceHistory(partial = !log)
    history[:abstol] = abstol
    history[:reltol] = reltol
    log && reserve!(history, :resnorm, maxiter)

    iterable = qmr_iterable!(x, A, b; abstol = abstol, reltol = reltol,
                             maxiter = maxiter, initially_zero = initially_zero)

    verbose && @printf("=== qmr ===\n%4s\t%7s\n","iter","resnorm")

    for (iteration, residual) = enumerate(iterable)
        if log
            nextiter!(history)
            push!(history, :resnorm, residual)
        end

        verbose && @printf("%3d\t%1.2e\n", iteration, residual)
    end

    verbose && println()
    log && setconv(history, converged(iterable))
    log && shrink!(history)

    log ? (x, history) : x
end


text:
using OutlierDetection
using OutlierDetectionData
using OutlierDetectionInterface
using OutlierDetectionNetworks
using OutlierDetectionNeighbors
using OutlierDetectionPython

using Documenter
using DocumenterMarkdown

DocMeta.setdocmeta!(OutlierDetection, :DocTestSetup, :(); recursive=true)
makedocs(;
    sitename="OutlierDetection.jl",
    authors="David Muhr <muhrdavid@gmail.com> and contributors",
    repo="https://github.com/OutlierDetectionJL/OutlierDetection.jl/blob/{commit}{path}#L{line}",
    format = Markdown(),
    modules = [
        OutlierDetection,
        OutlierDetectionData,
        OutlierDetectionInterface,
        OutlierDetectionNetworks,
        OutlierDetectionNeighbors,
        OutlierDetectionPython
    ],
    clean=true
)


text:
using SimilaritySearch, TextSearch, InvertedFiles
using Test, SparseArrays, LinearAlgebra, CategoricalArrays, StatsBase, Random
const fit = TextSearch.fit

const text0 = "@user;) #jello.world"
const text1 = "hello world!! @user;) #jello.world :)"
const text2 = "a b c d e f g h i j k l m n o p q"
const corpus = ["hello world :)", "@user;) excellent!!", "#jello world."]

@testset "individual tokenizers" begin
    m = Tokenizer(TextConfig(nlist=[1]))
    @test tokenize(m, text0) == hash.(["@user", ";)", "#jello", ".", "world"])

    m = Tokenizer(TextConfig(nlist=[2]))
    @test decode.(m, tokenize(m, text0)) == ["\tn@user ;)", "\tn;) #jello", "\tn#jello .", "\tn. world"]

    m = Tokenizer(TextConfig(nlist=[3]))
    @test decode.(m, tokenize(m, text0)) == ["\tn@user ;) #jello", "\tn;) #jello .", "\tn#jello . world"]

    m = Tokenizer(TextConfig(qlist=[3]))
    @test decode.(m, tokenize(m, text0)) == map(p -> "\tq" * p, [" @u", "@us", "use", "ser", "er;", "r;)", ";) ", ") #", " #j", "#je", "jel", "ell", "llo", "lo.", "o.w", ".wo", "wor", "orl", "rld", "ld "])
    m = Tokenizer(TextConfig(nlist=[1]))
    @test decode.(m, tokenize(m, text1)) == ["hello", "world", "!!", "@user", ";)", "#jello", ".", "world", ":)"]

    m = Tokenizer(TextConfig(slist=[Skipgram(2,1)]))
    @test decode.(m, tokenize(m, text1)) == map(p -> "\ts" * p, ["hello !!", "world @user", "!! ;)", "@user #jello", ";) .", "#jello world", ". :)"])
end

@testset "Normalize and tokenize" begin
    tok = Tokenizer(TextConfig(del_punc=true, group_usr=true, nlist=[1, 2, 3]))
    @test decode.(tok, tokenize(tok, text1)) == ["hello", "world", "_usr", "#jello", "world", "\tnhello world", "\tnworld _usr", "\tn_usr #jello", "\tn#jello world", "\tnhello world _usr", "\tnworld _usr #jello", "\tn_usr #jello world"]
end

@testset "Normalize and tokenize bigrams and trigrams" begin
    tok = Tokenizer(TextConfig(del_punc=true, group_usr=true, nlist=[2, 3]))
    @test decode.(tok, tokenize(tok, text1)) == ["\tnhello world", "\tnworld _usr", "\tn_usr #jello", "\tn#jello world", "\tnhello world _usr", "\tnworld _usr #jello", "\tn_usr #jello world"]
end

@testset "Tokenize skipgrams" begin
    tok = Tokenizer(TextConfig(del_punc=false, group_usr=false, slist=[Skipgram(3,1)]))
    tokens = tokenize(tok, text1)
    @show text1 tokens
    @test tokens == hash.(["\tshello !! ;)", "\tsworld @user #jello", "\ts!! ;) .", "\ts@user #jello world", "\ts;) . :)"])

    config = Tokenizer(TextConfig(del_punc=false, group_usr=false, nlist=[], slist=[Skipgram(3,1), Skipgram(2, 1)]))
    tokens = tokenize(config, text1)
    @show text1 tokens
    @test tokens == hash.(["\tshello !!", "\tsworld @user", "\ts!! ;)", "\ts@user #jello", "\ts;) .", "\ts#jello world", "\ts. :)", "\tshello !! ;)", "\tsworld @user #jello", "\ts!! ;) .", "\ts@user #jello world", "\ts;) . :)"])
end

@testset "Tokenizer, DVEC, and vectorize" begin
    tok = Tokenizer(TextConfig(group_usr=true, nlist=[1]))
    model = VectorModel(BinaryGlobalWeighting(), FreqWeighting(), compute_bow_corpus(tok, corpus))
    x = vectorize(tok, model, text1)
    @show x
    @show corpus
    @show text1
    @show text2
    @test nnz(x) == 8
    x = vectorize(tok, model, text2)
    @test nnz(x) == 1
end

const sentiment_corpus = ["me gusta", "me encanta", "lo lo odio", "odio esto", "me encanta esto LOL!"]
const sentiment_labels = categorical(["pos", "pos", "neg", "neg", "pos"])
const sentiment_msg = "lol, esto me encanta"

@testset "Tokenizer, DVEC, and vectorize" begin
    config = Tokenizer(TextConfig(group_usr=true, nlist=[1]))
    model = VectorModel(EntropyWeighting(), BinaryLocalWeighting(), compute_bow_corpus(config, sentiment_corpus), sentiment_labels)
    @info model.tokens
    @info sum(t.weight for t in values(model.tokens))
    @test (7.059714 - sum(t.weight for t in values(model.tokens))) < 1e-5
end

@testset "Weighting schemes" begin
    tok = Tokenizer(TextConfig(group_usr=true, nlist=[1]))
    for (gw, lw, dot_) in [
            (BinaryGlobalWeighting(), FreqWeighting(), 0.3162),
            (BinaryGlobalWeighting(), TfWeighting(), 0.3162),
            (BinaryGlobalWeighting(), TpWeighting(), 0.3162),
            (IdfWeighting(), BinaryLocalWeighting(), 0.40518),
            (IdfWeighting(), TfWeighting(), 0.23334),

            (EntropyWeighting(), FreqWeighting(), 0.44456),
            (EntropyWeighting(), TfWeighting(), 0.44456),
            (EntropyWeighting(), TpWeighting(), 0.44456),
            (EntropyWeighting(), BinaryLocalWeighting(), 0.70290)
        ]

        if gw isa EntropyWeighting
            model = VectorModel(gw, lw, compute_bow_corpus(tok, sentiment_corpus), sentiment_labels)
        else
            model = VectorModel(gw, lw, compute_bow_corpus(tok, sentiment_corpus))
        end

        x = vectorize(model, compute_bow(tok, sentiment_corpus[3]))
        y = vectorize(model, compute_bow(tok, sentiment_corpus[4]))
        @show gw, lw, dot_, dot(x, y), x, y
        @test abs(dot(x, y) - dot_) < 1e-3
    end

    for (gw, lw, dot_, p) in [
            (EntropyWeighting(), BinaryLocalWeighting(), 0.70711, 0.9),
            (IdfWeighting(), TfWeighting(), 0.23334, 0.9),
        ]
        if gw isa EntropyWeighting
            model = VectorModel(gw, lw, compute_bow_corpus(tok, sentiment_corpus), sentiment_labels)
        else
            model = VectorModel(gw, lw, compute_bow_corpus(tok, sentiment_corpus))
        end

        model = prune_select_top(model, p)

        x = vectorize(model, compute_bow(tok, sentiment_corpus[3]))
        y = vectorize(model, compute_bow(tok, sentiment_corpus[4]))
        @show gw, lw, dot_, dot(x, y), x, y
        @test abs(dot(x, y) - dot_) < 1e-3
    end
end

@testset "distances" begin
    u = Dict(:el => 0.9, :hola => 0.1, :mundo => 0.2) |> normalize!
    v = Dict(:el => 0.4, :hola => 0.2, :mundo => 0.4) |> normalize!
    w = Dict(:xel => 0.4, :xhola => 0.2, :xmundo => 0.4) |> normalize!

    dist = AngleDistance()
    @test evaluate(dist, u, v) ≈ 0.5975474808029686
    @test evaluate(dist, u, u) <= eps(Float32)
    @test evaluate(dist, w, u) ≈ 1.5707963267948966
end

@testset "operations" begin
    u = Dict(:el => 0.1, :hola => 0.2, :mundo => 0.4)
    v = Dict(:el => 0.2, :hola => 0.4, :mundo => 0.8)
    w = Dict(:el => 0.1^2, :hola => 0.2^2, :mundo => 0.4^2)
    y = Dict(:el => 0.1/9, :hola => 0.2/9, :mundo => 0.4/9)
    @test u == u
    @test u != v
    @test u + u == v
    @test u * u == w
    @test u * (1/9) == y
    @test (1/9) * u == y
    @test dot(normalize!(u + v - v), normalize!(u)) > 0.99
end

_corpus = [
    "la casa roja",
    "la casa verde",
    "la casa azul",
    "la manzana roja",
    "la pera verde esta rica",
    "la manzana verde esta rica",
    "la hoja verde",
]

@testset "transpose bow" begin
    tok = Tokenizer(TextConfig(nlist=[1]))
    model = VectorModel(BinaryGlobalWeighting(), FreqWeighting(), compute_bow_corpus(tok, _corpus))
    X = vectorize_corpus(tok, model, _corpus)
    dX = transpose(X)
end

function are_posting_lists_sorted(invindex)
    for (k, lists) in invindex.lists
        if issorted([p.id for p in lists]) == false
            return false
        end
    end

    true
end

@testset "invindex" begin
    config = Tokenizer(TextConfig(nlist=[1]))
    model = VectorModel(IdfWeighting(), TfWeighting(), compute_bow_corpus(config, _corpus))
    ψ(text) = vectorize(model, compute_bow(config, text))
    invindex = append!(InvertedFile(), vectorize_corpus(config, model, _corpus))
    begin # searching
        q = ψ("la casa roja")
        res = search(invindex, q, KnnResult(4))
        @info res
        @test sort([id for (id, dist) in res]) == [1, 2, 3, 4]
    end
end

@testset "centroid computing" begin
    tok = Tokenizer(TextConfig(nlist=[1]))
    model = VectorModel(BinaryGlobalWeighting(), FreqWeighting(), compute_bow_corpus(tok, _corpus))
    X = vectorize_corpus(tok, model, _corpus)
    x = sum(X) |> normalize!
    vec = bow(tok, x)
    expected = Dict("la" => 0.7366651330405098, "verde" => 0.39921969741172364, "azul" => 0.11248181187626208, "pera" => 0.08712803682959973, "esta" => 0.17425607365919946, "roja" => 0.22496362375252416, "hoja" => 0.11248181187626208, "casa" => 0.33744543562878626, "rica" => 0.17425607365919946, "manzana" => 0.19960984870586182)
    @test 0.999 < dot(vec, expected)
end

@testset "neardup" begin
    function create_corpus()
        alphabet = Char.(97:100)
        corpus = String[]
        L = Char[]
        for i in 1:10000
            resize!(L, 0)
            for j in 1:7
                append!(L, rand(alphabet, 2))
                push!(L, ' ')
            end

            push!(corpus, join(L))
        end

        corpus
    end

    tok = Tokenizer(TextConfig(nlist=[1]))
    randcorpus = create_corpus()
    model = VectorModel(BinaryGlobalWeighting(), TfWeighting(), compute_bow_corpus(tok, randcorpus))
    @show randcorpus[1:10]
    X = vectorize_corpus(tok, model, randcorpus)
    N = neardup(X, 0.2)
    @test length(X) > length(unique(N.nn))
end


text:
@mpi_handle Datatype

const DATATYPE_NULL = _Datatype(MPI_DATATYPE_NULL)
Datatype() = Datatype(DATATYPE_NULL.val)

const MPIInteger = Union{Int8, UInt8, Int16, UInt16, Int32, UInt32, Int64, UInt64}
const MPIFloatingPoint = Union{Float32, Float64}
const MPIComplex = Union{ComplexF32, ComplexF64}

const MPIDatatype = Union{Char,
                            Int8, UInt8, Int16, UInt16, Int32, UInt32, Int64,
                            UInt64,
                            Float32, Float64, ComplexF32, ComplexF64}
MPIBuffertype{T} = Union{Ptr{T}, Array{T}, SubArray{T}, Ref{T}}

MPIBuffertypeOrConst{T} = Union{MPIBuffertype{T}, SentinelPtr}

fieldoffsets(::Type{T}) where {T} = Int[fieldoffset(T, i) for i in 1:length(fieldnames(T))]

# Define a function mpitype(T) that returns the MPI datatype code for
# a given type T. In the case the the type does not exist, it is created and
# then returned. The dictonary is defined in __init__ so the module can be
# precompiled

# accessor and creation function for getting MPI datatypes
function mpitype(::Type{T}) where T
    get!(mpitype_dict, T) do
        if !isbitstype(T)
            throw(ArgumentError("Type must be isbitstype()"))
        end

        # get the data from the type
        fieldtypes = T.types
        offsets = fieldoffsets(T)
        nfields = Cint(length(fieldtypes))

        if nfields == 0  # primitive type
            if sizeof(T) == 0
                error("Can't convert 0-size type to MPI")
            end
            nfields, blocklengths, displacements, types = factorPrimitiveType(T)
        else  # struct
            # put data in MPI format
            blocklengths = ones(Cint, nfields)
            displacements = zeros(Cptrdiff_t, nfields)  # size_t == MPI_Aint ?
            types = Array{MPI_Datatype}(undef, nfields)
            for i=1:nfields
                displacements[i] = offsets[i]
                # create an MPI_Datatype for the current field if it does not exist yet
                types[i] = mpitype(fieldtypes[i])
            end

        end
        newtype = Type_Create_Struct(nfields, blocklengths, displacements, types)

        # commit the datatatype
        Type_Commit!(newtype)
        return newtype.val
    end
end

function factorPrimitiveType(::Type{T}) where {T}

  tsize = sizeof(T)  # size in bytes
  displacements = zeros(Cptrdiff_t, 0)  # size_t == MPI_Aint ?
  types = MPI_Datatype[]

  # put largest sizes first
  mpi_types = [mpitype(Int64), mpitype(Int32), mpitype(Int16), mpitype(Int8)]
  mpi_sizes = [8, 4, 2, 1]
  curr_disp = 0

  while curr_disp != tsize
    remsize = tsize - curr_disp
    for i=1:length(mpi_types)
      size_i = mpi_sizes[i]
      # because each size is a multiple of the smaller sizes, taking the largest
      # size always results in the smallest number of types and doesn't result
      # in avoidable small remainders
      if remsize >= size_i
        push!(types, mpi_types[i])
        push!(displacements, curr_disp)
        curr_disp += size_i
        break
      end
    end
  end

  nfields = length(types)
  blocklengths = ones(Cint, nfields)

  return nfields, blocklengths, displacements, types
end

function Type_Create_Struct(nfields::Integer, blocklengths::MPIBuffertype{Cint},
                            displacements::MPIBuffertype{Cptrdiff_t},
                            types::MPIBuffertype{MPI_Datatype})

    newtype = Datatype()

    # int MPI_Type_create_struct(int count, const int array_of_blocklengths[],
    #                            const MPI_Aint array_of_displacements[],
    #                            const MPI_Datatype array_of_types[], MPI_Datatype *newtype)
    @mpichk ccall((:MPI_Type_create_struct, libmpi), Cint,
                  (Cint, Ptr{Cint}, Ptr{Cptrdiff_t}, Ptr{MPI_Datatype}, Ptr{MPI_Datatype}),
                  nfields, blocklengths, displacements, types, newtype)
    return newtype
end

"""
    Type_Create_Subarray(ndims::Integer, array_of_sizes::MPIBuffertype{Cint},
                         array_of_subsizes::MPIBuffertype{Cint},
                         array_of_starts::MPIBuffertype{Cint}, order::Integer, oldtype)

Creates a data type describing an `ndims`-dimensional subarray of size `array_of_subsizes`
of an `ndims-dimensional` array of size `array_of_sizes` and element type `oldtype`,
starting at the top-left location `array_of_starts`. The parameter `order` refers to
the memory layout of the parent array, and can be either `MPI_ORDER_C` or
`MPI_ORDER_FORTRAN`. Note that, like other MPI data types, the type returned by this
function should be committed with `MPI_Type_commit`.
"""
function Type_Create_Subarray(ndims::Integer,
                              array_of_sizes::MPIBuffertype{Cint},
                              array_of_subsizes::MPIBuffertype{Cint},
                              array_of_starts::MPIBuffertype{Cint},
                              order::Integer,
                              oldtype)

    newtype = Datatype()
    @mpichk ccall((:MPI_Type_create_subarray, libmpi), Cint,
        (Cint, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Cint, MPI_Datatype, Ptr{MPI_Datatype}),
        ndims, array_of_sizes, array_of_subsizes, array_of_starts,
        order, mpitype(oldtype), newtype)
    return newtype
end

function Type_Contiguous(count::Integer, oldtype)
    newtype = Datatype()
    @mpichk ccall((:MPI_Type_contiguous, libmpi), Cint,
                  (Cint, MPI_Datatype, Ptr{MPI_Datatype}),
                  count, oldtype, newtype)
    return newtype
end

function Type_Commit!(newtype::Datatype)
    # int MPI_Type_commit(MPI_Datatype *datatype)
    @mpichk ccall((:MPI_Type_commit, libmpi), Cint,
                  (Ptr{MPI_Datatype},), newtype)
end

"""
  Setter function for mpitype_dict and mpitype_dict_inverse
"""
function recordDataType(T::DataType, mpiT::MPI_Datatype)

  if !haskey(mpitype_dict, T)
      mpitype_dict[T] = mpiT
  end

  if !haskey(mpitype_dict_inverse, mpiT)
      mpitype_dict_inverse[mpiT] = T
  end

  return nothing
end
recordDataType(T::DataType, dtyp::Datatype) = recordDataType(T, dtyp.val)

const mpitype_dict = Dict{DataType, MPI_Datatype}()
const mpitype_dict_inverse = Dict{MPI_Datatype, DataType}()

function init_datatypes()
    if Sys.iswindows() || MPI_VERSION >= v"2.2"
        # use specific-width types if available
        for (T, c) in [
            Int8 => MPI_INT8_T,
            UInt8 => MPI_UINT8_T,
            Int16 => MPI_INT16_T,
            UInt16 => MPI_UINT16_T,
            Int32 => MPI_INT32_T,
            UInt32 => MPI_UINT32_T,
            Int64 => MPI_INT64_T,
            UInt64 => MPI_UINT64_T,
            ComplexF32 => MPI_C_FLOAT_COMPLEX,
            ComplexF64 => MPI_C_DOUBLE_COMPLEX,
        ]
            recordDataType(T, MPI_Datatype(c))
        end
    end

    for (T, c) in [
        UInt8 => MPI_BYTE,
        Cshort => MPI_SHORT,
        Cushort => MPI_UNSIGNED_SHORT,
        Cint => MPI_INT,
        Cuint => MPI_UNSIGNED,
        Clong => MPI_LONG,
        Culong => MPI_UNSIGNED_LONG,
        Cchar => MPI_CHAR,
        Cchar => MPI_SIGNED_CHAR,
        Cuchar => MPI_UNSIGNED_CHAR,
        Cwchar_t => MPI_WCHAR,
        Float32 => MPI_FLOAT,
        Float64 => MPI_DOUBLE,

        Char => MPI_UNSIGNED,
    ]
        recordDataType(T, MPI_Datatype(c))
    end
end

push!(mpi_init_hooks, init_datatypes)

function Get_address(location::MPIBuffertype{T}) where T
    addr = Ref{Cptrdiff_t}(0)
    @mpichk ccall((:MPI_Get_address, libmpi), Cint, (Ptr{T}, Ref{Cptrdiff_t}), location, addr)
    return addr[]
end


text:
using KernelFunctions

k = SqExponentialKernel()

@testset "Generic functions to test" begin
    @test length(k) == 1
    @test iterate(k) == (k,nothing)
    @test iterate(k,1) == nothing
end


text:
module EMLmultiply
	using DanaTypes
	using DotPlusInheritance
	using Reexport
	@reexport using ...types.EMLtypes
	import EMLtypes.length
	include("multiply/Multiply.jl")
end


text:
module CompositeLayupModule

using FinEtools
using LinearAlgebra: norm, Transpose, mul!, I
using FinEtoolsDeforLinear.MatDeforLinearElasticModule: tangentmoduli!
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures.TransformerModule: QEQTTransformer

"""
    cartesian_csys(axes)

Create a material Cartesian coordinate system.

- `axes` = tuple of signed labels of the axes. For instance, `(1, 2, 3)` creates
  a coordinate system identical to the global cartesian coordinate system. `
  (2, 1, -3)` creates a coordinate system so that the first material basis
  vector is along the second global basis vector, the second material basis
  vector is along the first global basis vector, and the third material basis
  vector is opposite to the third global basis vector.
"""
function cartesian_csys(axes)
    function   cartesian!(csmatout::FFltMat, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
        csmatout[:] .= 0.0
        for j in 1:3
            aj = abs(axes[j]); sj = sign(axes[j])
            csmatout[aj, j] = sj * 1.0
        end
        return csmatout
    end
    return  CSys(3, 3, cartesian!)
end

abstract type AbstractPly end

"""
    Ply{M} <: AbstractPly

Type of a ply.
"""
struct Ply{M} <: AbstractPly
    name::String # name of the ply
    material::M # material of the ply
    thickness::FFlt # thickness of the ply
    angle::FFlt # angle in degrees
    _Dps::FFltMat # plane stress stiffness matrix
    _Dts::FFltMat # transverse shear stiffness matrix
end

"""
    Ply{M} <: AbstractPly

Create a ply.

Provide name, material of the ply, thickness of the ply, an angle between the
first bases vector of the layup coordinate system and the first direction of
the ply material coordinate system.
"""
function Ply(name, material::M, thickness, angle) where {M}
    # First we extract the full three dimensional stiffness matrix of the
    # material
    D = fill(0.0, 6, 6)
    # What if the material changes composition from point to point? What do we
    # do now?
    t::FFlt, dt::FFlt, loc::FFltMat, label::FInt = 0.0, 0.0, [0.0 0.0 0.0], 0
    tangentmoduli!(material,  D,  t, dt, loc, label)
    # Next we reduce the 3d stiffness to plane stress
    # Both stiffness matrices are expressed on the reference coordinate system
    # for the ply: direction 1 = fiber direction, direction 3 = normal to the
    # ply
    Dps = fill(0.0, 3, 3)
    Dps[1:2, 1:2] = D[1:2, 1:2] -  (reshape(D[1:2,3], 2, 1) * reshape(D[3,1:2], 1, 2))/D[3, 3]
    ix = [1, 2, 4];
    for i = 1:3
        Dps[3,i] = Dps[i,3] = D[4, ix[i]];
    end
    # And we also extract the transverse shear matrix
    Dts = fill(0.0, 2, 2)
    ix = [5, 6];
    for i = 1:2
        Dts[i,i] = D[ix[i], ix[i]];
    end
    return Ply(name, material, thickness, FFlt(angle), Dps, Dts)
end

"""
    lamina_material(E1, E2, nu12, G12, G13, G23)

Create a transversely isotropic lamina material with default (zero) mass
density.
"""
function lamina_material(E1, E2, nu12, G12, G13, G23)
    rho = 0.0
    return MatDeforElastOrtho(DeforModelRed3D, rho, E1, E2, E2, nu12, nu12, 0.0, G12, G13, G23, 0.0, 0.0, 0.0)
end

"""
    lamina_material(rho, E1, E2, nu12, G12, G13, G23)

Create a transversely isotropic lamina material.
"""
function lamina_material(rho, E1, E2, nu12, G12, G13, G23)
    return MatDeforElastOrtho(DeforModelRed3D, rho, E1, E2, E2, nu12, nu12, 0.0, G12, G13, G23, 0.0, 0.0, 0.0)
end

"""
    lamina_material(E, nu)

Create an isotropic lamina material with default (zero) mass
density.
"""
function lamina_material(E, nu)
    rho = 0.0
    return MatDeforElastIso(DeforModelRed3D, rho, E, nu, 0.0)
end

"""
    lamina_material(rho, E, nu)

Create an isotropic lamina material.
"""
function lamina_material(rho, E, nu)
    return MatDeforElastIso(DeforModelRed3D, rho, E, nu, 0.0)
end

"""
    CompositeLayup

Type for composite layup.
"""
struct CompositeLayup
    name::String
    offset::FFlt # offset of the reference surface from the mid surface, negative when the offset is against the normal
    plies::Vector{AbstractPly} # vector of plies; the first ply is at the bottom of the shell (SNEG surface), the last ply is at the top (SPOS)
    csys::CSys # updater of the material orientation matrix
end

"""
    CompositeLayup(name, plies, mcsys)

Create a composite layup.

Provide the name, the array of plies, and the coordinate system that defines the
orientation of the composite layup. The first base spector of this coordinate
system is the reference direction for the layup.
"""
function CompositeLayup(name, plies, mcsys)
    offset = 0.0
    return CompositeLayup(name, offset, plies, mcsys)
end

"""
    thickness(cl::CompositeLayup)

Compute the thickness of the layup (some of the thicknesses of the plies).
"""
function thickness(cl::CompositeLayup)
    return sum(p.thickness for p in cl.plies)
end

"""
    laminate_stiffnesses!(cl::CompositeLayup, A, B, C)

Compute the laminate stiffnesses, membrane, extension-bending coupling, and bending.
"""
function laminate_stiffnesses!(cl::CompositeLayup, A, B, C)
    # Aij coefficients represent in-plane stiffness of the laminate, the Cij
    # coefficients represent bending stiffness, the Bij represent
    # bending-extension coupling, and the Hij represent intralaminar shear
    # stiffness.
    A .= zero(eltype(A))
    B .= zero(eltype(A))
    C .= zero(eltype(A))
    Dps = deepcopy(A)
    T = deepcopy(A)
    tf = QEQTTransformer(Dps)
    # Transform into the composite layup coordinate system.
    layup_thickness = thickness(cl)
    zs = -layup_thickness/2 - cl.offset
    for p in cl.plies
        ze = zs + p.thickness
        plane_stress_T_matrix!(T, -p.angle/180*pi)
        # Transform the plane stress matrix into the layup coordinates
        @. Dps = p._Dps
        Dps = tf(Dps, T)
        # Compute the in-plane stiffness
        @. A += (ze - zs) * Dps
        # Compute the extension-bending coupling stiffness
        @. B += (ze^2 - zs^2)/2 * Dps
        # Compute the bending stiffness
        @. C += (ze^3 - zs^3)/3 * Dps
        zs += p.thickness
    end
    return A, B, C
end

"""
    laminate_transverse_stiffness!(cl::CompositeLayup, H)

Computed the laminate transverse stiffness.
"""
function laminate_transverse_stiffness!(cl::CompositeLayup, H)
    # Hij represent intralaminar shear stiffness.
    H .= zero(eltype(H))
    Dts = deepcopy(H)
    T = deepcopy(H)
    tf = QEQTTransformer(Dts)
    # Transform into the composite layup coordinate system.
    layup_thickness = sum(p.thickness for p in cl.plies)
    zs = -layup_thickness/2 - cl.offset
    for p in cl.plies
        ze = zs + p.thickness
        transverse_shear_T_matrix!(T, p.angle/180*pi)
        # Transform the plane stress matrix into the layup coordinates
        @. Dts = p._Dts
        Dts = tf(Dts, T)
        # Compute the transverse shear stiffness The correction factor
        # accounting for the parable distribution of the shear stress is due to
        # Vinson, Sierakowski, The Behavior of Structures Composed of Composite
        # Materials, 2008
        @. H += 5/4*(ze - zs - 4/3*(ze^3-zs^3)/layup_thickness^2) * Dts
        zs += p.thickness
    end
    return H
end

function laminate_inertia!(cl::CompositeLayup)
    layup_thickness = sum(p.thickness for p in cl.plies)
    zs = -layup_thickness/2 - cl.offset
    mass_density = 0.0
    moment_of_inertia_density = 0.0
    for p in cl.plies
        rho = massdensity(p.material); # mass density
        ze = zs + p.thickness
        mass_density += (ze - zs) * rho
        moment_of_inertia_density += (ze^3 - zs^3) * rho / 3
        zs += p.thickness
    end
    return mass_density, moment_of_inertia_density
end

"""
    plane_stress_T_matrix!( Tinvm::Array{T, 2},
     angle) where {T}

Compute the transformation matrix between strain tensor components on the layup
coordinate system into the ply coordinate system.

`angle` = angle between the first basis vector of the layup coordinate system
    and the first basis vector of the ply coordinate system
"""
function  plane_stress_T_matrix!(Tm::Array{T, 2}, angle) where {T}
    m=cos(angle);
    n=sin(angle);
    return plane_stress_T_matrix!(Tm, m, n)
end

function  plane_stress_T_matrix!(Tm::Array{T, 2}, m, n) where {T}
    Tm[1, 1] =  (m^2);  Tm[1, 2] = (n^2);   Tm[1, 3] = (2*m*n)
    Tm[2, 1] = (n^2);   Tm[2, 2] = (m^2);   Tm[2, 3] = (-2*m*n)
    Tm[3, 1] = (-m*n);  Tm[3, 2] = (m*n);   Tm[3, 3] = (m*m-n*n)
    return Tm
end

"""
    plane_stress_Tinv_matrix!(Tinvm::Array{T, 2}, angle) where {T}

Compute the transformation matrix between strain tensor components on the ply
coordinate system into the layout coordinate system.

`angle` = angle between the first basis vector of the layup coordinate system
    and the first basis vector of the ply coordinate system
"""
function  plane_stress_Tinv_matrix!(Tinvm::Array{T, 2}, angle) where {T}
    m=cos(angle);
    n=sin(angle);
    Tinvm[1, 1] = (m^2); Tinvm[1, 2] = (n^2);  Tinvm[1, 3] = (-2*m*n)
    Tinvm[2, 1] = (n^2); Tinvm[2, 2] = (m^2);  Tinvm[2, 3] = (2*m*n)
    Tinvm[3, 1] = (m*n); Tinvm[3, 2] = (-m*n); Tinvm[3, 3] = (m*m-n*n)
    return Tinvm
end

"""
    plane_stress_T_matrix_eng!( Tinvm::Array{T, 2},
     angle) where {T}

Compute the transformation matrix between strain engineering components on the
layup coordinate system into the ply coordinate system.

`angle` = angle between the first basis vector of the layup coordinate system
    and the first basis vector of the ply coordinate system

The matrix of transformation is `Tme = R * Tm / R`, where `Tm`
is the transformation matrix in tensor components, and `R` is the Reuter matrix,
`R = [1 0 0; 0 1 0; 0 0 2]`.
"""
function  plane_stress_T_matrix_eng!(Tm::Array{T, 2}, angle) where {T}
    m=cos(angle);
    n=sin(angle);
    Tm[1, 1] =  (m^2);  Tm[1, 2] = (n^2);   Tm[1, 3] = (m*n)
    Tm[2, 1] = (n^2);   Tm[2, 2] = (m^2);   Tm[2, 3] = (-m*n)
    Tm[3, 1] = (-2*m*n);  Tm[3, 2] = (2*m*n);   Tm[3, 3] = (m*m-n*n)
    return Tm
end

"""
    plane_stress_Tinv_matrix!(Tinvm::Array{T, 2}, angle) where {T}

Compute the transformation matrix between strain engineering components on the
ply coordinate system into the layout coordinate system.

`angle` = angle between the first basis vector of the layup coordinate system
    and the first basis vector of the ply coordinate system

The matrix of transformation is `Tme = R * Tm / R`, where `Tm`
is the transformation matrix in tensor components, and `R` is the Reuter matrix,
`R = [1 0 0; 0 1 0; 0 0 2]`.
"""
function  plane_stress_Tinv_matrix_eng!(Tinvm::Array{T, 2}, angle) where {T}
    m=cos(angle);
    n=sin(angle);
    Tinvm[1, 1] = (m^2); Tinvm[1, 2] = (n^2);  Tinvm[1, 3] = (-m*n)
    Tinvm[2, 1] = (n^2); Tinvm[2, 2] = (m^2);  Tinvm[2, 3] = (m*n)
    Tinvm[3, 1] = (2*m*n); Tinvm[3, 2] = (-2*m*n); Tinvm[3, 3] = (m*m-n*n)
    return Tinvm
end

"""
    transverse_shear_T_matrix!(Tm::Array{T, 2}, angle) where {T}

Compute the transformation matrix for the transverse shear stiffness.
"""
function  transverse_shear_T_matrix!(Tm::Array{T, 2}, angle) where {T}
    m=cos(angle);
    n=sin(angle);
    return transverse_shear_T_matrix!(Tm, m, n)
end

"""
    transverse_shear_T_matrix!(Tm::Array{T, 2}, m, n) where {T}

Compute the transformation matrix for the transverse shear stiffness.
"""
function  transverse_shear_T_matrix!(Tm::Array{T, 2}, m, n) where {T}
    # Barbero, Introduction, a in equation 5.7
    Tm[1, 1] =  m;  Tm[1, 2] = n;
    Tm[2, 1] = -n;  Tm[2, 2] = m;
    return Tm
end

end # module


text:
function get_k(obj::RiccatiEquationSolver)
    return jcall(obj, "getK", RealMatrix, ())
end

function get_p(obj::RiccatiEquationSolver)
    return jcall(obj, "getP", RealMatrix, ())
end


text:
struct EMR_MSM_Prediction{T<:AbstractFloat, S<:MSM_Timeseries{T}}
    est::EMR_MSM_Estimate
    pred_timeseries::AbstractArray{S,1}
    start_ind::Array{Integer,1}
end

# TODO: enlarge existing prediction
function EMR_MSM_Prediction(
    est::EMR_MSM_Estimate{T, M, S},
    timesteps::Array{A,1},
    num_samples::Integer = 1,
    start_ind::Array{I,1} = length.(est.timeseries),
    last_layer_residuals::Union{Nothing, Array{B,1}} = nothing;
    rand_last_layer_res = true,
    use_last_layer = true
    ) where {T<:AbstractFloat, M<:EMR_MSM_Model_Estimate{T},
            S<:MSM_Timeseries{T}, A<:AbstractArray{T,1},
            B<:AbstractArray{T,1}, I<:Integer}

    pred_timeseries = [MSM_Timeseries(est.pred_timeseries[ts], timesteps[ts],
                        num_samples)
                        for ts in 1:length(est.pred_timeseries)]
    EMR_MSM_Prediction(pred_timeseries, est, start_ind, last_layer_residuals;
        rand_last_layer_res = rand_last_layer_res,
        use_last_layer = use_last_layer)
end

function EMR_MSM_Prediction(
    est::EMR_MSM_Estimate{T, M, S},
    num_preds::Integer,
    timestep::Array{T,1} = [one(T) for i in 1:length(est.timeseries)],
    num_samples::Integer = 1,
    start_ind::Array{I,1} = length.(est.timeseries),
    last_layer_residuals::Union{Nothing, Array{B,1}} = nothing;
    rand_last_layer_res = true,
    use_last_layer = true
    ) where {T<:AbstractFloat, M<:EMR_MSM_Model_Estimate{T},
            B<:AbstractArray{T,1}, S<:MSM_Timeseries{T}, I<:Integer}

    pred_timeseries = [MSM_Timeseries(est.pred_timeseries[ts], num_preds, timestep[ts],
        num_samples)
        for ts in 1:length(est.pred_timeseries)]

    EMR_MSM_Prediction(pred_timeseries, est, start_ind, last_layer_residuals;
        rand_last_layer_res = rand_last_layer_res,
        use_last_layer = use_last_layer)
end

# TODO use est.pred_ts not est.timeseries !!

function EMR_MSM_Prediction(pred_timeseries::Array{MSM_Timeseries_Point{T},1},
        est::EMR_MSM_Estimate{T, EMR_MSM_Model_PointEstimate{T}, S,
            MSM_Timeseries_Dist{T}},
        start_ind::Array{I,1} = length.(est.timeseries),
        last_layer_residuals::Union{Nothing, Array{B,1}} = nothing;
        rand_last_layer_res = true,
        use_last_layer = true
        ) where {T<:AbstractFloat, S<:MSM_Timeseries{T},
            B<:AbstractArray{T,1}, I<:Integer}

    EMR_MSM_Prediction(pred_timeseries,
        EMR_MSM_Estimate(est.model,est.timeseries,
        [rand(ps.ts_points) for ps in est.pred_timeseries]),
        start_ind,
        last_layer_residuals;
        rand_last_layer_res = rand_last_layer_res,
        use_last_layer = use_last_layer
        )
end

function EMR_MSM_Prediction(pred_timeseries::Array{MSM_Timeseries_Point{T},1},
        est::EMR_MSM_Estimate{T, EMR_MSM_Model_PointEstimate{T}, S,
            MSM_Timeseries_Point{T}},
        start_ind::Array{I,1} = length.(est.timeseries),
        last_layer_residuals::Union{Nothing, Array{B,1}} = nothing;
        rand_last_layer_res = true,
        use_last_layer = true
        ) where {T<:AbstractFloat, S<:MSM_Timeseries{T},
            B<:AbstractArray{T,1}, I<:Integer}

    num_layers = layers(pred_timeseries[1])
    num_params = params(pred_timeseries[1])

    L_mats = est.model.RCorrs
    A_mat = est.model.Lin_mat
    B_mats = est.model.Quad_mats
    F = est.model.F
    μ = est.model.μ
    Σ = est.model.Σ

    for ts in 1:length(pred_timeseries)
        r = residuals(pred_timeseries[ts])
        x = values(pred_timeseries[ts])

        x[1,:] = values(est.timeseries[ts])[start_ind[ts], :]
        r[1,:,:] = residuals(est.pred_timeseries[ts])[start_ind[ts], :, :]
        pred_timesteps = timesteps(pred_timeseries[ts])

        for i=2:length(pred_timeseries[ts])
            r[i,:,num_layers] = if isnothing(last_layer_residuals)
                ifelse(use_last_layer,rand(MvNormal(μ, Σ)), zeros(T,num_params))
            elseif rand_last_layer_res
                last_layer_residuals[ts][rand(1:size(rand_last_layer_res,1)),:]
            else
                last_layer_residuals[ts][i,:]
            end

            for l = (num_layers-1):-1:1
                r[i,:,l] = pred_timesteps[i]*(L_mats[l]*vcat(x[i-1,:],vec(r[i-1,:,1:l]))) +
                    r[i-1,:,l+1] + r[i-1,:,l]
            end

            x[i,:] = pred_timesteps[i]*(-A_mat*x[i-1,:] +
                [transpose(x[i-1,:])*B_mats[:,:,k]*x[i-1,:] for k in 1:num_params] +
                F) + r[i-1,:,1] + x[i-1,:]
        end

    end
    EMR_MSM_Prediction{T, MSM_Timeseries_Point{T}}(est, pred_timeseries, start_ind)
end

function EMR_MSM_Prediction(pred_timeseries::Array{MSM_Timeseries_Dist{T},1},
        est::EMR_MSM_Estimate{T, EMR_MSM_Model_DistEstimate{T}, S},
        start_ind::Array{I,1} = length.(est.timeseries),
        last_layer_residuals::Union{Nothing, Array{B,1}} = nothing;
        rand_last_layer_res = true,
        use_last_layer = true
        ) where {T<:AbstractFloat, S<:MSM_Timeseries{T}, B<:AbstractArray{T,1},  I<:Integer}

    for ts in 1:length(pred_timeseries)
        for which_ts in 1:samples(pred_timeseries[ts])
            r_ind = rand(1:size(est.model))
            EMR_MSM_Prediction(
                [MSM_Timeseries_Point{T}(
                    view(pred_timeseries[ts].x,:,:,which_ts),
                    view(pred_timeseries[ts].residuals,:,:,:,which_ts),
                    view(pred_timeseries[ts].timesteps,:))],
                EMR_MSM_Estimate(
                    est.model.estimates[r_ind],
                    [est.timeseries[ts]],
                    [est.pred_timeseries[ts].ts_points[r_ind]]
                    ),
                [start_ind[ts]],
                if !isnothing(last_layer_residuals)
                    if typeof(last_layer_residuals[ts]) <: AbstractArray{T,3}
                        [last_layer_residuals[ts][:,:,which_ts]]
                    elseif typeof(last_layer_residuals[ts]) <: AbstractArray{T,2}
                        [last_layer_residuals[ts]]
                    else
                        nothing
                    end
                else
                    nothing
                end ; rand_last_layer_res = rand_last_layer_res,
                use_last_layer = use_last_layer)
        end
    end

    EMR_MSM_Prediction{T, MSM_Timeseries_Dist{T}}(
        est, pred_timeseries, start_ind)
end

function write(output::String, pred::EMR_MSM_Prediction)
    mkpath(output)
    for i in 1:length(pred.pred_timeseries)
        write(joinpath(output,"ts-"*string(i)), pred.pred_timeseries[i])
    end
end


text:
using Base.Test
using VegaLite

@test_throws ErrorException config(abcd="trus")
@test_throws ErrorException config(viewport="trus")
@test isa(config(viewport=[400,330]), VegaLiteVis)

xs = collect(linspace(0,2,100))
ys = Float64[ rand()*0.1 + cos(x) for x in xs]

@test isa(data_values(time=xs, res=ys), VegaLiteVis)
@test isa(mark_line(), VegaLiteVis)
@test isa(encoding_x_quant(:time), VegaLiteVis)
@test isa(encoding_y_quant(:res), VegaLiteVis)

using RDatasets
mpg = dataset("ggplot2", "mpg")

@test isa(data_values(mpg), VegaLiteVis)
@test isa(mark_point(), VegaLiteVis)

@test isa(encoding_x_quant(:Cty, axis=false), VegaLiteVis)
@test isa(encoding_y_quant(:Hwy, scale=scale(zero=false)), VegaLiteVis)
@test isa(encoding_color_nominal(:Manufacturer), VegaLiteVis)
@test isa(config_cell(width=350, height=400), VegaLiteVis)

@test isa(encoding_x_ord(:Year,
             axis  = axis(labelAngle=-45, labelAlign="right"),
             scale = scale(bandSize=50)), VegaLiteVis)

@test isa(encoding_column_ord(:Cyl), VegaLiteVis)
@test isa(encoding_row_ord(:Year), VegaLiteVis)

@test isa(encoding_x_quant(:Displ), VegaLiteVis)
@test isa(encoding_y_quant(:Hwy), VegaLiteVis)

@test isa(encoding_size_quant(:Cty), VegaLiteVis)

@test isa(mark_text(), VegaLiteVis)

@test isa(encoding_text_quant(:Displ, aggregate="mean"), VegaLiteVis)
@test isa(config_mark(fontStyle="italic", fontSize=12, font="courier"), VegaLiteVis)

println("Finished")


text:
"""
aa_qr_update(Q, R, vnew, m, k, Qd)

Update the QR factorization for the Anderson acceleration optimization
problem.

Still need to make the allocation for Qd go away.
"""
function aa_qr_update!(Q, R, vnew, m, k, Qd)
(n,m)=size(Q)
    aaqr_dim_check(Q, R, vnew, m, k)
    if k == 0
        R[1, 1] = norm(vnew)
        @views Q[:, 1] .= vnew / norm(vnew)
    else
        if k > m - 1
            downdate_aaqr!(Q, R, m, Qd)
        end # inner if block
        kq = min(k, m - 1)
        update_aaqr!(Q, R, vnew, m, kq)
    end # outer if block
    return (Q, R)
end

function update_aaqr!(Q, R, vnew, m, k)
    (nq, mq) = size(Q)
    (k > m - 1) && error("Dimension error in Anderson QR")
    @views Qkm=Q[:,1:k]
    @views hv = vec(R[1:k+1, k+1])
    Orthogonalize!(Qkm, hv, vnew, "cgs2")
    @views R[1:k+1, k+1] .= hv
    @views Q[:, k+1] .= vnew
#    return (Q = Q, R = R)
end

function downdate_aaqr!(Q, R, m, Qd)
(nq,mq)=size(Q)
(pd,md)=size(Qd)
(md == m-1) || @error("dimension error in downdate")
    @views Rp = R[:, 2:m]
    G = qr!(Rp)
    Rd = Matrix(G.R)
    Qx = Matrix(G.Q)
    @views R[1:m-1, 1:m-1] .= Rd
    @views R[:,m].=0.0
if (pd==nq)
    mul!(Qd,Q,Qx)
    @views Q[:,1:m-1] .= Qd
else
    blocksize=pd
    (dlow,dhigh)=blockdim(nq,blocksize)
    blen=length(dlow)
    for il=1:blen
    asize=dhigh[il]-dlow[il]+1
    @views QZ=Qd[1:asize,:]
    @views Qsec=Q[dlow[il]:dhigh[il],:]
    @views mul!(QZ,Qsec,Qx)
    @views Qsec[:,1:m-1] .= QZ
    end
end
    @views Q[:, m] .= 0.0
    return (Q, R)
end

function aaqr_dim_check(Q, R, vnew, m, k)
    (mq, nq) = size(Q)
    (mr, nr) = size(R)
    n = length(vnew)
    dimqok = ((mq == n) && (nq == m))
    dimrok = ((mr == m) && (nr == m))
    dimok = (dimqok && dimrok)
    dimok || error("array size error in AA update")
end

function blockdim(n, block)
    p = Int(floor(n / block))
    res = n - p * block
    ilow = Int64[]
    ihigh = Int64[]
    for jb = 1:p
        lowval = (jb - 1) * block + 1
        push!(ilow, lowval)
        highval = ilow[jb] + block - 1
        push!(ihigh, highval)
    end
    if res > 0
        lowval = p * block + 1
        push!(ilow, lowval)
        push!(ihigh, n)
    end
    return (ilow, ihigh)
end