text:
function Base.:<(N::T, n::NT) where {T <: QuantitativeNetwork, NT <: Number}
   ReturnType = T <: AbstractBipartiteNetwork ? BipartiteNetwork : UnipartiteNetwork
   return ReturnType(N.edges .< n, _species_objects(N)...)
end

function Base.:<(N::T, n::NT) where {T <: ProbabilisticNetwork, NT <: AbstractFloat}
   ReturnType = T <: AbstractBipartiteNetwork ? BipartiteNetwork : UnipartiteNetwork
   return ReturnType(N.edges .< n, _species_objects(N)...)
end

function Base.:<(n::NT, N::T) where {T <: QuantitativeNetwork, NT <: Number}
   ReturnType = T <: AbstractBipartiteNetwork ? BipartiteNetwork : UnipartiteNetwork
   return ReturnType(n .< N.edges, _species_objects(N)...)
end

function Base.:<(n::NT, N::T) where {T <: ProbabilisticNetwork, NT <: AbstractFloat}
   ReturnType = T <: AbstractBipartiteNetwork ? BipartiteNetwork : UnipartiteNetwork
   return ReturnType(n .< N.edges, _species_objects(N)...)
end

function Base.:≤(N::T, n::NT) where {T <: QuantitativeNetwork, NT <: Number}
   ReturnType = T <: AbstractBipartiteNetwork ? BipartiteNetwork : UnipartiteNetwork
   return ReturnType(N.edges .≤ n, _species_objects(N)...)
end

function Base.:≤(N::T, n::NT) where {T <: ProbabilisticNetwork, NT <: AbstractFloat}
   ReturnType = T <: AbstractBipartiteNetwork ? BipartiteNetwork : UnipartiteNetwork
   return ReturnType(N.edges .≤ n, _species_objects(N)...)
end

function Base.:≤(n::NT, N::T) where {T <: QuantitativeNetwork, NT <: Number}
   ReturnType = T <: AbstractBipartiteNetwork ? BipartiteNetwork : UnipartiteNetwork
   return ReturnType(n .≤ N.edges, _species_objects(N)...)
end

function Base.:≤(n::NT, N::T) where {T <: ProbabilisticNetwork, NT <: AbstractFloat}
   ReturnType = T <: AbstractBipartiteNetwork ? BipartiteNetwork : UnipartiteNetwork
   return ReturnType(n .≤ N.edges, _species_objects(N)...)
end


text:
"""
$(SIGNATURES)
Returns a summary list of all landmarks for a given robot and session.
"""
function getLandmarks(robotId::String, sessionId::String)::Vector{NodeResponse}
    landmarkList = filter(n -> occursin(r"[l][0-9]+", n.label), getVariables().nodes)
    return landmarkList
end

"""
$(SIGNATURES)
Returns a summary list of all landmarks for a given robot and session.
"""
function getLandmarks()::Vector{NodeResponse}
    config = getGraffConfig()
    if config == nothing
        error("Graff config is not set, please call setGraffConfig with a valid configuration.")
    end
    if config.robotId == "" || config.sessionId == ""
        error("Your config doesn't have a robot or a session specified, please attach your config to a valid robot or session by setting the robotId and sessionId fields. Robot = $(config.robotId), Session = $(config.sessionId)")
    end

    return getLandmarks(config.robotId, config.sessionId)
end

"""
$(SIGNATURES)
BETA method - better if run on the server.
Gets the estimates (if available) for a session or a filtered list if one is provided.
"""
function getEstimates(robotId::String, sessionId::String, nodes::Union{Vector{Union{String, Symbol}}, Vector{NodeResponse}, Nothing}=nothing)::Dict{String, Union{Array{Float64}, Nothing}}
    if nodes == nothing
        # Easy solve
        nodes = getVariables(robotId, sessionId)
        nodes = nodes.nodes
        return Dict(map(n -> n.label, nodes) .=> map(n -> n.mapEst, nodes))
    end
    if typeof(nodes) == Vector{NodeResponse}
        nodes = map(n -> n.label, nodes)
    end
    # Cleaning up if symbols
    nodes = map(n -> String(n), nodes)
    # Now do a getVariables
    respNodes = getVariables(robotId, sessionId).nodes
    respDict = Dict(map(n -> n.label, respNodes) .=> map(n -> n.mapEst, respNodes))
    # And map them - keeping order.
    ests = Dict{String, Union{Array{Any}, Nothing}}()
    for n in nodes
        if haskey(respDict, n)
            ests[n] = respDict[n]
        end
    end
    return ests
end

"""
$(SIGNATURES)
BETA method - better if run on the server.
Gets the estimates (if available) for a session or a filtered list if one is provided.
"""
function getEstimates(nodes::Union{Vector{Union{String, Symbol}}, Vector{NodeResponse}, Nothing}=nothing)::Dict{String, Union{Array{Float64}, Nothing}}
    config = getGraffConfig()
    if config == nothing
        error("Graff config is not set, please call setGraffConfig with a valid configuration.")
    end
    if config.robotId == "" || config.sessionId == ""
        error("Your config doesn't have a robot or a session specified, please attach your config to a valid robot or session by setting the robotId and sessionId fields. Robot = $(config.robotId), Session = $(config.sessionId)")
    end

    return getEstimates(config.robotId, config.sessionId, nodes)
end

ext: jl


text:
using Caesar, DelimitedFiles
using Gadfly, Cairo, Fontconfig
# using MAT
# using AbstractPlotting, Makie

using DocStringExtensions
include(joinpath(@__DIR__,"..","..","examples","marine","asv","kayaks","slamUtils.jl"))

fFloor = 250;
fCeil = 1750;
nFFT_czt = 1024;
fSampling = 37500.0
nPhones = 10;
soundSpeed = 1481;
azimuthDivs = 180;
azimuths = range(0,360,length=azimuthDivs)*pi/180;

# Load test data of length 10
# winstart = 0;
# posFile = joinpath(ENV["HOME"],"data","sas","test_array_positions.csv");
# arrayElemPos = readdlm(posFile,',',Float64,'\n')
# arrayElemPos = arrayElemPos[:,1:2]
# dataFile = joinpath(ENV["HOME"],"data","sas","test_array_waveforms.csv");
# rawWaveData = readdlm(dataFile,',',Float64,'\n')

#Load Experimental Data - 10 frames
# winstart = 160;
# rawWaveData = zeros(8000,nPhones);
# arrayElemPos = zeros(nPhones,2);
# for ele in winstart:winstart+nPhones-1
#     dataFile = joinpath(ENV["HOME"],"data", "sas", "06_20_sample","waveform$(ele).csv");
#     posFile = joinpath(ENV["HOME"],"data", "sas", "06_20_sample","nav$(ele).csv");
#     #dataFile = joinpath(ENV["HOME"],"liljondir", "kayaks", "20_gps_pos","waveform$(ele).csv");
#     #posFile = joinpath(ENV["HOME"],"liljondir", "kayaks", "20_gps_pos","nav$(ele).csv");
#     tempRead = readdlm(dataFile,',',Float64,'\n') #first element only
#     rawWaveData[:,ele-winstart+1] = adjoint(tempRead[1,:]);
#     tempRead = readdlm(posFile,',',Float64,'\n');
#     arrayElemPos[ele-winstart+1,:] = tempRead;
# end
# beacongt = [17.0499;1.7832];

# chirpFile = joinpath(ENV["HOME"],"data","sas","chirp250.txt");

# Synthetic data
winstart = 5;
posFile = joinpath(ENV["HOME"],"data","sas","syntheticArray.csv");
arrayElemPos = readdlm(posFile,',',Float64,'\n')
arrayElemPos = arrayElemPos[:,1:2]
dataFile = joinpath(ENV["HOME"],"data","sas","synthetic_data.csv");
rawWaveData = readdlm(dataFile,',',Float64,'\n')
chirpFile = joinpath(ENV["HOME"],"data","sas","chirp_signal_synth.csv");
nPhones = 11;
beacongt = [100,125];

chirpIn = readdlm(chirpFile,',',Float64,'\n')
FFTfreqs = collect(LinRange(fFloor,fCeil,nFFT_czt))

#Leave One Out Full CBF
leaveout = 6;
elementset = setdiff(Int[1:nPhones;],[leaveout;]);
arrayLIE = arrayElemPos[elementset,:];

cfg = CBFFilterConfig(fFloor,fCeil,nFFT_czt,nPhones-1,azimuths,soundSpeed,FFTfreqs)
myCBF = zeros(Complex{Float64}, getCBFFilter2Dsize(cfg));
lm = zeros(2);
dataTempHolder = zeros(Complex{Float64},nPhones-1,nFFT_czt)
delaysHolder = FFTfreqs;
@time constructCBFFilter2D!(cfg, arrayLIE, myCBF, lm, dataTempHolder)

# MF and CZT
w = exp(-2im*pi*(fCeil-fFloor)/(nFFT_czt*fSampling))
a = exp(2im*pi*fFloor/fSampling)

chirpIn = readdlm(chirpFile,',',Float64,'\n')

#Matched Filter on Data In
nFFT_full = nextpow(2,size(rawWaveData,1))  # MF
mfData = zeros(Complex{Float64}, nFFT_full, nPhones)
mf = prepMF(chirpIn,nFFT_full,nPhones) # MF
mf(rawWaveData,mfData) # MF

cztData = zeros(Complex{Float64}, nFFT_czt,nPhones)
filterCZT = prepCZTFilter(nFFT_full,nPhones,w,nFFT_czt,a)
filterCZT(mfData,cztData)

#Try CBF over LIE
dataOut = zeros(length(azimuths));
cztDataLIE = cztData[:,elementset];
temp1 = zeros(Complex{Float64},nFFT_czt);
temp2 = zeros(Complex{Float64},size(cztDataLIE));
@time CBF2D_DelaySum!(cfg, cztDataLIE, dataOut,temp1,temp2,myCBF)

pl = Gadfly.plot(
 Gadfly.layer(x=azimuths,y=real(dataOut), Geom.line),
 Gadfly.Coord.cartesian(xmin=0, xmax=2*pi),
 Guide.ylabel(nothing),Guide.xlabel("Azimuths (rad)")
 ) ; pl |> PDF("/tmp/test_CBF_LIE.pdf")

# LIE CBF - Assume Known Source, azi phase shift
dataOutLIE = zeros(Complex{Float64}, 1);
dx = arrayLIE[1,1] - beacongt[1];
dy = arrayLIE[1,2] - beacongt[2];
azi = atan(dy,dx);
cfgLIE = CBFFilterConfig(fFloor,fCeil,nFFT_czt,nPhones-1,[azi;],soundSpeed,FFTfreqs)
myCBFLIE = zeros(Complex{Float64}, getCBFFilter2Dsize(cfgLIE));
@time constructCBFFilter2D!(cfgLIE, arrayLIE, myCBFLIE, lm, dataTempHolder)

#Try Correlation
cztDataHolder = zeros(Complex{Float64},size(cztDataLIE));
@time liebf!(dataOutLIE, cztDataLIE, myCBFLIE, 1, cztDataHolder, normalize=true)

#Try singlephaseshift
phaseshiftLOO = cztData[:,leaveout];
@time phaseShiftSingle!(lm, cfgLIE, azi, arrayElemPos[leaveout,1:2], phaseshiftLOO)

#Visualize residual at correct look angle
variation = -5:0.1:5;
allshifts = zeros(Complex{Float64},length(variation),length(variation),length(cztData[:,leaveout]));
dataOutRes = zeros(Complex{Float64}, nFFT_czt);

dxArr = arrayElemPos[1,1] - beacongt[1];
dyArr = arrayElemPos[1,2] - beacongt[2];
aziArr = atan(dyArr/dxArr)

cfgLIE = CBFFilterConfig(fFloor,fCeil,nFFT_czt,nPhones-1,[aziArr;],soundSpeed,FFTfreqs)
myCBFLIE = zeros(Complex{Float64}, getCBFFilter2Dsize(cfgLIE));
constructCBFFilter2D!(cfgLIE, arrayLIE, myCBFLIE, lm, dataTempHolder,delaysHolder)
cztDataHolder = zeros(Complex{Float64},size(cztDataLIE));
liebf!(dataOutRes, cztDataLIE, myCBFLIE, 1, cztDataHolder, normalize=true)

for xInd in 1:length(variation)
    for yInd in 1:length(variation)

        dx = arrayElemPos[leaveout,1] - beacongt[1] + variation[xInd];
        dy = arrayElemPos[leaveout,2] - beacongt[2] + variation[yInd];
        azi = atan(dy,dx);
        phaseshiftLOO = cztData[:,leaveout];
        phaseShiftSingle!(lm, cfgLIE, azi, arrayElemPos[leaveout,1:2], phaseshiftLOO)
        allshifts[xInd,yInd,:] = phaseshiftLOO .+ dataOutRes;
    end
end

# tmpfile = "/tmp/res$(winstart).mat";
# file = matopen(tmpfile, "w")
# write(file, "allshifts", allshifts)
# close(file)

#Visualize residual over many look angles
variation = range(-5,5,length=200);
uvec = arrayElemPos[leaveout,1:2]-beacongt
uvec ./= norm(uvec)
rvec = R(pi/2)*uvec;

allshiftsAziX = zeros(Complex{Float64},length(azimuths),length(variation),length(cztData[:,leaveout]));
dataOutRes = zeros(Complex{Float64}, nFFT_czt);
allshiftsAziY = zeros(Complex{Float64},length(azimuths),length(variation),length(cztData[:,leaveout]));
dataOutRes = zeros(Complex{Float64}, nFFT_czt);
allshiftsAzilin = zeros(Complex{Float64},length(azimuths),length(variation),length(cztData[:,leaveout]));
dataOutRes = zeros(Complex{Float64}, nFFT_czt);

zdx = zeros(length(azimuths),length(variation)); #plot vars
zdy = zeros(length(azimuths),length(variation));
zlin = zeros(length(azimuths),length(variation));

dataOutHolder = zeros(Complex{Float64}, nFFT_czt);

for aziInd in 1:length(azimuths)
    cfgLIE = CBFFilterConfig(fFloor,fCeil,nFFT_czt,nPhones-1,[azimuths[aziInd];],soundSpeed,FFTfreqs)
    myCBFLIE = zeros(Complex{Float64}, getCBFFilter2Dsize(cfgLIE));
    constructCBFFilter2D!(cfgLIE, arrayLIE, myCBFLIE, lm, dataTempHolder)
    cztDataHolder = zeros(Complex{Float64},size(cztDataLIE));
    liebf!(dataOutRes, cztDataLIE, myCBFLIE, 1, cztDataHolder, normalize=true)

    for xInd in 1:length(variation)
            newPos = arrayElemPos[leaveout,1:2]+[variation[xInd];0];
            dx = newPos[1] - beacongt[1];
            dy = newPos[2] - beacongt[2];
            azi = atan(dy,dx);
            phaseshiftLOO = cztData[:,leaveout];
            phaseShiftSingle!(lm, cfgLIE, azi,newPos , phaseshiftLOO)
            zdx[aziInd,xInd] = 0.0;
            copy!(dataOutHolder,dataOutRes)
            @inbounds for i in 1:length(phaseshiftLOO)
              dataOutHolder[i] += phaseshiftLOO[i]
              zdx[aziInd,xInd] += norm(dataOutHolder[i])
            end
    end

    for yInd in 1:length(variation)
            # newPos = arrayElemPos[leaveout,1:2]+[0;variation[yInd]];
            dlin = variation[yInd]*rvec;
            newPos = arrayElemPos[leaveout,1:2]+dlin;

            dx = newPos[1] - beacongt[1];
            dy = newPos[2] - beacongt[2];
            azi = atan(dy,dx);
            phaseshiftLOO = cztData[:,leaveout];
            phaseShiftSingle!(lm, cfgLIE, azi,newPos, phaseshiftLOO)
            zdy[aziInd,yInd] = 0.0;
            copy!(dataOutHolder,dataOutRes)
            @inbounds for i in 1:length(phaseshiftLOO)
              dataOutHolder[i] += phaseshiftLOO[i]
              zdy[aziInd,yInd] += norm(dataOutHolder[i])
            end
    end

    for linInd in 1:length(variation)
            dlin = variation[linInd]*uvec;
            newPos = arrayElemPos[leaveout,1:2]+dlin;
            dx = newPos[1] - beacongt[1];
            dy = newPos[2] - beacongt[2];
            azi = atan(dy,dx);
            phaseshiftLOO = cztData[:,leaveout];
            phaseShiftSingle!(lm, cfgLIE, azi,newPos, phaseshiftLOO)
            zlin[aziInd,linInd] = 0.0;
            copy!(dataOutHolder,dataOutRes)
            @inbounds for i in 1:length(phaseshiftLOO)
              dataOutHolder[i] += phaseshiftLOO[i]
              zlin[aziInd,linInd] += norm(dataOutHolder[i])
          end
    end
end

# scene = AbstractPlotting.surface(azimuths, variation, z./maximum(z))
# Makie.save("plot.png", scene)
#
# tmpfile = "/tmp/resAzis$(winstart).mat";
# file = matopen(tmpfile, "w")
# write(file, "zdx", zdx)
# write(file, "zdy", zdy)
# write(file, "zlin", zlin)
# close(file)


text:
## Optimisation Benchmarks from Jamil, Yang 2013
## =============================================

using DynamicAnalysis
import DynamicAnalysis: benchmark

using Lens
using MathProgBase
using JuMP
using DReal
using NLopt

e_ = convert(Float64, e)
pi_ = convert(Float64, π)

type OptProb <: Problem
  name::String
  n::Int
  lb::Float64
  ub::Float64
  sense::Symbol
  obj::Expr
  min::Float64
end

type OptAlgo <: Algorithm
  solver::MathProgBase.MathProgSolverInterface.AbstractMathProgSolver
  capture::Vector{Symbol}
end

# Benchmark an algorithm against a optimization problem
function benchmark(a::OptAlgo, p::OptProb)
  value, results = capture(()->eval(gen_f(a,p)), a.capture)
  reset_global_ctx!()
  results
end

function gen_f(a::OptAlgo, p::OptProb)
  quote
    m = Model(solver = $(a.solver))
    n = $(p.n)
    @defVar(m, $(p.lb) <= x[1:$(p.n)] <= $(p.ub))
    @setNLObjective(m, :Min , $(p.obj))
    @time solve(m)
    @show obj_val = getObjectiveValue(m)
    lens(:objective_val, obj_val)
  end
end

## Algorithms
## ==========
dreal = OptAlgo(DRealSolver(precision = 0.001), [:objective_val])
mma = OptAlgo(NLoptSolver(algorithm=:LD_MMA), [:objective_val])
cobyla = OptAlgo(NLoptSolver(algorithm=:LN_COBYLA), [:objective_val])

## Problems
## ========
d = 5
ackley1 = OptProb("ackley1", d, -35.0, 35.0, :Min,
  :(-20*exp(-0.02*sqrt((1/n)*sum{x[i]^2, i=1:n})) -
    exp((1/n)*sum{cos(2*pi_*x[i]), i=1:n}) + 20 + e), 0.0)

ackley2 = OptProb("ackley2", 2, -32.0, 32.0, :Min,
  :(200*exp(-0.02*sqrt(x[1]^2 + x[2]^2))), -200.0)

ackley3 = OptProb("ackley3", 2, -32.0, 32.0, :Min,
  :(200*exp(-0.02*sqrt(x[1]^2 + x[2]^2)) + 5*exp(cos(3*x[1])+sin(3*x[2]))), -219.1418)

ackley4 = OptProb("ackley4", d, -35.0, 35.0, :Min,
  :(sum{exp(-0.2*sqrt(x[1]^2 + x[2]^2) + 3*(cos(2x[1]) + sin(2x[1+1]))), i=1:n-1}), -3.917275)

alpine7 = OptProb("alpine7", d, 0.0, 10.0, :Min,
  :(prod{sqrt(x[i])*sin(x[i]), i=1:n}), 2.808)

# Skip due to unequal constraints
# brad8 = OptProb("brad8", )

# 9. Skip due to abs

beale10 = OptProb("beale10", 2, -4.5, 4.5, :Min,
  :((1.5 - x[1] +x[1]*x[2])^2 + (2.25 - x[1] + x[1]*x[2]^2)^2 + (2.625 - x[1] + x[1]*x[2]^3)^2), 1.0)

biggs11 = OptProb("biggs11", 2, 0.0, 20.0, :Min,
  :(sum{(exp(-0.1*i*x[1]) - 5*exp(-0.1*i*x[2]) - (exp(-0.1*i) - 5*exp(i)))^2, i=1:n}), 0.0)

biggs12 = OptProb("biggs12", 3, 0.0, 20.0, :Min,
  :(sum{(exp(-0.1*i*x[1]) - x[3]*exp(-0.1*i*x[2]) - (exp(-0.1*i) - 5*exp(i)))^2, i=1:10}), 0.0)

biggs13 = OptProb("biggs13", 4, 0.0, 20.0, :Min,
  :(sum{(x[3]*exp(-0.1*i*x[1]) - x[4]*exp(-0.1*i*x[2]) - (exp(-0.1*i) - 5*exp(i)))^2, i=1:10}), 0.0)

biggs14 = OptProb("biggs14", 6, -20.0, 20.0, :Min,
  :(sum{(x[3]*exp(-0.1i*x[1]) - x[4]*exp(-0.1i*x[2]) + 3*exp(-0.1i*x[5]) - (exp(-0.1i) - 5*exp(i) + 3*exp(-0.4i)))^2, i=1:11}), 0.0)

biggs15 = OptProb("biggs15", 6, -20.0, 20.0, :Min,
  :(sum{(x[3]*exp(-0.1i*x[1]) - x[4]*exp(-0.1i*x[2]) + x[6]*exp(-0.1i*x[5]) - (exp(-0.1i) - 5*exp(i) + 3*exp(-0.4i)))^2, i=1:13}), 0.0)

bird16 = OptProb("bird16", 2, -2*pi_, 2*pi_, :Min,
  :(sin(x[1])*exp((1-cos(x[2])^2)) +cos(x[2])*exp((1-sin(x[1]))^2) + (x[1] - x[2])^2), -106.764537)

bohachevsky17 = OptProb("bohachevsky17", 2, -100.0, 100.0, :Min,
  :(x[1]^2 + 2*x[2]^2 - 0.3*cos(3*pi_*x[1]) - 0.4*cos(4*pi_*x[2]) + 0.7), 0.0)

bohachevsky18 = OptProb("bohachevsky18", 2, -100.0, 100.0, :Min,
  :(x[1]^2 + 2*x[2]^2 - 0.3*cos(3*pi_*x[1]) * 0.4*cos(4*pi_*x[2]) + 0.3), 0.0)

bohachevsky19 = OptProb("bohachevsky19", 2, -100.0, 100.0, :Min,
  :(x[1]^2 +2*x[2]^2 - 0.3*cos(3*pi_*x[1] + 4*pi_*x[2]) + 0.3), 0.0)

booth20 = OptProb("booth20", 2, -10.0, 10.0, :Min,
  :((x[1] + 2*x[2] - 7)^2 + (2x[1] + x[2] -5)^2), 0.0)

# Uneven Bounds
# box_betts21 = OptProb("box_betts21", )

#uneven bounds branin

branin23 = OptProb("branin23", 3, -5.0, 15.0, :Min,
  :((x[2] - (5.1x[1]^2)/(4*pi_^2) + (5x[1])/pi_ - 6)^2 +
    10*(1 - 1/(8pi_))*cos(x[1])*cos(x[2])*log(x[1]^2+x[2]^2+1) + 10), 5.559037)

brent24 = OptProb("brent24", 2, -10.0, 10.0, :Min,
  :((x[1] + 10)^2 + (x[2] + 10)^2 +exp(-x[1]^2 - x[2]^2)), 0.0)

brown25 = OptProb("brown", 4, -4.0, 4.0, :Min,
  :(sum{x[i]^(2*x[i+1]^2+1) + x[i+1]^(2*x[i]^2+1), i=1:n-1}),
  0.0)

# bukin26 = OptProb("bukin26", 2, )
camel29 = OptProb("camel29", 2, -5.0, 5.0, :Min,
  :(2x[1]^2 - 1.05x[1]^4 + (x[1]^6)/6 +x[1]*x[2] +x[2]^2), 0.0)

camel30 = OptProb("camel30", 2, -5.0, 5.0, :Min,
  :((4-2.1x[1]^2 + (x[1]^4)/3)*x[1]^2 +x[1]*x[2] + (4*x[2]^2 - 4) * x[2]^2), -1.0316)

# Contains Floor

chichinadze33 = OptProb("chichinadze33", 2, -30.0, 30.0, :Min,
  :(x[1]^2 - 12x[1] + 11 + 10*cos(pi_*x[1]/2) + 8*sin(5pi_*x[1]/2) -
    (1/5)^0.5*exp(-0.5*(x[2] - 0.5)^2)), -43.3159)

chungreynolds34 = OptProb("chungreynolds34", d, -100.0, 100.0, :Min,
  :((sum{x[i]^2, i=1:n})^2), 0.0)

#35 .. matrices...

# abs
# corana37 = OptProb("corana37", d, -500.0, 500.0, :Min,
#   )

# ifelset = OptProb("ifelset", 1, 0.0, 1.0, :Min,
#   :(ifelse(x[1] > 0.5), x, 2x), 10)

cosinemix38 = OptProb("cosinemix38", 4, -1.0, 1.0, :Min,
  :(-0.1 * sum{cos(5pi_*x[i]), i=1:n} - sum{x[i]^2, i=1:n}), 0.4)

csendes40 = OptProb("csendes40", d, -1.0, 1.0, :Min,
  :(sum{x[i]^6 * (2 + sin(1/x[i])), i=1:n}), 0)

cube41 = OptProb("cube41", 2, -10.0, 10.0, :Min,
  :(100*(x[2] - x[1]^3)^2 + (1-x[1])^2), 0.0)

deb43 = OptProb("deb43", d, -1.0, 1.0, :Min,
  :(-1/n * sum{sin(5*pi_*x[i])^6, i=1:n}), -Inf)

deb44 = OptProb("deb44", d, -1.0, 1.0, :Min,
  :(-1/n * sum{sin(5pi_*(x[1]^.75 - 0.05))^6, i=1:n}), -Inf)

deckkers_arts45 = OptProb("deckkers_arts45", 2, -20.0, 20.0, :Min,
  :(10^5*x[1]^2 + x[2]^2 - (x[1]^2 + x[2]^2)^2 + 10.0^-5 * (x[1]^2 + x[2]^2)^4), -2477)

rosenbrock105 = OptProb("rosenbrock105", d, -30.0, 30.0, :Min,
  :(sum{100*(x[i+1] - x[i]^2)^2 + (x[i] -1)^2, i=1:n-1}), 0.0)

ripple = OptProb("ripple", 2, 0.0, 1.0, :Min, :(sum{x[i]^(2*x[i+1]^2+1) + x[i+1]^(2*x[i]^2+1), i=1:n-1}), -Inf)
whitley = OptProb("whitley", 5, -100.0, 100.0, :Min, :(sum{((100*(x[i]^2 -x[j])^2 + (1 - x[j])^2)^2)/4000 - cos(100*(x[i]^2 - x[j])^2 + (1 - x[j])^2) + 1, i = 1:n, j = 1:n}), -1.)

## Run Experiments
problems = [ackley1, ackley2, ackley3, ackley4, alpine7, beale10, biggs11, brown, ripple, whitley]
algos = [mma, dreal, cobyla]

# results = record(algos,problems; runname = "kl",savedb=false,exceptions=true)

## 171. Xin-She
# xinshe171 = OptProb(5, -10.0, 10.0, :Min, :(exp(-sum{(x[i]/15)^(2*5), i=1:n}) - 2*exp(-sum{x[i]^2, i=1:n}), -1.0)


text:
#=
Visual illustration of the law of large numbers.

@author : Spencer Lyon <spencer.lyon@nyu.edu>

References
----------

Based off the original python file illustrates_lln.py
=#
using PyPlot
using Distributions

n = 100
srand(42)  # reproducible results

# == Arbitrary collection of distributions == #
distributions = {"student's t with 10 degrees of freedom" => TDist(10),
                 "beta(2, 2)" => Beta(2.0, 2.0),
                 "lognormal LN(0, 1/2)" => LogNormal(0.5),
                 "gamma(5, 1/2)" => Gamma(5.0, 2.0),
                 "poisson(4)" => Poisson(4),
                 "exponential with lambda = 1" => Exponential(1)}

# == Create a figure and some axes == #
num_plots = 3
fig, axes = plt.subplots(num_plots, 1, figsize=(10, 10))

bbox = [0., 1.02, 1., .102]
legend_args = {:ncol => 2,
               :bbox_to_anchor => bbox,
               :loc => 3,
               :mode => "expand"}
subplots_adjust(hspace=0.5)

for ax in axes
    dist_names = collect(keys(distributions))
    # == Choose a randomly selected distribution == #
    name = dist_names[rand(1:length(dist_names))]
    dist = pop!(distributions, name)

    # == Generate n draws from the distribution == #
    data = rand(dist, n)

    # == Compute sample mean at each n == #
    sample_mean = Array(Float64, n)
    for i=1:n
        sample_mean[i] = mean(data[1:i])
    end

    # == Plot == #
    ax[:plot](1:n, data, "o", color="grey", alpha=0.5)
    axlabel = LaTeXString("\$\\bar{X}_n\$ for \$X_i \\sim\$ $name")
    ax[:plot](1:n, sample_mean, "g-", lw=3, alpha=0.6, label=axlabel)
    m = mean(dist)
    ax[:plot](1:n, ones(n)*m, "k--", lw=1.5, label=L"$\mu$")
    ax[:vlines](1:n, m, data, lw=0.2)
    ax[:legend](;legend_args...)
end


text:
# helper to create fake data from scores such that the detectors return the scores again
to_table(points) = reshape(repeat(points, length(points)), (length(points), length(points))) |> table
const base = [1., 2, 3]
const test_scores = [base, # train and test equal
                     [1.,2], # len(test) < len(train)
                     [-Inf, 0., 1, 2, 3, 4, Inf], # len(test) > len(train)
                     [1.,1,1],  # len(test) = len(train)
                     [2.,2,2],
                     [3.,3,3]]
const scores0, scores1, scores2, scores3, scores4, scores5 = map(t -> (base, t), test_scores)
const X_test0, X_test1, X_test2, X_test3, X_test4, X_test5 = map(to_table, test_scores)
const X = to_table(base)
const y = to_categorical(["normal", "normal", "outlier"])

# # fake detectors
supervised = MinimalSupervisedDetector()
unsupervised = MinimalUnsupervisedDetector()

# # fake machines
u = machine(unsupervised, X) |> fit!
s = machine(supervised, X, y) |> fit!
raw_machines = [u, s]

uc = machine(CompositeDetector(unsupervised), X) |> fit!
sc = machine(CompositeDetector(supervised), X, y) |> fit!
up = machine(ProbabilisticDetector(unsupervised), X) |> fit!
sp = machine(ProbabilisticDetector(supervised), X, y) |> fit!
ud = machine(DeterministicDetector(unsupervised), X) |> fit!
sd = machine(DeterministicDetector(supervised), X, y) |> fit!
composite_machines = [uc, sc, up, sp, ud, sd]

upc = machine(ProbabilisticDetector(CompositeDetector(unsupervised)), X) |> fit!
spc = machine(ProbabilisticDetector(CompositeDetector(supervised)), X, y) |> fit!
udc = machine(DeterministicDetector(CompositeDetector(unsupervised)), X) |> fit!
sdc = machine(DeterministicDetector(CompositeDetector(supervised)), X, y) |> fit!
wrapped_machines = [upc, spc, udc, sdc]

usc = machine(CompositeDetector(u=unsupervised, s=supervised), X, y) |> fit!
usp = machine(ProbabilisticDetector(u=unsupervised, s=supervised), X, y) |> fit!
usd = machine(DeterministicDetector(u=unsupervised, s=supervised), X, y) |> fit!
combined_machines = [usc, usp, usd]

machines = [composite_machines..., wrapped_machines..., combined_machines...]
probabilistic_machines = [up, sp, upc, spc, usp]
deterministic_machines = [ud, sd, udc, sdc, usd]

# create a simple detector learning network
Xs = source()
ys = source()
uscores = augmented_transform(machine(unsupervised, Xs), Xs)
sscores = augmented_transform(machine(supervised, Xs, ys), Xs)
upscores = predict(machine(ProbabilisticTransformer()), uscores)
spscores = predict(machine(ProbabilisticTransformer()), sscores)
udscores = predict(machine(DeterministicTransformer()), uscores)
sdscores = predict(machine(DeterministicTransformer()), sscores)

# TODO: The empty source for ys in unsupervised models is currently necessary to enable evaluation, but
# we could fix this in MLJBase
upmach = machine(OD.ProbabilisticUnsupervisedDetector(), Xs, source(); predict=upscores)
spmach = machine(OD.ProbabilisticSupervisedDetector(), Xs, ys; predict=spscores)
udmach = machine(OD.DeterministicUnsupervisedDetector(), Xs, source(); predict=udscores)
sdmach = machine(OD.DeterministicSupervisedDetector(), Xs, ys; predict=sdscores)

@testset "normalization, combination and classification" begin
    minmax_test = scores -> last(scale_minmax(scores)) # extract test scores
    unify_test = scores -> last(scale_unify(scores)) # extract test scores
    classify = scores -> last(classify_quantile(DEFAULT_THRESHOLD)(scale_minmax(scores))) # extract test scores

    @testset "scores" begin
        raw_proba(detector, data) = from_univariate_finite.(predict(detector, data))
        to_scores(data) = [[transform(m, data) for m in machines]..., # transform works for all machines
                           [raw_proba(m, data) for m in probabilistic_machines]...]
        test_scoring(scores, data, labels) = @test all((labels,) .== [minmax_test(scores), to_scores(data)...])

        # Tests expected results on simple score vectors
        test_scoring(scores0, X_test0, [0.0, 0.5, 1.0])
        test_scoring(scores1, X_test1, [0.0, 0.5])
        test_scoring(scores2, X_test2, [0.0, 0.0, 0.0, 0.5, 1.0, 1.0, 1.0])
        test_scoring(scores3, X_test3, [0.0, 0.0, 0.0])
        test_scoring(scores4, X_test4, [0.5, 0.5, 0.5])
        test_scoring(scores5, X_test5, [1.0, 1.0, 1.0])

        # Min-Max normalization and unify should not change sort order
        @test sortperm(minmax_test(scores0)) == sortperm(unify_test(scores0))
        @test sortperm(minmax_test(scores1)) == sortperm(unify_test(scores1))
        @test sortperm(minmax_test(scores2)) == sortperm(unify_test(scores2))
        @test sortperm(minmax_test(scores3)) == sortperm(unify_test(scores3))
        @test sortperm(minmax_test(scores4)) == sortperm(unify_test(scores4))
        @test sortperm(minmax_test(scores5)) == sortperm(unify_test(scores5))
    end

    @testset "classification" begin
        raw_class(detector, data) = from_categorical.(predict(detector, data))
        to_classes(data) = [[predict(m, data) for m in deterministic_machines]...,
                            [raw_class(m, data) for m in deterministic_machines]...]
        test_classification(scores, data, labels) = @test all((labels,) .== [classify(scores), to_classes(data)...])

        # Test expected classification results with given threshold
        test_classification(scores0, X_test0, ["normal", "normal", "outlier"])
        test_classification(scores1, X_test1, ["normal", "normal"])
        test_classification(scores2, X_test2, ["normal", "normal", "normal", "normal", "outlier", "outlier", "outlier"])
        test_classification(scores3, X_test3, ["normal", "normal", "normal"])
        test_classification(scores4, X_test4, ["normal", "normal", "normal"])
        test_classification(scores5, X_test5, ["outlier", "outlier", "outlier"])
    end

    @testset "combination of scores" begin
        @test combine_mean(minmax_test.([scores3, scores4])...) == [0.25, 0.25, 0.25]
        @test combine_mean(minmax_test.([scores3, scores5])...) == [0.5, 0.5, 0.5]
        @test combine_mean(minmax_test.([scores4, scores5])...) == [0.75, 0.75, 0.75]

        @test combine_median(minmax_test.([scores3, scores4])...) == [0.25, 0.25, 0.25]
        @test combine_median(minmax_test.([scores3, scores5])...) == [0.5, 0.5, 0.5]
        @test combine_median(minmax_test.([scores4, scores5])...) == [0.75, 0.75, 0.75]

        @test combine_max(minmax_test.([scores3, scores4])...) == [0.5, 0.5, 0.5]
        @test combine_max(minmax_test.([scores3, scores5])...) == [1.0, 1.0, 1.0]
        @test combine_max(minmax_test.([scores4, scores5])...) == [1.0, 1.0, 1.0]
    end
end

@testset "wrappers result in expected types" begin
    # unsupervised only
    @test CompositeDetector(unsupervised) isa UnsupervisedDetectorComposite
    @test CompositeDetector(detector=unsupervised) isa UnsupervisedDetectorComposite
    @test CompositeDetector(uns1=unsupervised, uns2=unsupervised) isa UnsupervisedDetectorComposite

    @test ProbabilisticDetector(unsupervised) isa ProbabilisticUnsupervisedDetectorComposite
    @test ProbabilisticDetector(detector=unsupervised) isa ProbabilisticUnsupervisedDetectorComposite
    @test ProbabilisticDetector(uns1=unsupervised, uns2=unsupervised) isa ProbabilisticUnsupervisedDetectorComposite

    @test DeterministicDetector(unsupervised) isa DeterministicUnsupervisedDetectorComposite
    @test DeterministicDetector(detector=unsupervised) isa DeterministicUnsupervisedDetectorComposite
    @test DeterministicDetector(uns1=unsupervised, uns2=unsupervised) isa DeterministicUnsupervisedDetectorComposite

    # supervised only
    @test CompositeDetector(supervised) isa SupervisedDetectorComposite
    @test CompositeDetector(detector=supervised) isa SupervisedDetectorComposite
    @test CompositeDetector(sup1=supervised, sup2=supervised) isa SupervisedDetectorComposite

    @test ProbabilisticDetector(supervised) isa ProbabilisticSupervisedDetectorComposite
    @test ProbabilisticDetector(detector=supervised) isa ProbabilisticSupervisedDetectorComposite
    @test ProbabilisticDetector(sup1=supervised, sup2=supervised) isa ProbabilisticSupervisedDetectorComposite

    @test DeterministicDetector(supervised) isa DeterministicSupervisedDetectorComposite
    @test DeterministicDetector(detector=supervised) isa DeterministicSupervisedDetectorComposite
    @test DeterministicDetector(sup1=supervised, sup2=supervised) isa DeterministicSupervisedDetectorComposite

    # mixed supervised and unsupervised
    @test CompositeDetector(sup=supervised, uns=unsupervised) isa SupervisedDetectorComposite
    @test ProbabilisticDetector(sup=supervised, uns=unsupervised) isa ProbabilisticSupervisedDetectorComposite
    @test DeterministicDetector(sup=supervised, uns=unsupervised) isa DeterministicSupervisedDetectorComposite

    # wrapped composites
    @test CompositeDetector(CompositeDetector(unsupervised)) isa UnsupervisedDetectorComposite
    @test ProbabilisticDetector(CompositeDetector(unsupervised)) isa ProbabilisticUnsupervisedDetectorComposite
    @test DeterministicDetector(CompositeDetector(unsupervised)) isa DeterministicUnsupervisedDetectorComposite

    @test CompositeDetector(CompositeDetector(supervised)) isa SupervisedDetectorComposite
    @test ProbabilisticDetector(CompositeDetector(supervised)) isa ProbabilisticSupervisedDetectorComposite
    @test DeterministicDetector(CompositeDetector(supervised)) isa DeterministicSupervisedDetectorComposite
end

@testset "wrappers property access" begin
    normalization_strategy = scale_unify
    combination_strategy = combine_max
    classification_strategy = classify_quantile(0.5)

    for m in getproperty.(machines, :model)
        initial_normalization_strategy = m.normalize
        initial_combination_strategy = m.combine
        m.normalize = normalization_strategy
        m.combine =  combination_strategy
        @test m.normalize == normalization_strategy
        @test m.combine == combination_strategy
        # reset back to original value
        m.normalize = initial_normalization_strategy
        m.combine = initial_combination_strategy
    end

    for m in getproperty.(deterministic_machines, :model)
        initial_classification_strategy = m.classify
        m.classify = classification_strategy
        @test m.classify == classification_strategy
        # reset back to original value
        m.classify = initial_classification_strategy
    end

    # wrappers throw an error if a property does not exist
    for m in getproperty.(machines, :model)
        @test_throws ErrorException m.foo
        @test_throws ErrorException m.foo = "bar"
    end
end

@testset "erroneous wrapper calls" begin
    # wrappers do not work with models other than detectors
    static_model = MLJBase.WrappedFunction(identity)
    @test_throws MethodError CompositeDetector(static_model)
    @test_throws MethodError ProbabilisticDetector(static_model)
    @test_throws MethodError DeterministicDetector(static_model)

    # wrappers do not work with multiple unnamed detectors
    @test_throws ArgumentError CompositeDetector(unsupervised, supervised)
    @test_throws ArgumentError ProbabilisticDetector(unsupervised, supervised)
    @test_throws ArgumentError DeterministicDetector(unsupervised, supervised)

    # wrappers warn if both arguments and named arguments are provided
    @test_logs (:warn, r"Wrapping the single detector") CompositeDetector(unsupervised, s=supervised)
    @test_logs (:warn, r"Wrapping the single detector") ProbabilisticDetector(unsupervised, s=supervised)
    @test_logs (:warn, r"Wrapping the single detector") DeterministicDetector(unsupervised, s=supervised)
end

@testset "correct augmented_transform calls" begin
    test_a(m) = @test augmented_transform(m) == (m.report.scores, transform(m))
    test_b(m, X) = @test augmented_transform(m, X) == (m.report.scores, transform(m, X))
    test_c(m; rows=:) = @test augmented_transform(m; rows=rows) == (m.report.scores, transform(m; rows=rows))

    # make sure augmented_transform works as expected on all kinds of machines
    for m in [raw_machines..., machines...]
        test_a(m)
        test_b(m, X)
        test_c(m; rows=1:2)
    end
end

@testset "erroneous augmented_transform calls" begin
    u_not_fitted = machine(unsupervised, X)
    s_not_fitted = machine(supervised, X, y)
    c_not_fitted = machine(CompositeDetector(unsupervised), X)

    # not-yet-fitted machines
    @test_throws ErrorException augmented_transform(u_not_fitted)
    @test_throws ErrorException augmented_transform(s_not_fitted)
    @test_throws ErrorException augmented_transform(c_not_fitted)
end

@testset "transformers expected results" begin
    Xs = source(X)
    ys = source(y)

    # prepare learning network machines
    um = machine(MinimalUnsupervisedDetector(), Xs)
    sm = machine(MinimalSupervisedDetector(), Xs, ys)

    # prepare the transformers
    score_transformer = machine(ScoreTransformer())
    probabilistic_transformer = machine(ProbabilisticTransformer())
    deterministic_transformer = machine(DeterministicTransformer())

    # get the augmented scores
    u_scores = augmented_transform(um, Xs)
    s_scores = augmented_transform(sm, Xs)

    fit_transform(transformer, scores) = fit!(transform(transformer, scores))()
    fit_predict(transformer, scores) = fit!(predict(transformer, scores))()

    # test the transformed scores
    @test fit_transform(score_transformer, u_scores) isa OD.Scores
    @test fit_transform(score_transformer, s_scores) isa OD.Scores
    @test fit_transform(probabilistic_transformer, u_scores) isa OD.Scores
    @test fit_transform(probabilistic_transformer, s_scores) isa OD.Scores
    @test fit_transform(deterministic_transformer, u_scores) isa OD.Scores
    @test fit_transform(deterministic_transformer, s_scores) isa OD.Scores

    # test the predicted probabilities
    @test fit_predict(probabilistic_transformer, u_scores) isa UnivariateFiniteVector
    @test fit_predict(probabilistic_transformer, s_scores) isa UnivariateFiniteVector

    # test the predicted classes
    @test fit_predict(deterministic_transformer, u_scores) isa OD.Labels
    @test fit_predict(deterministic_transformer, s_scores) isa OD.Labels
end

@testset "evaluation works as expected" begin
    train, test = [1, 2], [3]
    evaluate_detector(detector) =
        @test evaluate(detector, X, y; measure=misclassification_rate, resampling=[(train, test)]).measurement[1] == 0

    # evaluation of wrapped detectors
    evaluate_detector(ProbabilisticDetector(unsupervised))
    evaluate_detector(DeterministicDetector(unsupervised))

    # create the detectors from the network
    @from_network upmach mutable struct CustomUnsupervisedProbabilisticDetector end
    @from_network spmach mutable struct CustomSupervisedProbabilisticDetector end
    @from_network udmach mutable struct CustomUnsupervisedDeterministicDetector end
    @from_network sdmach mutable struct CustomSupervisedDeterministicDetector end

    evaluate_detector(CustomUnsupervisedProbabilisticDetector())
    evaluate_detector(CustomSupervisedProbabilisticDetector())
    evaluate_detector(CustomUnsupervisedDeterministicDetector())
    evaluate_detector(CustomSupervisedDeterministicDetector())

    # TODO: Enable composites with custom detectors after https://github.com/JuliaAI/MLJBase.jl/pull/644 is merged
    # CompositeDetector(CustomUnsupervisedProbabilisticDetector)
end


text:
#=
Copyright (c) 2015, Intel Corporation

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

    * Redistributions of source code must retain the above copyright notice,
      this list of conditions and the following disclaimer.
    * Redistributions in binary form must reproduce the above copyright
      notice, this list of conditions and the following disclaimer in the
      documentation and/or other materials provided with the distribution.
    * Neither the name of Intel Corporation nor the names of its contributors
      may be used to endorse or promote products derived from this software
      without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
=#

using Latte
srand(12345)

net = Net(100)
data, label = HDF5DataLayer(net, "data/train.txt", "data/test.txt")
conv1    = ConvolutionLayer(:conv1, net, data, 20, 5, 1, 1)
relu1    = ReLULayer(:relu1, net, conv1)
pool1    = MaxPoolingLayer(:pool1, net, relu1, 2, 2, 0)
conv2    = ConvolutionLayer(:conv2, net, pool1, 50, 5, 1, 1)
relu2    = ReLULayer(:relu2, net, conv2)
pool2    = MaxPoolingLayer(:pool2, net, relu2, 2, 2, 0)
conv3    = ConvolutionLayer(:conv3, net, pool1, 50, 3, 1, 1)
relu3    = ReLULayer(:relu3, net, conv3)
pool3    = MaxPoolingLayer(:pool3, net, relu3, 2, 2, 0)
fc4      = InnerProductLayer(:fc4, net, pool3, 512)
relu4    = ReLULayer(:relu4, net, fc4)
drop4    = DropoutLayer(:drop4, net, relu4, .5f0)
fc5      = InnerProductLayer(:fc5, net, drop4, 512)
relu5    = ReLULayer(:relu5, net, fc5)
drop5    = DropoutLayer(:drop5, net, relu5, .5f0)
fc6      = InnerProductLayer(:fc6, net, drop5, 10)
loss     = SoftmaxLossLayer(:loss, net, fc6, label)
accuracy = AccuracyLayer(:accuracy, net, fc6, label)

params = SolverParameters(
    LRPolicy.Inv(0.01, 0.0001, 0.75),
    MomPolicy.Fixed(0.9),
    10000,
    .0005,
    10000)
sgd = SGD(params)
solve(sgd, net)


text:
type PopulationDumpLogger <: Logger
  # The destination directory of the population dumps.
  destination::AbstractString

  PopulationDumpLogger() = new()
end

"""
The population dump logger writes the state of each member of the population
to a single log file, including their fitness and their representation(s), at
the end of each generation.
"""
population_dump() = PopulationDumpLogger()

function prepare!(l::PopulationDumpLogger, output::AbstractString)
  l.destination = joinpath(output, "population_dumps")
  mkdir(l.destination)
end

"""
Creates a new dump file for the current generation within the destination
folder specified by the logger.
"""
#function log!(l::PopulationDumpLogger, s::State)
# f = open(joinpath(l.destination, string(s.iterations)), "w")
# buffer = ""
#
# for (i, deme) in enumerate(s.population.demes)
#  buffer *= "Deme $(i):\n\n"
#  for ind in deme.members
#    buffer *= describe(ind) * "\n\n"
#  end
# end
#
# write(f, buffer)
# close(f)
#end


text:
module GenericRefFEs

using Gridap
using Gridap.Helpers

using LinearAlgebra

export GenericRefFE
export RTRefFE
export NedelecRefFE
export GenericDOFBasis

import Gridap: shfbasis
import Gridap: polytope
import Gridap: nfacedofs
import Gridap: dofbasis
import Gridap: evaluate!

import Base:length

# """
# It represents the DOF basis for a generic reference FE, in which we provide the vector of nodes
# in which a function must be evaluated in order to be able to compute a set of
# moments. The array of moments (functionals) take the vector of nodal values of
# the function and returns the evaluation of moments. Thus, we can define the
# set of nodes as $\{ xn_i \}$, and the set of moments as $m_i = \sum m_ij f(xn_j)$.
# We note that $m_ij$ is a field value, i.e., if $f$ is a scalar (resp., vector,
# or tensor) field then $m_ij$ is a scalar, (resp., vector, or tensor) value.
#
# We note that the nodes and moments arrays are stored n-face wise. Nodes and
# moments in reference FEs are owned by a particular n-face. By imposing continuity
# of these moments (degrees of freedom) across cells in the physical space we can
# prove a particular type of continuity (e.g., normal, tangent, or full continuity
# for div, curl, or grad-conforming FE spaces). Since the nodal values required
# for every n-face are different, to segregate these computations n-face-wise does
# not introduce any overhead. In fact, it is more effective for the storage of
# $m_ij$, since one does not have to store the zeros related to nodes of other
# n-faces.
# """
struct GenericDOFBasis{D,T,S} <: DOFBasis{D,T}
  nodes::Vector{Array{Point{D,S}}}
  moments::Vector{Array{T}}
  _cache_field#::Vector{T}
  _cache_basis#::Matrix{T}
end

nodes(b::GenericDOFBasis) = b.nodes

moments(b::GenericDOFBasis) = b.moments

length(b::GenericDOFBasis{D,T} where {D,T}) = sum([size(i,1) for i in b.moments])

# """
# It represents the reference FE related to the `GenericDOFBasis` (see
# documentation for more details). The rest of ingredients are standard. It
# involves a polytope, a basis generated as a prebasis and change-of-basis to
# get the canonical basis related to the DOF basis, and an array of nfacedofs
# that provide the ids of dofs for a given n-face.
# """
# @santiagobadia : Do we really need nfacedofs?
# @santiagobadia : We could probably put polytope reference in dofbasis since
# the dofbasis is n-face wise
struct GenericRefFE{D,T} <: RefFE{D,T}
  polytope::Polytope{D}
  dof_basis::GenericDOFBasis{D,T}
  shfbasis::Basis{D,T}
  nfacedofs::Vector{Vector{Int}}
end

dofbasis(this::GenericRefFE{D,T} where {D,T})::DOFBasis = this.dof_basis

polytope(this::GenericRefFE{D,T} where {D,T})::Polytope = this.polytope

shfbasis(this::GenericRefFE{D,T} where {D,T})::Basis = this.shfbasis

nfacedofs(this::GenericRefFE{D,T} where {D,T})::Vector{Vector{Int}} = this.nfacedofs

function _GenericRefFE(p::Polytope{D}, dof_basis::GenericDOFBasis{D,T,S}, shfbasis::Basis{D,T}, nfacedofs) where {D,T,S}
  GenericRefFE{D,T}(p,dof_basis,shfbasis,nfacedofs)
end

function GenericDOFBasis(nodes::Vector{Array{Point{D,S}}}, moments::Vector{Array{T}}) where {D,T,S}
  ndofs = sum([size(i,1) for i in moments])
  nnodes = [length(i) for i in nodes]
  cache_field = [ zeros(T,i) for i in nnodes]
  cache_basis = [ zeros(T,ndofs,i) for i in nnodes]

  GenericDOFBasis{D,T,S}(
  nodes,
  moments,
  cache_field,
  cache_basis)
end

function evaluate!(
  b::GenericDOFBasis{D,T},f::Field{D,T},dofs::AbstractVector{E}) where {D,T,E}
  for (n,v) in zip(b.nodes,b._cache_field)
    evaluate!(f,n,v)
  end
  k = 0
  for (m,v) in zip(b.moments,b._cache_field)
    l = size(m,1)
    if length(m) > 0
      dofs[k+1:k+l] = m*v
      k += l
    end
  end
  return dofs
end

function evaluate!(
  b::GenericDOFBasis{D,T},f::Basis{D,T},dofs::AbstractMatrix{E}) where {D,T,E}
  for (n,v) in zip(b.nodes,b._cache_basis)
    evaluate!(f,n,v)
  end
  k = 0
  for (m,v) in zip(b.moments,b._cache_basis)
    l = size(m,1)
    if length(m) > 0
      dofs[:,k+1:k+l] = v*m'
      k += l
    end
  end
  dofs
end

function _nfaces_evaluation_points_weights(p, fgeomap, fips, wips)
  nc = length(fgeomap)
  c_fips = ConstantCellValue(fips,nc)
  c_wips = ConstantCellValue(wips,nc)
  pquad = evaluate(fgeomap,c_fips)
  c_fips, pquad, c_wips
end

# Ref FE to faces geomaps
function _ref_face_to_faces_geomap(p,fp)
  cfvs = nfaces_vertices(p,dim(fp))
  nc = length(cfvs)
  freffe = LagrangianRefFE(Float64,fp,1)
  fshfs = shfbasis(freffe)
  cfshfs = ConstantCellValue(fshfs, nc)
  fgeomap = lincomb(cfshfs,cfvs)
end

function _broadcast(::Type{T},n,b) where T
  c = Array{T}(undef,size(b))
  for (ii, i) in enumerate(b)
    c[ii] = i*n
  end
  return c
end

function _broadcast_cross(::Type{T},n,b) where T
  c = Array{T}(undef,size(b))
  for (ii, i) in enumerate(b)
    c[ii] = T(cross(i.array,n.array))# cross product
  end
  return c
end

function _broadcast_extend(::Type{T},Tm,b) where T
  c = Array{T}(undef,size(b))
  for (ii,i) in enumerate(b)
    c[ii] = T(Tm*[i...])
  end
  return c
end

function _nfacedofs_basis(p,moments)
  ndofs = [size(moments[i],1) for i in 1:length(moments)]
  _nfacedofs = Vector{Vector{Int}}(undef,num_nfaces(p))
  _nfacedofs[1] = Int[1:ndofs[1]...]
  k = ndofs[1]
  for i in 2:length(ndofs)
    _nfacedofs[i] = [k+1:k+ndofs[i]...]
    k += ndofs[i]
  end
  return _nfacedofs
end

function _nfaces_array_dim!(p,dim,array,nf_vals)
  nfs = nfaces_dim(p,dim)
  array[nfs] = nf_vals
end

function _initialize_arrays(prebasis,p)

  # Field, point, and entry types
  ft = eltype(prebasis)
  et = eltype(ft)
  pt = Point{dim(p),Float64}

  # Arrays of moments (as its evaluation for all prebasis shape functions)
  # and evaluation points per n-face
  nf_moments = Vector{Array{ft}}(undef,num_nfaces(p))
  nf_nodes = Vector{Array{pt}}(undef,num_nfaces(p))
  nshfs = length(prebasis)
  pb_moments = zeros(et,nshfs,0)

  # Initialize to zero arrays
  zero_moments = zeros(ft,0,0)
  zero_nodes = zeros(pt,0)
  for dim in 0:dim(p)
    for inf in nfaces_dim(p,dim)
      nf_moments[inf] = zero_moments
      nf_nodes[inf] = zero_nodes
    end
  end
  return nf_nodes, nf_moments, pb_moments

end

function _GenericRefFE(p::Polytope,prebasis::Basis,
                      nf_nodes,nf_moments,pb_moments)

  # Change of basis matrix, inv([DF,DC])
  cob = inv(hcat(pb_moments))
  basis = change_basis(prebasis,cob)

  nfacedofs = _nfacedofs_basis(p,nf_moments)

  # Build DOFBasis and RefFE with all this
  dof_basis = GenericDOFBasis(nf_nodes, nf_moments)

  divreffe = _GenericRefFE(p,dof_basis,basis,nfacedofs)

end

function _insert_nface_values!(nf_nodes, nf_moments, pb_moments,
                               prebasis, fcips, fmoments, p, dim)

  # Evaluate basis in faces points, i.e., S(Fi)_{ab} = ϕ^a(xgp_Fi^b)
  pbasis_fcips = [evaluate(prebasis,ps) for ps in fcips]

  # Face moments evaluated for basis, i.e., DF = [S(F1)*M(F1)^T, …, S(Fn)*M(Fn)^T]
  fms_preb = [bps*ms' for (bps,ms) in zip(pbasis_fcips,fmoments)]

  _nfaces_array_dim!(p,dim,nf_moments,fmoments)
  _nfaces_array_dim!(p,dim,nf_nodes,fcips)
  pb_moments = hcat(pb_moments,fms_preb...)

  return nf_nodes, nf_moments, pb_moments

end

include("DivRefFEs.jl")

include("CurlRefFEs.jl")

end # module


text:
# Autogenerated wrapper script for VideoCaptureWrap_jll for x86_64-apple-darwin14-cxx11
export libvideocapture

using libcxxwrap_julia_jll
using OpenCVQt_jll
## Global variables
PATH = ""
LIBPATH = ""
LIBPATH_env = "DYLD_FALLBACK_LIBRARY_PATH"
LIBPATH_default = "~/lib:/usr/local/lib:/lib:/usr/lib"

# Relative path to `libvideocapture`
const libvideocapture_splitpath = ["lib", "libvideocapture.dylib"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libvideocapture_path = ""

# libvideocapture-specific global declaration
# This will be filled out by __init__()
libvideocapture_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libvideocapture = "@rpath/libvideocapture.dylib"

# Inform that the wrapper is available for this platform
wrapper_available = true

"""
Open all libraries
"""
function __init__()
    # This either calls `@artifact_str()`, or returns a constant string if we're overridden.
    global artifact_dir = find_artifact_dir()

    global PATH_list, LIBPATH_list
    # Initialize PATH and LIBPATH environment variable listings
    # From the list of our dependencies, generate a tuple of all the PATH and LIBPATH lists,
    # then append them to our own.
    foreach(p -> append!(PATH_list, p), (libcxxwrap_julia_jll.PATH_list, OpenCVQt_jll.PATH_list,))
    foreach(p -> append!(LIBPATH_list, p), (libcxxwrap_julia_jll.LIBPATH_list, OpenCVQt_jll.LIBPATH_list,))

    global libvideocapture_path = normpath(joinpath(artifact_dir, libvideocapture_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libvideocapture_handle = dlopen(libvideocapture_path, RTLD_LAZY | RTLD_DEEPBIND)
    push!(LIBPATH_list, dirname(libvideocapture_path))

    # Filter out duplicate and empty entries in our PATH and LIBPATH entries
    filter!(!isempty, unique!(PATH_list))
    filter!(!isempty, unique!(LIBPATH_list))
    global PATH = join(PATH_list, ':')
    global LIBPATH = join(vcat(LIBPATH_list, [joinpath(Sys.BINDIR, Base.LIBDIR, "julia"), joinpath(Sys.BINDIR, Base.LIBDIR)]), ':')

end  # __init__()


text:
# Kernel Type
abstract type AbstractKernelType end
struct ReadOnly <: AbstractKernelType end
struct WriteOnly <: AbstractKernelType end
struct ReadWrite <: AbstractKernelType end

# LoadStoreType
abstract type AbstractLoadStoreType end
struct Standard <: AbstractLoadStoreType end
struct Nontemporal <: AbstractLoadStoreType end

# Iteration Style
abstract type AbstractIterationStyle end
struct Sequential <: AbstractIterationStyle end
struct PseudoRandom <: AbstractIterationStyle end

"""
    KernelParam{K,L,S,I,T,N,U}

Top level type governing the creation of a microbenchmark kernel.
Construct using Keywords:

    KernelParam(; [kw...])

The available keywords are:

* `kernel`: Controls the read/write behavior of the kernel. Options are:
    - `ReadOnly`: Only generate loads.
    - `WriteOnly`: Only generate stores.
    - `ReadWrite`: Generate a load followed by a store.
    Default: `ReadOnly`

* `loadtype`: Controls whether standard or nontemporal loads are emitted. Options are:
    - `Standard`: Use standard loads.
    - `Nontemporal`: Use nontemporal loads.
    Default: `Standard`

* `storetype`: Controls whether standard or nontemporal stores are emitter. Options are:
    - `Standard`: Use standard stores.
    - `Nontemporal`: Use nontemporal stores.
    Default: `Standard`

* `iterator`: Controls how the underlying array is accessed. Optiosn are:
    - `Sequential`: Touch every index sequentially.
    - `PseudoRandom`: Use a MaxLFSR to touch each index in a pseudo-random order.
    Default: `Sequential`

* `eltype`: The element type of the underlying array. Default: `Float32`.

* `vectorsize`: The size of the vector to use. This determines whether 128 bit, 256 bit,
    or 512 bit vector instructions are used. The default is `16` to enable AVX512
    instructions for `Float32` entries.

* `unroll`: The number of times the innermost load/store operation is performed.
    Especially useful for the pseudo-random benchmarks to get read/write granularity larger
    than 64 B. Default: `1`.
"""
struct KernelParam{
        # Kernel Type
        K <: AbstractKernelType,
        # Load Type
        L <: AbstractLoadStoreType,
        # Store Type
        S <: AbstractLoadStoreType,
        # Iteration Style
        I <: AbstractIterationStyle,
        # Vector element Type
        T,
        # Vector Width (bytes)
        N,
        # Unroll factor
        U,
        # Loop over unrolled code
        O,
    }

    function KernelParam(;
            kernel = ReadOnly,
            loadtype = Standard,
            storetype = Standard,
            iterator = Sequential,
            eltype = Float32,
            vectorsize = 8,
            unroll = 1,
            loops = 1,
        )

        return new{kernel,loadtype,storetype,iterator,eltype,vectorsize,unroll,loops}()
    end
end

# Convenience accessors
kernel(::Type{KernelParam{K,L,S,I,T,N,U,O}}) where {K,L,S,I,T,N,U,O} = K
kernel(K::KernelParam) = kernel(typeof(K))

loadtype(::Type{KernelParam{K,L,S,I,T,N,U,O}}) where {K,L,S,I,T,N,U,O} = L
loadtype(K::KernelParam) = loadtype(typeof(K))

storetype(::Type{KernelParam{K,L,S,I,T,N,U,O}}) where {K,L,S,I,T,N,U,O} = S
storetype(K::KernelParam) = storetype(typeof(K))

iterator(::Type{KernelParam{K,L,S,I,T,N,U,O}}) where {K,L,S,I,T,N,U,O} = I
iterator(K::KernelParam) = iterator(typeof(K))

Base.eltype(::Type{KernelParam{K,L,S,I,T,N,U,O}}) where {K,L,S,I,T,N,U,O} = T
Base.eltype(K::KernelParam) = eltype(typeof(K))

vectorsize(::Type{KernelParam{K,L,S,I,T,N,U,O}}) where {K,L,S,I,T,N,U,O} = N
vectorsize(K::KernelParam) = vectorsize(typeof(K))

unroll(::Type{KernelParam{K,L,S,I,T,N,U,O}}) where {K,L,S,I,T,N,U,O} = U
unroll(K::KernelParam) = unroll(typeof(K))

loops(::Type{KernelParam{K,L,S,I,T,N,U,O}}) where {K,L,S,I,T,N,U,O} = O
loops(K::KernelParam) = loops(typeof(K))

#####
##### Upper level handling
#####

# The upper level functions reinterpret the top level array as an array of `Vec` elements.
# Performs alignment checking and length checking.
# Constructs the iterator
"""
    execute!(A::AbstractVector, K::KernelParam)

Run the kernel described by `K` on `A`.
The following restrictions apply to `K`:

* `eltype(A) == eltype(K)`
* `A` must be aligned to 64 bytes (i.e., aligned to a cache line)
* The lengthh of `A` must be a multiple of `unroll * vectorsize * Treads.nthreads()`
"""
function execute!(A::AbstractVector, K::KernelParam)
    # Sanity check
    @assert eltype(A) == eltype(K)
    # Alignment check
    @assert iszero(mod(convert(Int, pointer(A)), 64))
    # Size check
    @assert iszero(mod(length(A), unroll(K) * vectorsize(K)))

    V = reinterpret(K, A)
    itr = makeitr(V, K)

    # Inner high-performance function.
    return _execute!(V, itr, K)
end

Base.reinterpret(K::KernelParam, A::AbstractArray) = reinterpret(Vec{vectorsize(K),eltype(K)}, A)

# iterators
makeitr(A, K::KernelParam) = makeitr(A, iterator(K), unroll(K), loops(K))

makeitr(A, ::Type{Sequential}, unroll, loops) = 1:div(length(A), unroll * loops)
makeitr(A, ::Type{PseudoRandom}, unroll, loops) = MaxLFSR.LFSR(div(length(A), unroll * loops))

#####
##### Generate the kernel
#####

# Note: This is not exactly an extendable API since we're relying on `@generated` functions.
# In particular, in order for any changes made to the `hbf` function or any functions
# called within to take effect, we need to restart Julia.
#
# This is a limitation of the dependency tracking in Julia.
@generated function _execute!(A::AbstractArray{Vec{N,T}}, itr, K::KernelParam) where {N,T}
    return emit(A, K)
end

function emit(A::Type{<:AbstractArray{Vec{N,T}}}, K::Type{<:KernelParam}) where {N,T}
    @nospecialize

    # Get the header, body, and footer for this kernel
    header, body, footer = hbf(A, K)
    nloops = loops(K)

    # Emit the rest of the function
    return quote
        $(header...)
        base = Base.unsafe_convert(Ptr{$T}, pointer(A))
        @inbounds for i in itr
            for j in 1:$(nloops)
                ptr = base + $(unroll(K) * sizeof(Vec{N,T})) * ($nloops * (i-1) + (j-1))
                $(body...)
            end
        end
        $(footer...)
    end
end

# header, body, footer
function hbf(::Type{<:AbstractArray{Vec{N,T}}}, K::Type{<:KernelParam}) where {N,T}
    @nospecialize

    # Check invariants
    @assert N == vectorsize(K)
    @assert T == eltype(K)

    lt = loadtype(K)
    st = storetype(K)

    # defaults
    header = Any[]
    footer = Any[:(return nothing)]

    # Check if we need to emit loads/stores
    if kernel(K) == ReadOnly
        # Header
        header = push!(header, :(s = zero(eltype(A))))

        # Body
        body = emit_loads(Vec{N,T}, unroll(K), loadtype(K))
        reduction = unroll(K) == 1 ? :(s += $(symbols(1)[1])) : :(s += +($(symbols(unroll(K))...)))
        push!(body, reduction)

        # footer
        footer = Any[:(return s)]
    elseif kernel(K) == WriteOnly
        # body
        body = emit_stores(Vec{N,T}, unroll(K), storetype(K), false)
    elseif kernel(K) == ReadWrite
        # body
        loads = emit_loads(Vec{N,T}, unroll(K), loadtype(K))
        stores = emit_stores(Vec{N,T}, unroll(K), storetype(K), true)
        body = vcat(loads, stores)
    else
        error("Unknown Kernel: $(kernel(K))")
    end

    return header, body, footer
end

symbols(N) = [Symbol("v$i") for i in 1:N]
lower(::Type{Nontemporal}) = Val{true}()
lower(::Type{Standard}) = Val{false}()

function emit_loads(vectype::Type{<:Vec}, unroll, loadtype)
    @nospecialize

    syms = symbols(unroll)
    return map(1:unroll) do i
        lhs = syms[i]
        shift = sizeof(vectype) * (i-1)
        return :($lhs = vload($vectype, ptr + $shift, nothing, Val{true}(), $(lower(loadtype))))
    end
end

function emit_stores(vectype::Type{<:Vec}, unroll, storetype, follows_load::Bool)
    @nospecialize

    syms = symbols(unroll)
    o = one(vectype)
    return map(1:unroll) do i
        rhs = syms[i]
        shift = sizeof(vectype) * (i-1)
        vec = follows_load ? :($rhs + $o) : o
        return :(vstore($vec, ptr + $shift, nothing, Val{true}(), $(lower(storetype))))
    end
end


text:
function Base.reduce(f, a::IndicesArray; dims=:, kwargs...)
    d = to_dims(a, dims)
    return rebuild(a, reduce(f, parent(a); dims=d, kwargs...), reduce_axes(a, d))
end

function Base.mapreduce(f, op, a::IndicesArray; dims=:, kwargs...)
    d = to_dims(a, dims)
    return rebuild(a, mapreduce(f, op, parent(a); dims=d, kwargs...), reduce_axes(a, d))
end

function Base.mapslices(f, a::IndicesArray; dims, kwargs...)
    d = to_dims(a, dims)
    return rebuild(a, mapslices(f, parent(a); dims=d, kwargs...), reduce_axes(a, d))
end


text:
include("storage.jl")
include("chunk.jl")
include("schunk.jl")
include("unsafe_functions.jl")
include("functions.jl")


text:
using VortexDistributions
using JLD2
using GLMakie
using Interpolations
using NearestNeighbors
using LinearAlgebra
using Distances
## Include util functions
include("../util.jl")

## Load psi simulation
@load "data/box_rings.jld2" psi_ring1 psi_ring2 psi_ring3 X # make sure julia repl is in base directory
# @load "data/box_vorts.jld2" psi_chaos1 psi_chaos2 psi_knots1 psi_knots2 psi_knots3 psi_knots4 psi_tubes1 psi_tubes2 psi_ringtube X
# @load "data/slab_tubes.jld2" psi_tubes1 psi_tubes2 psi_tubes3 X

## Set psi
psi = psi_ring3;

## Plot iso surface
plot_iso(psi, X)

naivePlaquette(psi[:, :, 42], X, 0)
findvortices(Torus(psi[:, :, 42], X[1], X[2]))

## Find vortices
@time vorts_3d = findvortices3D_itp(psi, X, 2);
num_vorts = length(vcat(vorts_3d...)[:,1])
print("Number of vortices found: ", length(vcat(vorts_3d...)[:,1]))

## Scatter vortices found
scatterVortsOnIso(vcat(vorts_3d'...))

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

## Unique vortices
interp_depth = 2
rtol = 0.44 # Largest distance between two points

@time fils = uniqueVortices(rtol, psi, X, interp_depth);

length(fils) # number of filaments found

## Plotting unique filaments
plot_iso(psi, X)
scatterDemo(fils)

plot!()


text:
load("cached.jl")

# code = :(function f(c::Context, x::Int)
#     x^2
# end)

code = :(f(c::Context, x::Int) = x^2)

c = wrap_cached(code)
print(c)

@cached function cr(c::Context, x)
    randi(100)
end
@cached function cr2(c::Context, x)
    randi(100)
end

c = Cache()

println()
@show cr(c,1), cr2(c,1)
@show cr(c,2), cr2(c,2)
@show cr(c,1), cr2(c,1)

println()
@show cr(c,(2,3))
@show cr(c,"hej")

max_stars_repo_name: toivoh/julia-kernels

max_issues_repo_name: toivoh/julia-kernels

max_forks_repo_name: toivoh/julia-kernels


text:
type TextureAtlas
    rectangle_packer::RectanglePacker
    mapping         ::Dict{Any, Int} # styled glyph to index in sprite_attributes
    index           ::Int
    images          ::Texture{Float16, 2}
    attributes      ::Vector{Vec4f0}
    scale           ::Vector{Vec2f0}
    extent          ::Vector{FontExtent{Float64}}

    function TextureAtlas(initial_size=(4096, 4096))

        images = Texture(fill(Float16(0.0), initial_size...), minfilter=:linear, magfilter=:linear)

        new(
            RectanglePacker(SimpleRectangle(0, 0, initial_size...)),
            Dict{Any, Int}(),
            1,
            images,
            Vec4f0[],
            Vec2f0[],
            FontExtent{Float64}[]
        )
    end
end

begin #basically a singleton for the textureatlas
const local TEXTURE_ATLAS = TextureAtlas[]
reset_texture_atlas!() = empty!(TEXTURE_ATLAS)
function get_texture_atlas()
    not_initilized = isempty(TEXTURE_ATLAS)
    if not_initilized
        global DEFAULT_FONT_FACE = newface(assetpath("fonts", "DejaVuSans.ttf"))
        alternatives = [
            "DejaVuSans.ttf",
            "NotoSansCJKkr-Regular.otf",
            "NotoSansCuneiform-Regular.ttf",
            "NotoSansSymbols-Regular.ttf",
            "FiraMono-Medium.ttf"
        ]
        global ALTERNATIVE_FONTS = map(alternatives) do font
            newface(assetpath("fonts", font))
        end
        atlas = push!(TEXTURE_ATLAS, TextureAtlas())[] # initialize only on demand
        for c in '\u0000':'\u00ff' #make sure all ascii is mapped linearly
            insert_glyph!(atlas, c, DEFAULT_FONT_FACE)
        end
        return atlas
    else
        return TEXTURE_ATLAS[]
    end
end
end

function glyph_index!(atlas::TextureAtlas, c::Char, font)
    if FT_Get_Char_Index(font[], c) == 0
        for afont in ALTERNATIVE_FONTS
            if FT_Get_Char_Index(afont[], c) != 0
                font = afont
            end
        end
    end
    if c < '\u00ff' && font == DEFAULT_FONT_FACE # characters up to '\u00ff'(255), are directly mapped for default font
        Int(c)+1
    else #others must be looked up, since they're inserted when used first
        return insert_glyph!(atlas, c, font)
    end
end

glyph_scale!(c::Char) = glyph_scale!(get_texture_atlas(), c, DEFAULT_FONT_FACE)
glyph_uv_width!(c::Char) = glyph_uv_width!(get_texture_atlas(), c, DEFAULT_FONT_FACE)

function glyph_uv_width!(atlas::TextureAtlas, c::Char, font)
    atlas.attributes[glyph_index!(atlas, c, font)]
end
function glyph_scale!(atlas::TextureAtlas, c::Char, font)
    atlas.scale[glyph_index!(atlas, c, font)]
end
function glyph_extent!(atlas::TextureAtlas, c::Char, font)
    atlas.extent[glyph_index!(atlas, c, font)]
end

function bearing(extent)
     Point2f0(extent.horizontal_bearing[1], -(extent.scale[2]-extent.horizontal_bearing[2]))
end
function glyph_bearing!{T}(atlas::TextureAtlas, c::Char, font, scale::T)
    T(bearing(atlas.extent[glyph_index!(atlas, c, font)])) .* scale
end
function glyph_advance!{T}(atlas::TextureAtlas, c::Char, font, scale::T)
    T(atlas.extent[glyph_index!(atlas, c, font)].advance) .* scale
end

insert_glyph!(atlas::TextureAtlas, glyph::Char, font) = get!(atlas.mapping, (glyph, font)) do
    uv, rect, extent, width_nopadd = render(atlas, glyph, font)
    tex_size       = Vec2f0(size(atlas.images))
    uv_start       = Vec2f0(uv.x, uv.y)
    uv_width       = Vec2f0(uv.w, uv.h)
    real_heightpx  = width_nopadd[2]
    halfpadding    = (uv_width - width_nopadd) / 2f0
    real_start     = uv_start + halfpadding # include padding
    relative_start = real_start ./ tex_size # use normalized texture coordinates
    relative_width = (real_start+width_nopadd) ./ tex_size

    uv_offset_width = Vec4f0(relative_start..., relative_width...)
    i               = atlas.index
    push!(atlas.attributes, uv_offset_width)
    push!(atlas.scale, Vec2f0(width_nopadd))
    push!(atlas.extent, extent)
    atlas.index = i+1
    return i
end

function sdistancefield(img, restrict_steps=2)
    w, h = size(img)
    w1, h1 = w, h

    halfpad = 24*(2^restrict_steps) # padd so that after restrict it comes out as roughly 48 pixel
    w, h = w+2halfpad, h+2halfpad #pad this, to avoid cuttoffs
    in_or_out = Bool[begin
        x, y = i-halfpad, j-halfpad
        if checkbounds(Bool, img, x,y)
            img[x,y] >= 1.0
        else
            false
        end
    end for i=1:w, j=1:h]

    sd = sdf(in_or_out)
    for i=1:restrict_steps
        w1, h1 = Images.restrict_size(w1), Images.restrict_size(h1)
        sd = Images.restrict(sd) #downsample
    end
    sz = Vec2f0(size(img))
    maxlen = norm(sz)
    sw, sh = size(sd)

    Float16[clamp(sd[i,j]/maxlen, -1, 1) for i=1:sw, j=1:sh], Vec2f0(w1, h1), (2^restrict_steps)
end

function GLAbstraction.render(atlas::TextureAtlas, glyph::Char, font)
    #select_font_face(cc, font)
    if glyph == '\n' # don't render  newline
        glyph = ' '
    end
    bitmap, extent = renderface(font, glyph, (164, 164))
    restrict_steps=2
    sd, width_nopadd, scaling_factor = sdistancefield(bitmap, restrict_steps)
    extent = extent ./ Vec2f0(2^restrict_steps)
    rect = SimpleRectangle(0, 0, size(sd)...)
    uv   = push!(atlas.rectangle_packer, rect) #find out where to place the rectangle
    uv == nothing && error("texture atlas is too small. Resizing not implemented yet. Please file an issue at GLVisualize if you encounter this") #TODO resize surface
    atlas.images[uv.area] = sd
    uv.area, rect, extent, width_nopadd
end


text:
##

"""
    $SIGNATURES

Hacky version to return which factor type to use between two variables of types T1 and T2.

Pending work on RoME.jl #244 and IIF.jl #1010 and others, after which this will be refactored.
"""
selectFactorType(T1::InstanceType{ContinuousScalar}, T2::InstanceType{ContinuousScalar}) = LinearRelative{1}
selectFactorType(T1::Type{<:ContinuousEuclid{N}}, T2::Type{<:ContinuousEuclid{N}}) where N = LinearRelative{N}
selectFactorType(T1::InferenceVariable, T2::InferenceVariable) = selectFactorType(typeof(T1), typeof(T2))
selectFactorType(dfg::AbstractDFG, s1::Symbol, s2::Symbol) = selectFactorType( getVariableType(dfg, s1), getVariableType(dfg, s2) )


text:
module MetaLearning

end # module


text:
using Parameters
include("common.jl")

# Creatine kinase and creatine shuttle parameters
@with_kw struct CKParams
    KCK_IN= 1.4E-4
    KCK_MT = 1.33E-6
    KTR_CR = 2E-3
    KEQ_CR = 0.0095
    # Basal consumption of cytosolic ATP (mM/ms)
    V_ATPASE_CYTO = 1E-5
    # Pool of creatine and creatine phosphate (mM)
    CR_T = 25.0
end

# CK shuttle rates
function ck_system(atp_i, atp_ic, crp_i, crp_ic, V_ATPASE_CYTO, KTR_CR, KCK_IN,
                  KCK_MT, KEQ, CR_T, AXP_I_T)
    cr_ic = CR_T - crp_ic
    cr_i = CR_T - crp_i
    adp_ic = AXP_I_T - atp_ic
    adp_i = AXP_I_T - atp_i

    vck_cyto = KCK_IN * (atp_ic * cr_ic - adp_ic * crp_ic / KEQ)
    vck_mito = KCK_MT * (atp_i * cr_i - adp_i * crp_i / KEQ)
    vtr_crp = KTR_CR * (crp_i - crp_ic)

    d_atp_in_ck = -vck_mito
    d_atp_ic = -vck_cyto - V_ATPASE_CYTO
    d_crp_i = vck_mito - vtr_crp
    d_crp_ic = vtr_crp + vck_cyto
    return (d_atp_in_ck, d_atp_ic, d_crp_i, d_crp_ic)
end

function ck_system(atp_i, atp_ic, crp_i, crp_ic, pCK::CKParams, AXP_I_T)
    @unpack V_ATPASE_CYTO, KTR_CR, KCK_IN, KCK_MT, KEQ_CR, CR_T = pCK
    ck_system(atp_i, atp_ic, crp_i, crp_ic, V_ATPASE_CYTO, KTR_CR, KCK_IN, KCK_MT, KEQ_CR, CR_T, AXP_I_T)
end


text:
## Do inference in parallel

using Distributed

## add number of processors desired

addprocs(4)

## Load relevant code into scope for each worker

@everywhere begin
    using Pkg;Pkg.activate(".")
end
@everywhere using JLD2,LinearAlgebra,TransformVariables,LogDensityProblems,MCMCChains,Dates
@everywhere import KenyaCoVSD

## Load data

@load("data/linelist_data_with_pos_neg_20feb_to_27apr_c_age.jld2")
@load("data/serological_data_with_20feb_to_31dec_age.jld2")
@load("data/deaths_data20210521.jld2")
@load("data/N_kenya.jld2")

## Parameter transformations and initial KE guess for HMC

trans = as((ct_min1 = as(Real, 0.4, 0.999),
        R₀ = as(Real, 0.0, 7.0),
        ϵ = as(Real, 0.7, 0.999),
        χ₁=as(Real, 0.1, 30.0),χ₂=as(Real, 0.1, 30.0),
        p_test₁=as(Real, 5e-6, 5e-3),p_test₂ = as(Real, 5e-6, 5e-3),
        p_choose1=as(Real, 0.1, 0.999),
        P_eff=as(Real, 0.5, 0.999),
        schooleffect = as(Real, 0.5, 0.999),
        extra_transmissibility = as(Real, 0.5, 2.5),
        influx_exposed_new_variant = as(Real, 0.0, 10e3),
        p_test_boost = as(Real, 1.0, 3.0),
        χ_boost= as(Real, 1.0, 3.0),
        E₀=as(Real, 0.0, 10e3)))

D = Diagonal(0.1*ones(15))

## Name for each county with substantial serological data, and matching priors

counties_with_serodata = ["Nairobi","Mombasa",
                        "Nakuru","Uasin Gishu","Embu","Kisumu","Siaya","Kisii","Nyeri",
                            "Kilifi","Kwale"]

priors = [fill(KenyaCoVSD.priors_twogroup_newvariant_cities,2);
           fill(KenyaCoVSD.priors_twogroup_newvariant_semi_urbanrural,7);
            fill(KenyaCoVSD.priors_twogroup_newvariant_rural,2)]

## parallelised loop for fitting each county
@distributed for i = 1:11
    name = counties_with_serodata[i]
    prior = priors[i]
    model = KenyaCoVSD.CoVAreaModel(name,KenyaCoVSD.ll_twogroup_newvariant_infboost,prior;
                case_data = linelist_data_with_pos_neg,
                sero_data = serological_data,
                death_data = deaths_data,
                pop_data = N_kenya)

    KenyaCoVSD.inferparameters!(model,2000,trans,0.05,D;serowaningrate = 0.0)

    @save("modelfits/$(name)_model.jld2",model)
end


text:
# This file contains functions for calculating spatial lags

"""
    slag(W, x)
Calculate the spatial lag of `x` using the spatial weights `W`.
"""
function slag(W::SpatialWeights, x::AbstractVector{T}) where {T <:Real}
    n = W.n
    nx = length(x)

    n == nx || throw(DimensionMismatch("dimensions must match: W has ($(n)), x has ($nx)"))

    sx = zeros(n)
    @inbounds @simd for i in 1:n
        ni = length(W.neighs[i])
        for j in 1:ni
            sx[i] = sx[i] + W.weights[i][j] * x[W.neighs[i][j]]
        end
    end

    return sx

end

Base.:*(W::SpatialWeights, x::AbstractVector{T}) where {T <:Real} = slag(W, x);


text:
module Backup

include("module_common.jl")

const _libcloud_backup_types = PyCall.PyNULL()
const _libcloud_backup_providers = PyCall.PyNULL()
const _libcloud_backup_base = PyCall.PyNULL()

immutable BackupDriver
    o::PyObject

    function BackupDriver(provider::Compat.String, args...; kwargs...)
        cls = _libcloud_backup_providers[:get_driver](provider)
        BackupDriver(cls(args...; kwargs...))
    end
    BackupDriver(driver::PyObject) = new(driver)
end
show(io::IO, c::BackupDriver) = print(io, c.o[:__str__]())

immutable BackupTarget
    o::PyObject

    id::Nullable{Compat.String}
    name::Compat.String
    address::Compat.String
    typ::Compat.String
    extra::Dict
    driver::BackupDriver

    function BackupTarget(o::PyObject)
        new(o, o[:id], o[:name], o[:address], o[:type], o[:extra], BackupDriver(o[:driver]))
    end
end
PyObject(o::BackupTarget) = o.o
convert(::Type{BackupTarget}, o::PyObject) =  BackupTarget(o)
show(io::IO, o::BackupTarget) = print(io, o.o[:__str__]())
@doc LazyHelp(_libcloud_backup_base, "BackupTarget") BackupTarget

immutable BackupTargetJob
    o::PyObject

    id::Nullable{Compat.String}
    status::Compat.String
    progress::Int
    target::BackupTarget
    extra::Dict
    driver::BackupDriver

    function BackupTargetJob(o::PyObject)
        new(o, o[:id], o[:status], o[:progress], BackupTarget(o[:target]), o[:extra], BackupDriver(o[:driver]))
    end
end
PyObject(o::BackupTargetJob) = o.o
convert(::Type{BackupTargetJob}, o::PyObject) =  BackupTargetJob(o)
show(io::IO, o::BackupTargetJob) = print(io, o.o[:__str__]())
@doc LazyHelp(_libcloud_backup_base, "BackupTargetJob") BackupTargetJob

immutable BackupTargetRecoveryPoint
    o::PyObject

    id::Nullable{Compat.String}
    date::DateTime
    target::BackupTarget
    extra::Dict
    driver::BackupDriver

    function BackupTargetRecoveryPoint(o::PyObject)
        new(o, o[:id], o[:date], BackupTarget(o[:target]), o[:extra], BackupDriver(o[:driver]))
    end
end
PyObject(o::BackupTargetRecoveryPoint) = o.o
convert(::Type{BackupTargetRecoveryPoint}, o::PyObject) =  BackupTargetRecoveryPoint(o)
show(io::IO, o::BackupTargetRecoveryPoint) = print(io, o.o[:__str__]())
@doc LazyHelp(_libcloud_backup_base, "BackupTargetRecoveryPoint") BackupTargetRecoveryPoint

const _backupdriver_fns = [
    :get_supported_target_types, :list_targets,
    :create_target, :create_target_from_node, :create_target_from_storage_container,
    :update_target, :delete_target,
    :list_recovery_points, :recover_target, :recover_target_out_of_place,
    :get_target_job, :list_target_jobs, :create_target_job, :resume_target_job, :suspend_target_job, :cancel_target_job
]

const _backup_target_fns = [ :update, :delete ]
const _backup_target_job_fns = [ :cancel, :suspend, :resume ]
const _backup_target_recovery_point_fns = [ :recover, :recover_to ]

for f in union(Set(_backupdriver_fns), Set(_base_driver_fns))
    sf = string(f)
    @eval @doc LazyHelp(_libcloud_backup_base, "BackupDriver", $sf) $(f)(b::BackupDriver, args...; kwargs...) = b.o[$(sf)](args..., kwargs...)
end

for f in _backup_target_fns
    sf = string(f)
    @eval @doc LazyHelp(_libcloud_backup_base, "BackupTarget", $sf) $(f)(b::BackupTarget, args...; kwargs...) = b.o[$(sf)](args..., kwargs...)
end

for f in _backup_target_job_fns
    sf = string(f)
    @eval @doc LazyHelp(_libcloud_backup_base, "BackupTargetJob", $sf) $(f)(b::BackupTargetJob, args...; kwargs...) = b.o[$(sf)](args..., kwargs...)
end

for f in _backup_target_recovery_point_fns
    sf = string(f)
    @eval @doc LazyHelp(_libcloud_backup_base, "BackupTargetRecoveryPoint", $sf) $(f)(b::BackupTargetRecoveryPoint, args...; kwargs...) = b.o[$(sf)](args..., kwargs...)
end

# Initialize cloud backup
function __init__()
    copy!(_libcloud_backup_types, pyimport("libcloud.backup.types"))
    copy!(_libcloud_backup_providers, pyimport("libcloud.backup.providers"))
    copy!(_libcloud_backup_base, pyimport("libcloud.backup.base"))

    global const BackupProvider = pywrap(_libcloud_backup_types[:Provider])
    global const BackupTargetType = pywrap(_libcloud_backup_types[:BackupTargetType])
    global const BackupTargetJobStatusType = pywrap(_libcloud_backup_types[:BackupTargetJobStatusType])

    _map_types(_libcloud_backup_base, (BackupTarget, BackupTargetJob, BackupTargetRecoveryPoint))
end

# types
export BackupProvider, BackupDriver, BackupTargetType, BackupTarget, BackupTargetJob, BackupTargetRecoveryPoint
# backup driver functions
export get_supported_target_types, list_targets,
       create_target, create_target_from_node, create_target_from_storage_container,
       update_target, delete_target,
       list_recovery_points, recover_target, recover_target_out_of_place,
       get_target_job, list_target_jobs, create_target_job, resume_target_job, suspend_target_job, cancel_target_job
# backup target functions
export update, delete
# backup target job functions
export cancel, suspend, resume
# backup target recovery point functions
export recover, recover_to

end # module Backup


text:
using SafeTestsets

@safetestset "Nystrom mesh" begin include("nystrommesh_test.jl") end

@safetestset "Kernels" begin include("kernels_test.jl") end

@safetestset "Density" begin include("density_test.jl") end

@safetestset "Potentials" begin include("potential_test.jl") end

@safetestset "Integral operators" begin include("integraloperator_test.jl") end

@safetestset "Density interpolation method" begin include("dim_test.jl") end

@safetestset "Discrete operators" begin include("discreteoperator_test.jl") end

# @safetestset "IFGF" begin include("ifgf_test.jl") end

max_stars_repo_name: WaveProp/Nystrom

max_issues_repo_name: WaveProp/Nystrom

max_forks_repo_name: WaveProp/Nystrom


text:
# # Interpolating `Point`s

# I made this doc with the command
# ```julia
# using Literate; Literate.markdown("examples/interpolating.jl",".";documenter=false)
# ```
# This document shows how to interpolate a `Vector{Point{T}}` instead of two standard `Vector{T}`:

module Points
using Interpolations
using BenchmarkTools

import Base: +, -, *, /

struct Point{T}
	x::T
	y::T
end

# Having defined the `Point` struct, we next define some arithmetic operations that will allow us to use the type with the `Interpolations.jl` package:

(+)(p1::Point, p2::Point) = Point(p1.x+p2.x, p1.y+p2.y)
(-)(p1::Point, p2::Point) = Point(p1.x-p2.x, p1.y-p2.y)
(*)(n::Number, p::Point) = Point(n*p.x, n*p.y)
(*)(p::Point, n::Number) = n*p
(/)(p::Point, n::Number) = Point(p.x/n, p.y/n)

Base.zero(::Type{Point{T}}) where {T} = Point(zero(T),zero(T))
Base.promote_rule(::Type{Point{T1}}, ::Type{T2}) where {T1,T2<:Number} = Point{promote_type(T1,T2)}
Base.promote_op(::typeof(*), ::Type{Point{T1}}, ::Type{T2}) where {T1,T2<:Number} = Point{promote_type(T1,T2)}
Base.promote_op(::typeof(*), ::Type{T1}, ::Type{Point{T2}}) where {T1<:Number,T2} = Point{promote_type(T1,T2)}

## (+)(p1::Point, p2::Point) = p1.y+p2.y
## (-)(p1::Point, p2::Point) = p1.y-p2.y
## (*)(n::Number, p::Point) = n*p.y
## (*)(p::Point, n::Number) = n*p
## (/)(p::Point, n::Number) =  p.y/n

## @noinline checksame(p1, p2) = p1.x == p2.x || throw(ArgumentError("arithmetic supported only for constant x, got $p1 and $p2"))
## function (+)(p1::Point, p2::Point)
##     # checksame(p1, p2)
##     Point(p1.x, p1.y + p2.y)
## end
## function (-)(p1::Point, p2::Point)
##     # checksame(p1, p2)
##     Point(p1.x, p1.y - p2.y)
## end
## (*)(n::Number, p::Point) = Point(n*p.x, n*p.y)
## (/)(n::Number, p::Point) = Point(n/p.x, n/p.y)
## (*)(p::Point, n::Number) = n*p

# Then we can move on to test the performance of this new way of interpolating. Notice that there are two concerns in the speed comparison:
# 1. how long does it take to interpolate a `Vector{Point{T}}` vs 2 `Vector{T}`s?
# 1. how long does it take to convert data in `Vector{Point{T}}` into 2 `Vector{T}`s?

n = 1_000_000
m = 10_000
# cast data as vector of Point
data = reinterpret(Point{Float64},hcat(collect(1:n),rand(n))',(n,))
# data to test interpolation on
rdata = rand(m)
# convert vec of points to separate x and y vectors
x = [i.x for i in data]
y = [i.y for i in data]

# my favourite version
# but also interpolates the x-grid :-(
function itp_points(d,rd)
	## itp = scale(interpolate(d, BSpline(Quadratic(Flat())),OnGrid()),linspace(d[1].x,d[end].x,length(d)))
	itp = scale(interpolate(d, BSpline(Linear()),OnGrid()),linspace(d[1].x,d[end].x,length(d)))
	out = similar(d)
	for i in 1:length(rd)
		out[i] = itp[rd[i]]
	end
	out
end

# traditional version
# only interpolate the y vector
function itp_x_y(x,y,rd)
	itp = scale(interpolate(y, BSpline(Linear()),OnGrid()),linspace(x[1],x[end],length(x)))
	out = similar(rd)
	for i in 1:length(rd)
		out[i] = itp[rd[i]]
	end
	out
end

# actual cost of traditional:
# need to first convert Point to vector!
# an in fact that conversion would have to happen FOR EACH interpolation
# in the below loop... so much more costly in the end!
function itp_pointsxy(d,rd)
	x = [i.x for i in data]
	y = [i.y for i in data]
	itp = scale(interpolate(y, BSpline(Linear()),OnGrid()),linspace(x[1],x[end],length(x)))
	out = similar(rd)
	for i in 1:length(rd)
		out[i] = itp[rd[i]]
	end
	out
end

# compile
itp_points(data[1:10],rdata[1:5])
itp_x_y(x[1:10],y[1:10],rdata[1:5])
itp_pointsxy(data[1:10],rdata[1:5])

# BenchmarkTools
t1 = @benchmark itp_points($data,$rdata)
t2 = @benchmark itp_x_y($x,$y,$rdata)
t3 = @benchmark itp_pointsxy($data,$rdata)

1+2

end

## We see here the timings for t1 and t2:

# ```julia
# julia v0.6.2> Points.t1
# BenchmarkTools.Trial:
#   memory estimate:  30.52 MiB
#   allocs estimate:  8
#   --------------
#   minimum time:     6.088 ms (9.38% GC)
#   median time:      8.402 ms (34.32% GC)
#   mean time:        8.675 ms (34.34% GC)
#   maximum time:     26.954 ms (28.22% GC)
#   --------------
#   samples:          576
#   evals/sample:     1
#
# julia v0.6.2> Points.t2
# BenchmarkTools.Trial:
#   memory estimate:  7.71 MiB
#   allocs estimate:  8
#   --------------
#   minimum time:     2.007 ms (0.00% GC)
#   median time:      2.213 ms (0.00% GC)
#   mean time:        3.514 ms (37.21% GC)
#   maximum time:     13.390 ms (77.77% GC)
#   --------------
#   samples:          1419
#   evals/sample:     1
#   ```

# # a hacked arithmetic for this type

# now lets redefine the module, in particular let's change the arithmetics:

module Points
using Interpolations
using BenchmarkTools

import Base: +, -, *, /

struct Point{T}
	x::T
	y::T
end

## (+)(p1::Point, p2::Point) = Point(p1.x+p2.x, p1.y+p2.y)
## (-)(p1::Point, p2::Point) = Point(p1.x-p2.x, p1.y-p2.y)
## (*)(n::Number, p::Point) = Point(n*p.x, n*p.y)
## (*)(p::Point, n::Number) = n*p
## (/)(p::Point, n::Number) = Point(p.x/n, p.y/n)

Base.zero(::Type{Point{T}}) where {T} = Point(zero(T),zero(T))
Base.promote_rule(::Type{Point{T1}}, ::Type{T2}) where {T1,T2<:Number} = Point{promote_type(T1,T2)}
Base.promote_op(::typeof(*), ::Type{Point{T1}}, ::Type{T2}) where {T1,T2<:Number} = Point{promote_type(T1,T2)}
Base.promote_op(::typeof(*), ::Type{T1}, ::Type{Point{T2}}) where {T1<:Number,T2} = Point{promote_type(T1,T2)}

(+)(p1::Point, p2::Point) = p1.y+p2.y
(-)(p1::Point, p2::Point) = p1.y-p2.y
(*)(n::Number, p::Point) = n*p.y
(*)(p::Point, n::Number) = n*p
(/)(p::Point, n::Number) =  p.y/n

n = 1_000_000
m = 10_000
data = reinterpret(Point{Float64},hcat(collect(1:n),rand(n))',(n,))
rdata = rand(m)
x = [i.x for i in data]
y = [i.y for i in data]

function itp_points(d,rd)
	## itp = scale(interpolate(d, BSpline(Quadratic(Flat())),OnGrid()),linspace(d[1].x,d[end].x,length(d)))
	itp = scale(interpolate(d, BSpline(Linear()),OnGrid()),linspace(d[1].x,d[end].x,length(d)))
	out = similar(rd)
	for i in 1:length(rd)
		out[i] = itp[rd[i]]
	end
	out
end

function itp_x_y(x,y,rd)
	itp = scale(interpolate(y, BSpline(Linear()),OnGrid()),linspace(x[1],x[end],length(x)))
	out = similar(rd)
	for i in 1:length(rd)
		out[i] = itp[rd[i]]
	end
	out
end

function itp_pointsxy(d,rd)
	x = [i.x for i in data]
	y = [i.y for i in data]
	itp = scale(interpolate(y, BSpline(Linear()),OnGrid()),linspace(x[1],x[end],length(x)))
	out = similar(rd)
	for i in 1:length(rd)
		out[i] = itp[rd[i]]
	end
	out
end

itp_points(data[1:10],rdata[1:5])
itp_x_y(x[1:10],y[1:10],rdata[1:5])
itp_pointsxy(data[1:10],rdata[1:5])

t1 = @benchmark itp_points($data,$rdata)
t2 = @benchmark itp_x_y($x,$y,$rdata)
t3 = @benchmark itp_pointsxy($data,$rdata)

end

# This shows that the *hacked* method is a bit faster, and has much less memory allocs.

# ```julia
# julia v0.6.2> Points.t1
# BenchmarkTools.Trial:
#   memory estimate:  15.34 MiB
#   allocs estimate:  8
#   --------------
#   minimum time:     3.760 ms (0.00% GC)
#   median time:      6.865 ms (13.17% GC)
#   mean time:        5.885 ms (26.24% GC)
#   maximum time:     20.627 ms (32.43% GC)
#   --------------
#   samples:          848
#   evals/sample:     1
#
# julia v0.6.2> Points.itp_points(Points.data[1:10],Points.rdata[1:5])
# 5-element Array{Float64,1}:
#  0.178303
#  0.171897
#  0.184147
#  0.179542
#  0.159741
#
# julia v0.6.2> Points.itp_x_y(Points.x[1:10],Points.y[1:10],Points.rdata[1:5])
# 5-element Array{Float64,1}:
#  0.178303
#  0.171897
#  0.184147
#  0.179542
#  0.159741
# ```