text:
"""
  plotlandmarksmatching(outdir)

Script for some plots after calling landmarksmatching.
"""
function plotlandmarksmatching(outdir)
  @warn "R packages tidyverse and ggforce should be installed."

  @rput outdir
  R"""
  setwd(outdir)
  widthfig <- 7.5
  library(tidyverse)
  library(ggforce)
  theme_set(theme_light())

  ########  read observations and process
  obs0df <- read_delim("obs0.csv", ";", escape_double = FALSE, trim_ws = TRUE)
  obsTdf <- read_delim("obsT.csv", ";", escape_double = FALSE, trim_ws = TRUE) %>% spread(key=pos,value=value) %>%
      mutate(shape=as.factor(shape))
  v0 <- bind_rows(obs0df, obs0df)
  vT <- bind_rows(obsTdf, obsTdf)

  #######  read noisefields
  nfsdf <- read_delim("noisefields.csv",";", escape_double = FALSE, trim_ws = TRUE)
  ######## read and analyse bridge iterates
  d <- read_table2("iterates.csv") %>%
    mutate(landmarkid=as.factor(landmarkid)) %>%
    gather(key="iteratenr",value="value",contains("iter"))
  d <- d %>% mutate(i = rep(1:(nrow(d)/4),each=4)) # add column to identify a particular combination (pos1,pos2,mom1,mom2) of one landmark at a specific time of a specific shape
  d <- d %>% spread(key=pqtype, value=value) %>% dplyr::select(-i) %>%  # can drop identifier column
    mutate(iteratenr=str_replace(iteratenr,"iter",""),iterate=as.numeric(str_replace(iteratenr,"iter","")))
  ####### read parameter updates
  parsdf <- read_delim("parameters.csv", ";", escape_double = FALSE, trim_ws = TRUE)
  ####### read acceptance information in mcmc-steps
  accdf <- read_delim("accdf.csv", ";", escape_double = FALSE, trim_ws = TRUE)

  # plots shapes and noisefields
  shapes <- ggplot() +
      geom_path(data=v0, aes(x=pos1,y=pos2), colour='black')+
      geom_path(data=vT, aes(x=pos1,y=pos2,group=shape), colour='orange') +
      geom_point(data=nfsdf, aes(x=locx, y=locy), color="Grey")+
      geom_circle(aes(x0 = locx, y0 = locy, r = nfstd), data = nfsdf,color="Grey",linetype="dotted")+
        theme(axis.title.x=element_blank(), axis.title.y=element_blank())
  pdf("shapes-noisefields.pdf",width=widthfig,height=4)
    show(shapes)
  dev.off()

  # plot overlaid landmark bridges
  bridges <- d %>%
   ggplot() +
      geom_path(aes(pos1,y=pos2,group=interaction(landmarkid,iteratenr),colour=iterate),size=0.5) +
      scale_colour_gradient(low="grey",high="darkblue")+
      geom_point(data=v0, aes(x=pos1,y=pos2), colour='black')+geom_point(data=vT, aes(x=pos1,y=pos2), colour='orange')+
      geom_path(data=v0, aes(x=pos1,y=pos2), colour='black',size=1.0)+geom_path(data=vT, aes(x=pos1,y=pos2,group=shape), colour='orange',size=1.0) +
      theme(axis.title.x=element_blank(), axis.title.y=element_blank()) +
      geom_point(data=nfsdf, aes(x=locx, y=locy), color="Grey")+
    geom_circle(aes(x0 = locx, y0 = locy, r = nfstd), data = nfsdf,color="Grey",linetype="dotted")

  pdf("bridges-overlaid.pdf",width=widthfig,height=4)
    show(bridges)
  dev.off()

  # traceplot of parameters
  tracepars <- parsdf %>% gather(key=par, value=value, a,  gamma)  %>%
    ggplot(aes(x=iterate, y=value)) + geom_path() + facet_wrap(~par, scales="free_y",ncol=1) +
   xlab("iterate") + ylab("")
  pdf("trace-pars.pdf",width=widthfig,height=4)
    show(tracepars)
  dev.off()

  # plot shape evolution
  ut <- unique(d$time)
  ind <- round(quantile(1:length(ut), (1:4)/4))
  times <- ut[ind]
  dhalfend <- bind_rows(d,d) %>%  dplyr::filter(time %in% times)

  shapes_evolution <-
    ggplot() +
    geom_path(data=dhalfend, aes(x=pos1,y=pos2,colour=iterate,group=iterate),alpha=0.5,size=0.5)+
      geom_point(data=v0, aes(x=pos1,y=pos2), colour='black')+
    geom_point(data=vT, aes(x=pos1,y=pos2), colour='orange')+
    geom_path(data=v0, aes(x=pos1,y=pos2), colour='black',size=0.6)+
    geom_path(data=vT, aes(x=pos1,y=pos2,group=shape), colour='orange',size=0.6)+
    facet_wrap(~time,ncol=2)+
    theme(axis.title.x=element_blank(), axis.title.y=element_blank())+
    scale_colour_gradient(low="grey",high="darkblue")+
  pdf("shapes_evolution.pdf", width = widthfig, height=4)
  show(shapes_evolution)
  dev.off()

  # plot acceptance rates
  accfig <- accdf %>% gather(key='kernel',value='acc',-iteration) %>%
        ggplot(aes(x=iteration,y=acc)) +
      geom_point(shape=124)+
      facet_wrap(~kernel)+xlab("iteration nr")+ylab("accept")
  pdf("acceptance.pdf",width=widthfig,height=4)
    show(accfig)
  dev.off()
  """
end

"""
  plottemplate_estimation(outdir)

Script for some plots after calling template_estimation
"""
function plottemplate_estimation(outdir)
  @warn "R packages tidyverse and ggforce should be installed."

  @rput outdir
  R"""
  setwd(outdir)
  widthfig <- 7.5
  library(tidyverse)
  library(ggforce)
  theme_set(theme_light())

  ########  read observations and process
  obsTdf <- read_delim("obsT.csv", ";", escape_double = FALSE, trim_ws = TRUE) %>% spread(key=pos,value=value) %>%
    mutate(shape=as.factor(shape))
  vT <- bind_rows(obsTdf, obsTdf)

  #######  read noisefields
  nfsdf <- read_delim("noisefields.csv",";", escape_double = FALSE, trim_ws = TRUE)

  ######## read and analyse bridge iterates
  d <- read_table2("iterates.csv") %>%
    mutate(landmarkid=as.factor(landmarkid)) %>%
    gather(key="iteratenr",value="value",contains("iter"))
  d <- d %>% mutate(i = rep(1:(nrow(d)/4),each=4)) # add column to identify a particular combination (pos1,pos2,mom1,mom2) of one landmark at a specific time of a specific shape
  d <- d %>% spread(key=pqtype, value=value) %>% dplyr::select(-i) %>%  # can drop identifier column
    mutate(iteratenr=str_replace(iteratenr,"iter",""),iterate=as.numeric(str_replace(iteratenr,"iter","")))
  d1 <- d %>% dplyr::filter(shapes==1)

  ####### read parameter updates
  parsdf <- read_delim("parameters.csv", ";", escape_double = FALSE, trim_ws = TRUE)

  ####### read acceptance information in mcmc-steps
  accdf <- read_delim("accdf.csv", ";", escape_double = FALSE, trim_ws = TRUE)

  # plots shapes and noisefields
  shapes <- ggplot() +
    geom_path(data=vT, aes(x=pos1,y=pos2,group=shape), colour='orange') +
    geom_point(data=nfsdf, aes(x=locx, y=locy), color="Grey")+
    geom_circle(aes(x0 = locx, y0 = locy, r = nfstd), data = nfsdf,color="Grey",linetype="dotted")+
    theme(axis.title.x=element_blank(), axis.title.y=element_blank())

  pdf("shapes-noisefields.pdf",width=widthfig,height=4)
  show(shapes)
  dev.off()

  # plot overlaid landmark bridges
  bridges <- d1 %>%
    ggplot() +
    geom_path(aes(pos1,y=pos2,group=interaction(landmarkid,iteratenr),colour=iterate),size=0.5) +
    scale_colour_gradient(low="grey",high="darkblue")+
    theme(axis.title.x=element_blank(), axis.title.y=element_blank()) +
    geom_point(data=nfsdf, aes(x=locx, y=locy), color="Grey")+
    geom_circle(aes(x0 = locx, y0 = locy, r = nfstd), data = nfsdf,color="Grey",linetype="dotted")+
    theme(legend.position='none')
  pdf("bridges-overlaid.pdf",width=widthfig,height=4)
  show(bridges)
  dev.off()

  # traceplot of parameters
  tracepars <- parsdf %>% gather(key=par, value=value, a, c, gamma)  %>%
    ggplot(aes(x=iterate, y=value)) + geom_path() + facet_wrap(~par, scales="free_y",ncol=1) +
    xlab("iterate") + ylab("")
  pdf("trace-pars.pdf",width=widthfig,height=4)
  show(tracepars)
  dev.off()

  # plot shape evolution
  ut <- unique(d$time)
  ind <- round(quantile(1:length(ut), (1:4)/4))
  times <- ut[ind]
  dhalfend <- bind_rows(d,d) %>%  dplyr::filter(time %in% times)

  shapes_evolution <-
    ggplot() +
    geom_path(data=dhalfend, aes(x=pos1,y=pos2,colour=iterate,group=iterate),alpha=0.5,size=0.5)+
    facet_wrap(~time,ncol=2)+
    geom_path(data=vT, aes(x=pos1,y=pos2,group=shape), colour='orange',size=0.4)+
    theme(axis.title.x=element_blank(), axis.title.y=element_blank())+
    scale_colour_gradient(low="grey",high="darkblue")

  pdf("shapes_evolution.pdf", width = widthfig, height=4)
    show(shapes_evolution)
  dev.off()

  # plot acceptance rates
  accfig <- accdf %>% gather(key='kernel',value='acc',-iteration) %>%
    ggplot(aes(x=iteration,y=acc)) +
    geom_point(shape=124)+
    facet_wrap(~kernel)+xlab("iteration nr")+ylab("accept")

  pdf("acceptance.pdf",width=widthfig,height=4)
    show(accfig)
  dev.off()

  # plot initial positions for all shapes (only interesting in case initial shape is unobserved)
  dtime0 <- d %>% dplyr::filter(time==0)
  # add factor for determining which phase of sampling
  dtime0 <-  dtime0 %>% mutate(phase =
      case_when(iterate < quantile(iterate,1/3) ~ "initial",
                between(iterate,quantile(iterate,1/3),quantile(iterate,2/3)) ~ "middle",
                iterate >= quantile(iterate,1/3) ~ "end")  ) %>% # reorder factor levels
    mutate(phase = fct_relevel(phase, "initial", "middle"))
  dtimeT <- d %>% dplyr::filter(time==1)

  initshapes0 <- ggplot()  +
    geom_point(data=vT, aes(x=pos1,y=pos2), colour='grey')+
    geom_path(data=vT, aes(x=pos1,y=pos2,group=shape), colour='grey', linetype="dashed",size=0.4) +
    geom_path(data=bind_rows(dtime0,dtime0),aes(x=pos1,y=pos2,colour=iterate)) +
    facet_wrap(~phase,ncol=1)

  pdf("initial_shapes.pdf")
    show(initshapes0)
  dev.off()

  observed_shapes = vT %>% dplyr::filter(landmark<=9) %>%
    ggplot(aes(x=pos1,y=pos2,label=landmark)) + geom_path() +
    geom_label() + facet_wrap(~shape) + xlab("") + ylab("")
  pdf("observed_shapes.pdf")
      show(observed_shapes)
  dev.off()
  """
  nothing
end


text:
# should be reimplement as macros by using Julia 1.5's @ccall macro
function igText(text)
    ccall((:igText, libcimgui), Cvoid, (Cstring, Cstring), "%s", text)
end

function igTextColored(col, text)
    ccall((:igTextColored, libcimgui), Cvoid, (ImVec4, Cstring), col, text)
end

function igTextDisabled(text)
    ccall((:igTextDisabled, libcimgui), Cvoid, (Cstring,), text)
end

function igTextWrapped(text)
    ccall((:igTextWrapped, libcimgui), Cvoid, (Cstring,), text)
end

function igLabelText(label, text)
    ccall((:igLabelText, libcimgui), Cvoid, (Cstring, Cstring), label, text)
end

function igBulletText(text)
    ccall((:igBulletText, libcimgui), Cvoid, (Cstring,), text)
end

function igTreeNodeStrStr(str_id, text)
    ccall((:igTreeNodeStrStr, libcimgui), Bool, (Cstring, Cstring), str_id, text)
end

function igTreeNodePtr(ptr_id, text)
    ccall((:igTreeNodePtr, libcimgui), Bool, (Ptr{Cvoid}, Cstring), ptr_id, text)
end

function igTreeNodeExStrStr(str_id, flags, text)
    ccall((:igTreeNodeExStrStr, libcimgui), Bool, (Cstring, ImGuiTreeNodeFlags, Cstring), str_id, flags, text)
end

function igTreeNodeExPtr(ptr_id, flags, text)
    ccall((:igTreeNodeExPtr, libcimgui), Bool, (Ptr{Cvoid}, ImGuiTreeNodeFlags, Cstring), ptr_id, flags, text)
end

function igSetTooltip(text)
    ccall((:igSetTooltip, libcimgui), Cvoid, (Cstring,), text)
end

function igOpenPopup(str_id)
    ccall((:igOpenPopup, libcimgui), Cvoid, (Cstring,), str_id)
end

function igLogText(text)
    ccall((:igLogText, libcimgui), Cvoid, (Cstring,), text)
end

ext: jl


text:
cases = [
# Calls
:(f(x)) "(call f x)"
:(f(x, y=1, z=2)) "(call f x (kw y 1) (kw z 2))"
:(f(x; y=1)) "(call f (parameters (kw y 1)) x)"
:(f(x...)) "(call f (... x))"
:(f(x) do a,b
    body
end) "(do (call f x) (-> (tuple a b) (block body)))"
# Operators
:(x+y) "(call + x y)"
:(a+b+c+d) "(call + a b c d)"
:(2x) "(call * 2 x)"
:(a&&b) "(&& a b)"
:(x += 1) "(+= x 1)"
:(a ? 1 : 2) "(if a 1 2)"
:(a:b) "(call : a b)" #
:(a:b:c) "(call : a b c)" #
:(a,b) "(tuple a b)"
:(a==b) "(call == a b)"
:(1<i<=n) "(comparison 1 < i <= n)"
:(a.b) "(. a (quote b))"
:(a.(b)) "(. a (tuple b))"
# Bracketed forms
:(a[i]) "(ref a i)"
:(t[i;j]) "(typed_vcat t i j)"
:(t[i j]) "(typed_hcat t i j)"
:(t[a b; c d]) "(typed_vcat t (row a b) (row c d))"
:(a{b}) "(curly a b)"
:(a{b;c}) "(curly a (parameters c) b)"
:([x]) "(vect x)"
:([x,y]) "(vect x y)"
:([x;y]) "(vcat x y)"
:([x y]) "(hcat x y)"
:([x y; z t]) "(vcat (row x y) (row z t))"
:([x for y in z, a in b]) "(comprehension (generator x (= y z) (= a b)))" #
:(T[x for y in z]) "(typed_comprehension T (generator x (= y z)))" #
:((a, b, c)) "(tuple a b c)"
:((a; b; c)) "(block a b c)" #
# Macros
:(@m x y) "(macrocall @m x y)"
:(Base.@m x y) "(macrocall (. Base (quote @m)) x y)"
:(@Base.m x y) "(macrocall (. Base (quote @m)) x y)"
# Strings
:("a") "\"a\""
:(x"y") "(macrocall @x_str \"y\")"
:(x"y"z) "(macrocall @x_str \"y\" \"z\")"
:("x = $x") "(string \"x = \" x)"
:(`a b c`) "(macrocall (|.| Core @cmd) \"a b c\")"
(quote
"some docs"
f(x) = x
end).args[2] "(macrocall (|.| Core @doc) \"some docs\" (= (call f x) (block x)))"
# Imports and such
:(import a) "(import (. a))"
:(import a.b.c) "(import (. a b c))"
:(import ...a) "(import (. . . . a))"
:(import a.b, c.d) "(import (. a b) (. c d))"
:(import Base: x) "(import (: (. Base) (. x)))"
:(import Base: x, y) "(import (: (. Base) (. x) (. y)))"
:(export a, b) "(export a b)"
# Numbers
:(11111111111111111111) "(macrocall (|.| Core @int128_str) nothing \"11111111111111111111\")"
:(0xfffffffffffffffff) "(macrocall (|.| Core @uint128_str) nothing \"0xfffffffffffffffff\")"
:(1111111111111111111111111111111111111111) "(macrocall (|.| Core @big_str) nothing \"1111111111111111111111111111111111111111\")"
# Blocks
(quote
if a
    b
elseif c
    d
else
    e
end
end).args[2] "(if a (block b) (elseif (block c) (block d) (block e)))"
:(while condition; body; end) "(while condition (block body))" #
:(for var = iter; body; end) "(for (= var iter) (block body))" #
:(for v1=iter1, v2=iter2; body; end) "(for (block (= v1 iter1) (= v2 iter2)) (block body))" #
:(break) "(break)"
:(continue) "(continue)"
:(let var=iter; body; end) "(let (= var iter) (block body))"
:(let v1=iter1, v2=iter2; body; end) "(let (block (= v1 iter1) (= v2 iter2)) (block body))"
:(f(x) = body) "(= (call f x) (block body))" #
:(function f(x::T; k = 1) where T
    return x+1
end) "(function (where (call f (parameters (kw k 1)) (:: x T)) T) (block (return (call + x 1))))"
:(mutable struct Foo{T<:S}
    x::T
end) "(struct true (curly Foo (<: T S)) (block (:: x T)))"
:(try; try_block;catch var; catch_block; finally finally_block; end) "(try (block try_block) var (block catch_block) (block finally_block))"
:(try; try_block;catch; catch_block; finally finally_block; end) "(try (block try_block) false (block catch_block) (block finally_block))"
:(try; try_block;catch var; catch_block; end) "(try (block try_block) var (block catch_block))"
:(try; try_block;catch; catch_block; end) "(try (block try_block) false (block catch_block))" #
]

@testset "AST example $i" for i = 1:size(cases, 1)
    @test pexpressions(cases[i,1]) == cases[i,2]
end


text:
module AdaptTests
    using FactCheck
    using Adapt

    # run tests
    @fact l1_norm(1:10) --> 55
    @fact l2_norm(1:10) --> 19.621416870348583
    @fact l2_norm(1:10) --> 19.621416870348583
    @fact lp_norm(1:10,2) --> 19.621416870348583
    @fact lp_norm(1:10,3) --> 14.462447418811122
    @fact lp_norm(1:10,0.5) --> 504.82352465265495
    # check negative values are properly handled
    b = collect(1:10)
    b[2] = -2
    @fact lp_norm(b,1) --> 55

    FactCheck.exitstatus()
end


text:
struct BordingPass
    row::Int
    seat::Int
    id::Int
end

function createBinaryString(string, zero::Char, one::Char)
    bstring = ""
    for c in string
        if c == zero
            bstring = bstring*"0"
        else
            bstring = bstring*"1"
        end
    end
    return bstring
end

function rowNumber(rowString)
    binaryRow = createBinaryString(rowString, 'F', 'B')
    return parse(Int, binaryRow, base=2)
end

function seatNumber(seatString)
    binarySeat = createBinaryString(seatString, 'L', 'R')
    return parse(Int, binarySeat, base=2)
end

function createBordingPass(line)
    rowString::String = SubString(line, 1:7)
    seatString::String = SubString(line, 8:10)
    row = rowNumber(rowString)
    seat = seatNumber(seatString)
    id = row * 8 + seat
    return BordingPass(row,seat,id)
end

open("data.txt") do file
    maxId = 0
    for line in eachline(file)
        bordingPass = createBordingPass(line)
        if bordingPass.id > maxId
            maxId = bordingPass.id
        end
        println(bordingPass)
    end
    println("maxId=",maxId)
end


text:
"""
    make_widgets(x, label)

Widets are built recursively from any type, labelled using the dict or struct
label, AxisArray names, etc.

They can then be edited before retrieving signals and building an interface.
"""
make_widgets(x, label...) = begin
    widgets = []
    labels = []
    if length(fieldnames(x)) > 0
        for name in fieldnames(x)
            field = getfield(x, name)
            widget = make_widgets(field, string(name))
            if widget != nothing
                push!(widgets, widget)
                push!(labels, name)
            end
        end
        if length(widgets) > 0
            if all(w -> typeof(w) == typeof(widgets[1]), widgets)
                widgettype = typeof(widgets[1])
            else
                widgettype = Any
            end
            return OrderedDict{Any,widgettype}(zip(labels, widgets))
        else
            return nothing
        end
    else
        return nothing
    end
end
make_widgets(x::Associative, label...) = begin
    widgets = []
    labels = []
    for key in keys(x)
        field = x[key]
        widget = make_widgets(field, string(key))
        if widget != nothing
            push!(widgets, widget)
            push!(labels, key)
        end
    end
    if length(widgets) > 0
        if all(w -> typeof(w) == typeof(widgets[1]), widgets)
            widgettype = typeof(widgets[1])
        else
            widgettype = Any
        end
        return OrderedDict{Any,widgettype}(zip(labels, widgets))
    else
        return nothing
    end
end
make_widgets(xs::Tuple) = begin
    label = defaults[:missing_label]
    make_widgets(xs, label)
end
make_widgets(xs::Tuple, label) = begin
    labels = [string(label, n) for n = 1:length(xs)]
    tuple(make_widgets.(xs, labels)...)
end
make_widgets(xs::AbstractArray, label) =
    if length(xs) <= defaults[:maxlen_widget_array]
        labels = [string(label, n) for n = 1:length(xs)]
        make_widgets.(xs, labels)
    end
make_widgets(x::Number, label) = begin
    range = get_range(x, label)
    return make_widgets(x, label, range)
end
@require Interact begin
    make_widgets(x::Number, label, range) = begin
        widget = slider(range, value = x, label = label)
        widget.signal.name = label * "_slider"
        return widget
    end
    make_widgets(x::Number, label, range) = begin
        widget = slider(range, value = x, label = label)
        widget.signal.name = label * "_slider"
        return widget
    end
    make_widgets(x::Bool, label) = begin
        widget = togglebutton(label, value = x)
        widget.signal.name = label * "_toggle"
        return widget
    end
end
@require InteractNext begin
    make_widgets(x::Number, label, range) = begin
        widget = slider(range, value = x, label = label)
        return widget
    end
    make_widgets(x::Bool, label) = begin
        widget = checkbox(x, label=label)
        return widget
    end
end
make_widgets(x::Function, label) = nothing
@require DataFrames begin
make_widgets(x::DataFrame, label...) = begin
    array = []
    for name in names(x)
        field = x[name]
        widget = make_widgets(field, string(name))
        if widget != nothing
            push!(array, widget)
        end
    end
    # Return 2 dimensional array for column output
    return reshape(array, 1, length(array))
end
end
make_widgets(xs::AxisArray, label...) = nothing
@require StaticArrays begin
    make_widgets(xs::FieldVector, label...) = begin
        names = string.(fieldnames(xs))
        [make_widgets(x, n) for (x, n) in zip(xs, names)]
    end
end

"""
    make_plottables(x, label)
Generate 'plottable' toggles recursively, with optional label.
"""
function make_plottables end

make_plottables(x::Associative, label...) = begin
    dict = OrderedDict()
    for key in keys(x)
        field = x[key]
        widget = make_plottables(field, string(key))
        if widget != nothing
            dict[key] = widget
        end
    end
    return length(keys(dict)) > 0 ? dict : nothing
end
make_plottables(x::Any, label...) = begin
    dict = OrderedDict()
    if length(fieldnames(x)) > 0
        for name in fieldnames(x)
            field = getfield(x, name)
            widget = make_plottables(field, string(name))
            if widget != nothing
                dict[name] = widget
            end
        end
    end
    return length(keys(dict)) > 0 ? dict : nothing
end
make_plottables(xs::Tuple) = make_plottables(xs, "unknown")
make_plottables(xs::Tuple, label) = begin
    labels = [string(label, n) for n = 1:length(xs)]
    widgets = tuple(make_plottables.(xs, labels)...)
    if all(w -> w == nothing, widgets)
        return nothing
    else
        return widgets
    end
end
@require DataFrames begin
make_plottables(x::DataFrame, label...) = begin
    labels = []
    values = []
    for name in names(x)
        field = x[name]
        widget = make_widgets(false, string(name))
        if widget != nothing
            push!(labels, name)
            push!(values, widget)
        end
    end
    # Return 2 dimensional array for column output
    return AxisArray(reshape(values, 1, length(values)),
                     Axis{:row}([:row]), Axis{:col}(labels))
end
end
make_plottables(xs::AxisArray, label...) = make_toggles(xs)
# Arrays get a toggle to control subplots
make_plottables(x::Vector, label) =
    if length(x) >= defaults[:minlen_plottable_array]
        make_widgets(false, label)
    end

make_toggles(x::A) where A<:AxisArray = begin
    a = [[string(a[i]) for i in 1:length(a)] for a in axes(x)]
    num_axes = length(a)
    if num_axes == 2
        return make_toggles(a[1])
    elseif num_axes == 3
        return make_toggles(a[1], a[2])
    else
        return nothing
    end
end
@require StaticArrays begin
    make_toggles(x::FieldVector) = make_toggles(fieldnames(x))
end
@require Interact begin
    make_toggles(colnames::AbstractArray) =
        [togglebutton(string(colnames[c])) for c = 1:length(colnames)]
    make_toggles(rownames::AbstractArray{String}, colnames::AbstractArray{String}) = begin
        rows = length(rownames)
        cols = length(colnames)
        return [togglebutton(join([string(rownames[r]), string(colnames[c])], ",")) for r = 1:rows, c = 1:cols]
    end
end
@require InteractNext begin
    make_toggles(colnames::AbstractArray) =
        [checkbox(false, label=string(colnames[c])) for c = 1:length(colnames)]
    make_toggles(rownames::AbstractArray{String}, colnames::AbstractArray{String}) = begin
        rows = length(rownames)
        cols = length(colnames)
        return [checkbox(false, label = join([string(rownames[r]), string(colnames[c])], ",")) for r = 1:rows, c = 1:cols]
    end
end

"""
    make_interface(xs; box)

Interfaces are build recursively from any type containing widgets. This
mostly means anything returned by make_widgets, but also after having any
sub-components deleted custom components added:
"""
make_interface(xs::Tuple; box = get_box(xs)) = begin
    widgets = []
    for x in xs
        widget = make_interface(x)
        if widget != nothing
            widget = push!(widgets, widget)
        end
    end
    box(widgets...)
end

make_interface(xs::AbstractArray{T,1}; box = get_box(xs)) where T =
    spreadwidgets(make_interface.(xs); cols = 2)
make_interface(xs::AbstractArray{T,2}; box = get_box(xs)) where T = arrange_columns(xs; box = box)
@require Interact begin
    make_interface(xs::AbstractVector{Interact.Slider{T}}; box = get_box(xs)) where T =
        spreadwidgets(make_interface.(xs); cols = 2)
    make_interface(xs::AbstractVector{Interact.ToggleButton}; box = get_box(xs)) =
        box(make_interface.(xs)...)
    make_interface(x::Interact.Widget; box = get_box(x)) = x
    make_interface(x::OrderedDict{Any,Interact.Slider{T}}; box = get_box(x)) where T = begin
        widgets = []
        for (key, val) in x
            if val != nothing
                push!(widgets, make_interface(val))
            end
        end
        return spreadwidgets(widgets; cols = 3)
    end
end
@require InteractNext begin
    get_box(x::WebIO.Scope) = hbox
    make_interface(xs::AbstractVector{WebIO.Scope}; box=get_box(xs)) =
        spreadwidgets(make_interface.(xs); cols = 2)
    make_interface(x::WebIO.Scope; box=get_box(x)) = x
end
make_interface(x::OrderedDict; box=get_box(x)) = begin
    widgets = []
    for (key, val) in x
        if val != nothing
            push!(widgets, make_interface(val))
        end
    end
    return box(widgets...)
end
make_interface(xs::Void; box = vbox) = nothing

function arrange_columns(xs; box = get_box(xs))
    columns = []
    for col in 1:size(xs, 2)
        push!(columns, vbox(make_interface.(xs[:,col])...))
    end
    return box(columns...)
end

"""
    spreadwidgets(widgets; cols = 6)

Create multiple columns to spread widgets accross a page
"""
function spreadwidgets(widgets; cols = 6)
    vboxes = []
    widget_col = []
    colsize = ceil(Int,length(widgets)/cols)
    # Build vbox columns for widgets
    for i = 1:length(widgets)
        push!(widget_col, widgets[i])
        if rem(i, colsize) == 0
            push!(vboxes, vbox(widget_col...))
            widget_col = []
        end
    end

    # Push remaining widgets
    push!(vboxes, vbox(widget_col...))
    return hbox(vboxes...)
end


text:
# Script to generate part of the temporal splitting experiments from Chen and Mišić (see Section EC.5 of the paper) with IRI data.
#
# This script is configured to be run on Amazon EC2 instances. It is accompanied by two other scripts: a helperfn (helper function) script, R2_IRI_temporal_helperfn.jl,
# and an include script, R2_IRI_temporal_include.jl.
#
# The helper function script contains all the code necessary to load data for a single instance, run all the necessary methods on the instance,
# and then to record the results.
#
# The include script contains all the include() statements to make sure the helper function script is equipped with all the necessary functions.
#
# The present script, R2_IRI_temporal.jl, is effectively a custom distributed job scheduler. It creates a list of "jobs" (different instances; for this temporal splitting experiment,
# this is just one of the product categories).
# It then loops through these jobs, assigning batches of these jobs to newly created processors. Once the processors
# finish, the processors are killed and new ones are created for the next batch of jobs. This loop continues until all jobs have been completed. This architecture was necessary due
# to issues with running Julia v0.6's pmap() on Amazon EC2.
#
# To run it, navigate to the mdpart_exec_scripts/ directory, and type
#
# julia R2_IRI_temporal.jl
#
# in the command line. On Amazon EC2, a more useful way to run it is
#
# nohup julia R2_IRI_temporal.jl &
#
# The results will be placed in a folder under tcm_testing/R2_IRI_temporal_test/ . The out_*.csv result files contain the performance of the different methods.
#
# This script assumes that you have the IRI data files under tcm_data/IRI_data_R2_temporally_split/ and tcm_data/IRI_data_assortments/ . These files are not provided in this code
# due to a confidentiality agreement with IRI.
#
# In addition to this script, please refer to R2_IRI_temporalRCS.jl (which executes the ranking-based model with tuned consideration set size) and R2_IRI_temporalRTS.jl
# (which executes the randomized tree sampling method for the decision forest model).
#
# At the end of the script, run(`sudo shutdown -h now`) will initiate the shutdown behavior on the machine. For some (most?) Amazon EC2 instances, this should cause the
# instance to "stop", to prevent you from accruing more EC2 fees; you can always start the instance later to download the results to your computer.
# *** Double check to make sure this is the default behavior for you (you default may be that shutdown causes the EC2 instance to terminate!!!). ***

include("R2_IRI_temporal_include.jl")

machine = 2;
expdirpath = "../tcm_testing/R2_IRI_temporal_test/"
if (!isdir(expdirpath))
	mkdir(expdirpath);
end

target_exec_script = string(expdirpath, "exec_script.jl");
if (isfile(target_exec_script))
	rm(target_exec_script)
end
cp(@__FILE__, target_exec_script)

initial_helper_script = string( pwd(), "/R2_IRI_temporal_helperfn.jl")
target_helper_script = string(expdirpath, "helperfn.jl")
if (isfile(target_helper_script))
	rm(target_helper_script)
end
cp(initial_helper_script, target_helper_script)

# forest_cv_outcsvfilepath = string(expdirpath, "out_forest_cv.csv");
# forest_cv_outcsvhandle = open(forest_cv_outcsvfilepath, "w")

# row_string = string("data_name",  ",",
# 							"test_fold", ",",
# 							"num_train_assortments", ",",
# 							"num_test_assortments", ",",
# 							"num_train_transactions", ",",
# 							"num_test_transactions", ",",
# 							"proc_id", ",",
# 							"attempt", ",",
# 							"method", ",",
# 							"best_d", ",",
# 							"best_KL_avg_cv", ",",
# 							"elapsed_time")
# row_string = string(row_string, "\n");
# print(forest_cv_outcsvhandle, row_string);
# flush(forest_cv_outcsvhandle);

forest_outcsvfilepath = string(expdirpath, "out_forest.csv");
forest_outcsvhandle = open(forest_outcsvfilepath, "w")

row_string = string("data_name",  ",",
							"num_train_assortments", ",",
							"num_test_assortments", ",",
							"num_train_transactions", ",",
							"num_test_transactions", ",",
							"proc_id", ",",
							"attempt", ",",
							"method", ",",
							"depth", ",",
							"loglik", ",",
							"K", ",",
							"perTol", ",",
							"test_KL", ",",
							"noinf_count", ",",
							"test_KL_avg", ",",
							"elapsed_time", ",",
							"simple_avg_leaves", ",",
							"weighted_avg_leaves", ",",
							"max_leaves", ",",
							"simple_avg_depth", ",",
							"weighted_avg_depth", ",",
							"max_depth")
row_string = string(row_string, "\n");
print(forest_outcsvhandle, row_string);
flush(forest_outcsvhandle);

ranking_outcsvfilepath = string(expdirpath, "out_ranking.csv");
ranking_outcsvhandle = open(ranking_outcsvfilepath, "w")
row_string = string("data_name",  ",",
							"num_train_assortments", ",",
							"num_test_assortments", ",",
							"num_train_transactions", ",",
							"num_test_transactions", ",",
							"proc_id", ",",
							"attempt", ",",
							"method", ",",
							"loglik", ",",
							"K", ",",
							"perTol", ",",
							"test_KL", ",",
							"noinf_count", ",",
							"test_KL_avg", ",",
							"elapsed_time", ",",
							"avg_CS_size", ",",
							"weighted_avg_CS_size", ",",
							"max_CS_size")
row_string = string(row_string, "\n");
print(ranking_outcsvhandle, row_string);
flush(ranking_outcsvhandle);

MNL_outcsvfilepath = string(expdirpath, "out_MNL.csv");
MNL_outcsvhandle = open(MNL_outcsvfilepath, "w")
row_string = string("data_name",  ",",
							"num_train_assortments", ",",
							"num_test_assortments", ",",
							"num_train_transactions", ",",
							"num_test_transactions", ",",
							"proc_id", ",",
							"attempt", ",",
							"method", ",",
							"loglik", ",",
							"K", ",",
							"perTol", ",",
							"numInitializations", ",",
							"test_KL", ",",
							"noinf_count", ",",
							"test_KL_avg", ",",
							"elapsed_time")
row_string = string(row_string, "\n");
print(MNL_outcsvhandle, row_string);
flush(MNL_outcsvhandle);

# LCMNL_cv_outcsvfilepath = string(expdirpath, "out_LCMNL_cv.csv");
# LCMNL_cv_outcsvhandle = open(LCMNL_cv_outcsvfilepath, "w")

# row_string = string("data_name",  ",",
# 							"test_fold", ",",
# 							"num_train_assortments", ",",
# 							"num_test_assortments", ",",
# 							"num_train_transactions", ",",
# 							"num_test_transactions", ",",
# 							"proc_id", ",",
# 							"attempt", ",",
# 							"method", ",",
# 							"best_K", ",",
# 							"perTol", ",",
# 							"numInitializations", ",",
# 							"best_KL_avg_cv", ",",
# 							"elapsed_time")
# row_string = string(row_string, "\n");
# print(LCMNL_cv_outcsvhandle, row_string);
# flush(LCMNL_cv_outcsvhandle);

#nfolds = 5
datanames = [ "diapers","paptowl","peanbutr","coldcer","toitisu","fzdinent","mayo","hhclean",
			"carbbev","coffee","fzpizza","saltsnck","mustketc","margbutr", "deod","yogurt","factiss","toothpa","shamp","soup",
			"milk","spagsauc","sugarsub","photo",
			"toothbr","laundet","beer",
			"cigets","hotdog","blades"]

# datanames = ["photo", "blades", "diapers"]

instances_vec = Any[];

# nfolds = 5

for single_data_name in datanames
	# assortments = convert(Array{Int64}, readcsv( string("../tcm_data/IRI_data_assortments/", single_data_name, "_assortments.csv") ) )
	#M = size(assortments,1)
	# for fold in 1:nfolds
		push!(instances_vec, single_data_name)
	# end
end

numInstances = length(instances_vec);

par_size = 10;

todo = collect(1:numInstances);

current_jobs = zeros(Int64, 0);
current_jobs_futures = Future[];
current_procs = zeros(Int64, 0)
isfinished = zeros(Bool, numInstances);

results = zeros(Float64, numInstances);
failed_attempts = zeros(Int64,numInstances);

overall_start_time = time()

iter = 0;

while ( !isempty(todo) )
	iter += 1
	numRemaining = length(todo);
	println("Starting iteration $iter ... ($numRemaining left) ")

	for i in 1:par_size
		if (isempty(todo))
			break;
		end
		single_job = pop!(todo)
		push!(current_jobs, single_job)
		newproc = addprocs(1)[1];
		fetch(@spawnat newproc include("R2_IRI_temporal_include.jl"))

		temp = @spawnat newproc helper_fn(instances_vec[single_job], failed_attempts[single_job])
		push!(current_jobs_futures, temp )
		push!(current_procs, newproc)
	end

	for i in 1:par_size
		if (isempty(current_jobs_futures))
			break;
		end
		single_future = pop!(current_jobs_futures);
		single_job = pop!(current_jobs)
		single_proc = pop!(current_procs)

		# single_result = fetch(single_future)
		# @show single_result
		# results[single_job] = single_result
		# rmprocs(single_proc)

		try
			single_result = fetch(single_future)
			@show single_result
			results[single_job] = single_result
			rmprocs(single_proc)
		catch y
			println("\t Error!")
			@show y
			push!(todo, single_job)
			failed_attempts[single_job] += 1;
			rmprocs(single_proc)
		end
	end

	# if (iter == 5)
	# 	break;
	# end

end

println("Finished")
@show results
@show failed_attempts

@show nprocs()

@show sum(results)
total_helper_time = sum(results);

overall_end_time = time();

overall_elapsed_time = overall_end_time - overall_start_time

@show overall_elapsed_time

timingcsvfilepath = string(expdirpath, "out_timing.txt");
timingcsvhandle = open(timingcsvfilepath, "w")
print(timingcsvhandle, string("sumhelper_time_vec,",  total_helper_time, "\n") )
print(timingcsvhandle, string("elapsed_time,", overall_elapsed_time, "\n") )
close(timingcsvhandle)

close(MNL_outcsvhandle);
# close(LCMNL_cv_outcsvhandle);
close(forest_outcsvhandle);
# close(forest_cv_outcsvhandle);
close(ranking_outcsvhandle);

run(`sudo shutdown -h now`)

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_licenses: MIT Unlicense


text:
include("DSGEmoments.jl")
include("DSGEmodel.jl") # defines prior and log-likelihood
using LinearAlgebra, Statistics, DelimitedFiles
function main()
    data = readdlm("dsgedata.txt")
    tuning = [0.001, 0.15, 0.1, 0.006, 0.15, 0.004, 0.005] # fix this somehow
    include("parameters.jl")
    lb = lb_param_ub[:,1]
    ub = lb_param_ub[:,3]
    θtrue = lb_param_ub[:,2]
    θinit = (ub + lb) / 2.0 # prior mean to start
    lnL = θ -> logL(θ, data)
    Prior = θ -> prior(θ) # uniform, doesn't matter
    # define things for MCMC
    verbosity = true
    burnin = 10000
    ChainLength = 50000
    # initial proposal moves one at a time
    Proposal = θ -> proposal1(θ, tuning)
    tuning = [0.0001, 0.1, 0.01, 0.001, 0.05, 0.001, 0.001]
    chain = mcmc(θinit, ChainLength, burnin, Prior, lnL, Proposal, verbosity)
    # keep every 10th
    i = 1:size(chain,1)
    keep = mod.(i,10.0).==0
    chain = chain[keep,:]
    # now use a MVN random walk proposal
    Σ = cov(chain[:,1:7])
    tuning = 0.1
    for j = 1:10
        P = (cholesky(Σ)).U
        Proposal = θ -> proposal2(θ,tuning*P)
        θinit = vec(mean(chain[:,1:7],dims=1))
        if j == 10
            ChainLength = 200000
        end
        chain = mcmc(θinit, ChainLength, 0, Prior, lnL, Proposal, verbosity)
        accept = mean(chain[:,end])
        #println("tuning: ", tuning, "  acceptance rate: ", accept)
        if accept > 0.35
            tuning *= 1.5
        elseif accept < 0.25
            tuning *= 0.5
        end
        # keep every 10th
        i = 1:size(chain,1)
        keep = mod.(i,10.0).==0
        θinit = vec(mean(chain[:,1:7],dims=1))
        Σ = 0.5*Σ + 0.5*cov(chain[:,1:7])
    end
    # keep every 10th to reduce autocorrelation
    i = 1:size(chain,1)
    j = mod.(i,10.0).==0
    chain = chain[j,:]
    # plain MCMC fit
    posmean = vec(mean(chain[:,1:7],dims=1))
    inci = zeros(7)
    for i = 1:7
        lower = quantile(chain[:,i],0.05)
        upper = quantile(chain[:,i],0.95)
        inci[i] = θtrue[i] >= lower && θtrue[i] <= upper
    end
    p = npdensity(chain[:,1]) # example of posterior plot
    display(p)
    prettyprint([posmean inci])
    return chain
end
main();

max_stars_repo_name: nilshg/Econometrics

max_issues_repo_name: nilshg/Econometrics

max_forks_repo_name: nilshg/Econometrics


text:
export write_sxml

function full_resp(xe::XMLElement)
  resp = MultiStageResp(24)
  ns = 0
  nmax = 0
  gain = one(Float64)
  f0 = one(Float64)
  xr = child_elements(xe)
  units = ""

  for r in xr
    if name(r) == "InstrumentSensitivity"
      for ii in child_elements(r)
        nii = name(ii)
        if nii == "Value"
          gain = parse(Float64, content(ii))
        elseif nii == "Frequency"
          f0 = parse(Float64, content(ii))
        elseif nii == "InputUnits"
          units = units2ucum(fix_units(content(get_elements_by_tagname(ii, "Name")[1])))
        end
      end

    # stages
    elseif name(r) == "Stage"
      ns = parse(Int64, attribute(r, "number"))
      nmax = max(ns, nmax)
      resp_code = 0x00
      c = one(Float64)
      a0 = one(Float64)
      f0 = one(Float64)
      p = Array{Complex{Float64},1}(undef, 0)
      z = Array{Complex{Float64},1}(undef, 0)
      num = Array{Float64,1}(undef, 0)
      den = Array{Float64,1}(undef, 0)
      units_in = ""
      units_out = ""

      for ii in child_elements(r)
        nii = name(ii)
        if nii == "PolesZeros"
          resp_code = 0x02
          for iii in child_elements(ii)
            niii = name(iii)
            if niii == "PzTransferFunctionType"
              c = content(iii) == "LAPLACE (RADIANS/SECOND)" ? c : Float64(2pi)

            # Zero, Pole
            elseif niii in ["Zero", "Pole"]
              pzv = complex(
                parse(Float64, content(get_elements_by_tagname(iii, "Real")[1])),
                parse(Float64, content(get_elements_by_tagname(iii, "Imaginary")[1]))
              )
              if niii == "Zero"
                push!(z, pzv)
              else
                push!(p, pzv)
              end

            elseif niii == "InputUnits"
              resp.i[ns] = fix_units(content(get_elements_by_tagname(iii, "Name")[1]))

            elseif niii == "OutputUnits"
              resp.o[ns] = fix_units(content(get_elements_by_tagname(iii, "Name")[1]))

            elseif niii == "NormalizationFactor"
              a0 = parse(Float64, content(iii))

            elseif niii == "NormalizationFrequency"
              f0 = parse(Float64, content(iii))
            end

          end

        elseif nii == "StageGain"
          resp.gain[ns] = parse(Float64, content(get_elements_by_tagname(ii, "Value")[1]))
          resp.fg[ns] = parse(Float64, content(get_elements_by_tagname(ii, "Frequency")[1]))

        elseif nii == "Coefficients" || nii == "FIR"
          resp_code = 0x03
          for iii in child_elements(ii)
            niii = name(iii)
            if niii == "InputUnits"
              resp.i[ns] = fix_units(content(get_elements_by_tagname(iii, "Name")[1]))

            elseif niii == "OutputUnits"
              resp.o[ns] = fix_units(content(get_elements_by_tagname(iii, "Name")[1]))

            elseif niii == "NumeratorCoefficient" || niii == "Numerator"
              push!(num, parse(Float64, content(iii)))

            elseif niii == "DenominatorCoefficient" || niii == "Denominator"
              push!(den, parse(Float64, content(iii)))
            end
          end

        elseif nii == "Decimation"
          resp.fs[ns] = parse(Float64, content(get_elements_by_tagname(ii, "InputSampleRate")[1]))
          resp.fac[ns] = parse(Int64, content(get_elements_by_tagname(ii, "Factor")[1]))
          resp.os[ns] = parse(Int64, content(get_elements_by_tagname(ii, "Offset")[1]))
          resp.delay[ns] = parse(Float64, content(get_elements_by_tagname(ii, "Delay")[1]))
          resp.corr[ns] = parse(Float64, content(get_elements_by_tagname(ii, "Correction")[1]))
        end
      end

      # Post-process
      if resp_code == 0x02
        rmul!(z, c)
        rmul!(p ,c)
        resp.stage[ns] = PZResp64(z = z, p = p, a0 = a0, f0 = f0)
      elseif resp_code == 0x03
        resp.stage[ns] = CoeffResp(num, den)
      else
        resp.stage[ns] = nothing
      end
    end
    # end stages

  end
  L = length(resp.fs)
  for f in fieldnames(MultiStageResp)
    deleteat!(getfield(resp, f), (nmax+1):L)
  end
  return gain, units, resp
end

function FDSN_sta_xml(xmlf::String, msr::Bool, s::String, t::String, v::Integer)

  xdoc = LightXML.parse_string(xmlf)
  xroot = LightXML.root(xdoc)
  xnet = child_elements(xroot)
  S = SeisData()
  s = string_time(s, BUF.date_buf)
  t = string_time(t, BUF.date_buf)

  for net in xnet
    # network level
    if name(net) == "Network"
      net_s = string_time(has_attribute(net, "startDate") ? attribute(net, "startDate") : "0001-01-01T00:00:00", BUF.date_buf)
      net_t = string_time(has_attribute(net, "endDate") ? attribute(net, "endDate") : "9999-12-31T11:59:59", BUF.date_buf)

      # do string comparisons work with <, > in Julia when dates are correctly formatted?
      if net_s ≤ t && net_t ≥ s
        nn = has_attribute(net, "code") ? attribute(net, "code") : ""
        xsta = child_elements(net)

        # station level
        for sta in xsta

          if name(sta) == "Station"
            sta_s = string_time(has_attribute(sta, "startDate") ? attribute(sta, "startDate") : "0001-01-01T00:00:00", BUF.date_buf)
            sta_t = string_time(has_attribute(sta, "endDate") ? attribute(sta, "endDate") : "9999-12-31T11:59:59", BUF.date_buf)

            if sta_s ≤ t && sta_t ≥ s
              ss = has_attribute(sta, "code") ? attribute(sta, "code") : ""

              # loop over child tags
              xcha = child_elements(sta)
              s_lat = 0.0
              s_lon = 0.0
              s_el  = 0.0
              s_name = ""

              for c in xcha
                # station lat
                if name(c) == "Latitude"
                  s_lat = parse(Float64, content(c))

                  # station lon
                elseif name(c) == "Longitude"
                  s_lat = parse(Float64, content(c))

                  # station el
                elseif name(c) == "Elevation"
                  s_el = parse(Float64, content(c))

                # site description -- we only care about the name
                elseif name(c) == "Site"
                  cx = child_elements(c)
                  for i in cx
                    if name(i) == "Name"
                      s_name = content(i)
                    end
                  end

                # channel element
                elseif name(c) == "Channel"
                  cha_s = string_time(has_attribute(c, "startDate") ? attribute(c, "startDate") : "0001-01-01T00:00:00", BUF.date_buf)
                  cha_t = string_time(has_attribute(c, "endDate") ? attribute(c, "endDate") : "9999-12-31T11:59:59", BUF.date_buf)
                  ll = has_attribute(c, "locationCode") ? attribute(c, "locationCode") : ""
                  cc = has_attribute(c, "code") ? attribute(c, "code") : ""

                  if cha_s ≤ t && cha_t ≥ s
                    c_id = "...."
                    c_name = identity(s_name)
                    c_units = ""
                    c_sensor = ""

                    # location defaults
                    c_lat = s_lat
                    c_lon = s_lon
                    c_el  = s_el
                    c_dep = 0.0
                    c_az  = 0.0
                    c_inc = 0.0

                    # sampling and response defaults
                    c_fs        = 0.0
                    c_gain      = 1.0
                    c_normfreq  = 1.0
                    c_drift     = 0.0
                    c_resp      = PZResp()

                    cx = child_elements(c)
                    for i in cx

                      # location params
                      if name(i) == "Latitude"
                        c_lat = parse(Float64, content(i))
                      elseif name(i) == "Longitude"
                        c_lon = parse(Float64, content(i))
                      elseif name(i) == "Elevation"
                        c_el = parse(Float64, content(i))
                      elseif name(i) == "Depth"
                        c_dep = parse(Float64, content(i))
                      elseif name(i) == "Azimuth"
                        c_az = parse(Float64, content(i))
                      elseif name(i) == "Dip"
                        c_inc = parse(Float64, content(i)) - 90.0

                      # sampling, drift
                      elseif name(i) == "SampleRate"
                        c_fs = parse(Float64, content(i))
                      elseif name(i) == "ClockDrift"
                        c_drift = parse(Float64, content(i))

                      # sensor type
                      elseif name(i) == "Sensor"
                        rx = child_elements(i)
                        for r in rx
                          if name(r) == "Description"
                            c_sensor = content(r)
                          end
                        end

                      # response
                      elseif name(i) == "Response"
                        if msr
                          c_gain, c_units, c_resp = full_resp(i)
                        else
                          rx = child_elements(i)
                          for r in rx
                            if name(r) == "InstrumentSensitivity"
                              for y in child_elements(r)
                                if name(y) == "Value"
                                  # gain
                                  c_gain = parse(Float64, content(y))
                                  # println("gain = ", c_gain)
                                elseif name(y) == "Frequency"
                                  # normfreq
                                  # c_normfreq = parse(Float64, content(y))
                                elseif name(y) == "InputUnits"
                                  # calibrationunits
                                  c_units = units2ucum(fix_units(content(get_elements_by_tagname(y, "Name")[1])))
                                end
                              end

                            # stages
                            elseif name(r) == "Stage"
                              for sx in child_elements(r)
                                if name(sx) == "PolesZeros"
                                  p = Array{Complex{Float32},1}(undef, 0)
                                  z = Array{Complex{Float32},1}(undef, 0)
                                  c = 1.0f0
                                  for pzx in child_elements(sx)
                                    npzx = name(pzx)
                                    if npzx == "PzTransferFunctionType"
                                      c = content(pzx) == "LAPLACE (RADIANS/SECOND)" ? c : Float32(2pi)
                                      #= ...this assumes that "DIGITAL (Z-TRANSFORM)"
                                      means what I think it does -- like all terms in
                                      FDSN, there is no documentation =#

                                    # Zero, Pole
                                    elseif npzx in ["Zero", "Pole"]
                                      pzv = complex(
                                        parse(Float32, content(get_elements_by_tagname(pzx, "Real")[1])),
                                        parse(Float32, content(get_elements_by_tagname(pzx, "Imaginary")[1]))
                                      )
                                      if npzx == "Zero"
                                        push!(z, pzv)
                                      else
                                        push!(p, pzv)
                                      end

                                    elseif npzx == "NormalizationFactor"
                                      c_resp.a0 = parse(Float64, content(pzx))

                                    elseif npzx == "NormalizationFrequency"
                                      c_resp.f0 = parse(Float64, content(pzx))

                                    end
                                  end
                                  if length(z) > 0
                                    rmul!(z, c)
                                    append!(c_resp.z, z)
                                  end
                                  if length(p) > 0
                                    rmul!(p, c)
                                    append!(c_resp.p, p)
                                  end
                                end
                              end
                            end
                            # end stages
                          end
                        end
                      end
                      # end response
                    end

                    # channel id
                    c_id = join([nn, ss, ll, cc], ".")
                    if findid(c_id, S) > 0 && (v > 0)
                      @warn(string("Channel ", c_id, " has multiple sets of parameters in time range ", s, " - ", t))
                    end

                    # channel location
                    c_loc = GeoLoc( lat = c_lat,
                                    lon = c_lon,
                                    el  = c_el,
                                    dep = c_dep,
                                    az  = c_az,
                                    inc = c_inc )

                    # build SeisChannel object
                    C = SeisChannel(id    = c_id,
                                    name  = c_name,
                                    gain  = c_gain,
                                    fs    = c_fs,
                                    loc   = c_loc,
                                    units = c_units,
                                    resp  = c_resp )
                    C.misc["ClockDrift"] = c_drift
                    C.misc["startDate"] = cha_s
                    C.misc["endDate"] = cha_t
                    if !isempty(c_sensor)
                      C.misc["SensorDescription"] = c_sensor
                    end
                    push!(S, C)

                  end
                  # done with channel
                end
                # end channel element
              end
            end
          end
        end
        # end station element
      end
    end
    # end network level
  end
  free(xdoc)
  return S
end

function read_sxml(fpat::String, s::String, t::String, memmap::Bool, msr::Bool, v::Integer)

  if safe_isfile(fpat)
    io = memmap ? IOBuffer(Mmap.mmap(fpat)) : open(fpat, "r")
    xsta = read(io, String)
    close(io)
    S = FDSN_sta_xml(xsta, msr, s, t, v)
  else
    files = ls(fpat)
    if length(files) > 0
      io = memmap ? IOBuffer(Mmap.mmap(files[1])) : open(files[1], "r")
      xsta = read(io, String)
      close(io)
      S = FDSN_sta_xml(xsta, msr, s, t, v)
      if length(files) > 1
        for i = 2:length(files)
          io = memmap ? IOBuffer(Mmap.mmap(files[i])) : open(files[i], "r")
          xsta = read(io, String)
          close(io)
          T = FDSN_sta_xml(xsta, msr, s, t, v)
          append!(S, T)
        end
      end
    else
      error("file(s) not found!")
    end
  end
  return S
end

function sxml_mergehdr!(S::GphysData, T::GphysData, app::Bool, nofs::Bool, v::Integer)

  relevant_fields = nofs ? (:name, :loc, :gain, :resp, :units) : (:name, :loc, :fs, :gain, :resp, :units)
  k = Int64[]
  for i = 1:length(T)

    # Match on ID, si, ei
    si = isempty(T.t[i]) ? get(T.misc[i], "startDate", typemin(Int64)) : starttime(T.t[i], T.fs[i])
    ei = isempty(T.t[i]) ? get(T.misc[i], "endDate", typemax(Int64)) : endtime(T.t[i], T.fs[i])
    id = T.id[i]
    c = 0
    for j = 1:length(S.id)
      if S.id[j] == id
        sj = isempty(S.t[j]) ? get(S.misc[j], "startDate", si) : starttime(S.t[j], S.fs[j])
        ej = isempty(S.t[j]) ? get(S.misc[j], "endDate", ei) : endtime(S.t[j], S.fs[j])
        if min(si ≤ ej, ei ≥ sj) == true
          c = j
          break
        end
      end
    end

    # Overwrite S[j] headers on match
    if c != 0
      (v > 2) && println("id/time match! id = ", S.id[c], ". Overwriting S[", c, "] headers from T[", i, "]")
      for f in relevant_fields
        # S.(f)[c] = T.(f)[i]
        setindex!(getfield(S, f), getindex(getfield(T, f), i), c)
      end
      note!(S, c, string("sxml_mergehdr!, overwrote ", relevant_fields))
      S.misc[c] = merge(T.misc[i], S.misc[c])

      # Flag k for delete from T
      push!(k, i)
    end
  end

  # Delete channels that were already used to overwrite headers in S
  if isempty(k) == false
    deleteat!(T, k)
  end

  # Append remainder of T
  (v > 1) && println("remaining entries in T = ", T.n)
  if app && (isempty(T) == false)
    append!(S, T)
  end

  return nothing
end

function read_station_xml!(S::GphysData, file::String, s::String, t::String, msr::Bool, v::Integer)
  if sizeof(file) < 256
    io = open(file, "r")
    xsta = read(io, String)
    close(io)
  else
    xsta = file
  end
  T = FDSN_sta_xml(xsta, msr, s, t, v)
  sxml_mergehdr!(S, T, true, false, v)
  return nothing
end

function read_station_xml(file::String, s::String, t::String, msr::Bool, v::Integer)

  S = SeisData()
  read_station_xml!(S, file, s, t, msr, v)
  return S
end

function msr_to_xml(io::IO, r::MultiStageResp, gain::Float64, units::String)
  Nstg = length(r.stage)
  sens_f = 0.0
  if Nstg > 0
    sens_f = r.fg[1]
  end

  write(io, "          <InstrumentSensitivity>\n            <Value>")
  print(io, gain)
  write(io, "</Value>\n            <Frequency>")
  print(io, sens_f)
  write(io, "</Frequency>\n            <InputUnits>\n              <Name>")
  write(io, get(ucum_to_seed, units, units))
  write(io, "</Name>\n            </InputUnits>\n          </InstrumentSensitivity>\n")

  u1 = "              <InputUnits>\n                <Name>"
  u2 = "</Name>\n              </InputUnits>\n              <OutputUnits>\n                <Name>"
  u3 =   "</Name>\n              </OutputUnits>\n"

  @inbounds for i in 1:Nstg
    ui = get(ucum_to_seed, r.i[i], r.i[i])
    uo = get(ucum_to_seed, r.o[i], r.o[i])

    write(io, "          <Stage number=\"")
    print(io, i)
    write(io, "\">\n")

    if typeof(r.stage[i]) in (PZResp64, PZResp)
      write(io, "            <PolesZeros>\n")
      write(io, u1)
      write(io, ui)
      write(io, u2)
      write(io, uo)
      write(io, u3)
      write(io, "              <PzTransferFunctionType>LAPLACE (RADIANS/SECOND)</PzTransferFunctionType>\n              <NormalizationFactor>")
      print(io, r.stage[i].a0)
      write(io, "</NormalizationFactor>\n              <NormalizationFrequency>")
      print(io, r.stage[i].f0)
      write(io, "</NormalizationFrequency>\n")
      for j = 1:length(r.stage[i].z)
        write(io, "              <Zero number = \"")
        print(io, j)
        write(io, "\">\n                <Real minusError=\"0\" plusError=\"0\">")
        print(io, real(r.stage[i].z[j]))
        write(io, "</Real>\n                <Imaginary minusError=\"0\" plusError=\"0\">")
        print(io, imag(r.stage[i].z[j]))
        write(io, "</Imaginary>\n              </Zero>\n")
      end
      for j = 1:length(r.stage[i].p)
        write(io, "              <Pole number = \"")
        print(io, j)
        write(io, "\">\n                <Real minusError=\"0\" plusError=\"0\">")
        print(io, real(r.stage[i].p[j]))
        write(io, "</Real>\n                <Imaginary minusError=\"0\" plusError=\"0\">")
        print(io, imag(r.stage[i].p[j]))
        write(io, "</Imaginary>\n              </Pole>\n")
      end
      write(io, "            </PolesZeros>\n")
    else
      if typeof(r.stage[i]) == CoeffResp
        write(io, "            <Coefficients>\n")
        write(io, u1)
        write(io, ui)
        write(io, u2)
        write(io, uo)
        write(io, u3)
        write(io, "              <CfTransferFunctionType>DIGITAL</CfTransferFunctionType>\n")
        for j = 1:length(r.stage[i].b)
          write(io, "              <Numerator minusError=\"0\" plusError=\"0\">")
          print(io, r.stage[i].b[j])
          write(io, "</Numerator>\n")
        end
        for j = 1:length(r.stage[i].a)
          write(io, "              <Denominator minusError=\"0\" plusError=\"0\">")
          print(io, r.stage[i].a[j])
          write(io, "</Denominator>\n")
        end
        write(io, "            </Coefficients>\n")
      end

      if r.fac[i] > 0
        write(io, "            <Decimation>\n              <InputSampleRate>")
        print(io, r.fs[i])
        write(io, "</InputSampleRate>\n              <Factor>")
        print(io, r.fac[i])
        write(io, "</Factor>\n              <Offset>")
        print(io, r.os[i])
        write(io, "</Offset>\n              <Delay>")
        print(io, r.delay[i])
        write(io, "</Delay>\n              <Correction>")
        print(io, r.corr[i])
        write(io, "</Correction>\n            </Decimation>\n")
      end
    end
    write(io, "            <StageGain>\n              <Value>")
    print(io, r.gain[i])
    write(io, "</Value>\n              <Frequency>")
    print(io, r.fg[i])
    write(io, "</Frequency>\n            </StageGain>\n          </Stage>\n")
  end
  return nothing
end

function mk_xml!(io::IO, S::GphysData, chans::Array{Int64,1})
  blank_t0 = round(Int64, d2u(now())*μs)

  write(io, "<?xml version=\"1.0\" encoding=\"ISO-8859-1\"?>\n\n<FDSNStationXML xmlns=\"http://www.fdsn.org/xml/station/1\" schemaVersion=\"1.0\" xsi:schemaLocation=\"http://www.fdsn.org/xml/station/1 http://www.fdsn.org/xml/station/fdsn-station-1.0.xsd\" xmlns:xsi=\"http://www.w3.org/2001/XMLSchema-instance\">\n  <Source>SeisIO</Source>\n\n  <Created>")
  print(io, now())
  write(io, "</Created>\n")

  nch = length(chans)
  nets = Array{String,1}(undef, nch)
  stas = Array{String,1}(undef, nch)
  locs = Array{String,1}(undef, nch)
  chas = Array{String,1}(undef, nch)
  t0 = zeros(Int64, nch)
  t1 = zeros(Int64, nch)

  done = falses(nch)

  # Fill in nets, stas, locs, chas
  fill!(t0, blank_t0)
  fill!(t1, xml_endtime)
  @inbounds for j = 1:nch
    i = chans[j]
    id = split_id(S.id[i])
    nets[j] = id[1]
    stas[j] = id[2]
    locs[j] = id[3]
    chas[j] = id[4]

    # precedence for start time: S.misc[i]["startDate"] > starttime(S.t[i], S.fs[i]) > blank_t0
    if haskey(S.misc[i], "startDate")
      c_start = get(S.misc[i], "startDate", blank_t0)
      if typeof(c_start) == Int64
        t0[j] = c_start
      end
    end

    # precedence for end time: S.misc[i]["startDate"] > endtime(S.t[i],S.fs[i]) > 19880899199000000
    if haskey(S.misc[i], "endDate")
      c_end = get(S.misc[i], "endDate", xml_endtime)
      if typeof(c_end) == Int64
        t1[j] = c_end
      end
    end

    if isempty(S.t[i]) == false
      if t0[j] == blank_t0
        t0[j] = starttime(S.t[i], S.fs[i])
      end
      if t1[j] == xml_endtime
        t1[j] = endtime(S.t[i], S.fs[i])
      end
    end
  end

  # Loop over all channels in cha that are not done
  @inbounds for j = 1:nch
    (done[j]) && continue

    # Get all selected channels from same network
    nn = nets[j]
    ss = stas[j]
    net = findall(nets.==nn)
    ts = minimum(t0[net])*μs
    te = maximum(t1[net])*μs
    ds = u2d(ts)
    de = u2d(te)

    # Open net node
    write(io, "  <Network code=\"")
    write(io, nn)
    write(io, "\" startDate=\"")
    print(io, ds)
    write(io, "\" endDate=\"")
    print(io, de)
    write(io, "\">\n")

    # Loop over stations in network
    ci = Int64[]
    sta_name = ""
    sta_lat = 0.0
    sta_lon = 0.0
    sta_el = 0.0

    # Fill channel index array, ci; set station location and name
    for j = 1:nch
      if stas[j] == ss && nets[j] == nn
        i = chans[j]
        push!(ci, j)

        if isempty(sta_name)
          if !isempty(S.name[i])
            sta_name = S.name[i]
          end
        end

        if (sta_lat == 0.0) || (sta_lon == 0.0) || (sta_el == 0.0)
          if typeof(S.loc[i]) == GeoLoc
            loc = S.loc[i]
            if loc.lat != 0.0
              sta_lat = loc.lat
            end

            if loc.lon != 0.0
              sta_lon = loc.lon
            end

            if loc.el != 0.0
              sta_el = loc.el
            end
          end
        end
      end
    end
    ts = minimum(t0[ci])*μs
    te = maximum(t1[ci])*μs
    ds = u2d(ts)
    de = u2d(te)

    # Open sta node
    write(io, "    <Station code=\"")
    write(io, ss)
    write(io, "\" startDate=\"")
    print(io, ds)
    write(io, "\" endDate=\"")
    print(io, de)
    write(io, "\">\n      <Latitude>")
    print(io, sta_lat)
    write(io, "</Latitude>\n      <Longitude>")
    print(io, sta_lon)
    write(io, "</Longitude>\n      <Elevation>")
    print(io, sta_el)
    write(io, "</Elevation>\n      <Site>\n        <Name>")
    write(io, sta_name)
    write(io, "</Name>\n      </Site>\n")

    # Loop over cha nodes
    for j in ci
      i = chans[j]

      # Channel location
      if typeof(S.loc[i]) == GeoLoc
        loc = S.loc[i]
      else
        loc = GeoLoc()
      end

      # Instrument response
      if typeof(S.resp[i]) == MultiStageResp
        resp = S.resp[i]
      else
        if typeof(S.resp[i]) in (PZResp, PZResp64)
          resp = MultiStageResp(2)
          resp.fac[2] = 1
        else
          resp = MultiStageResp(1)
        end
        resp.stage[1] = S.resp[i]
      end

      # Open cha node
      write(io, "      <Channel code=\"")
      write(io, chas[j])
      write(io, "\" locationCode=\"")
      write(io, locs[j])
      write(io, "\" startDate=\"")
      print(io, u2d(t0[j]*μs))
      write(io, "\" endDate=\"")
      print(io, u2d(t1[j]*μs))
      write(io, "\">\n", "        <Latitude>")
      print(io, loc.lat)
      write(io, "</Latitude>\n        <Longitude>")
      print(io, loc.lon)
      write(io, "</Longitude>\n        <Elevation>")
      print(io, loc.el)
      write(io, "</Elevation>\n        <Depth>")
      print(io, loc.dep)
      write(io, "</Depth>\n        <Azimuth>")
      print(io, loc.az)
      write(io, "</Azimuth>\n        <Dip>")
      print(io, 90.0 + loc.inc)
      write(io, "</Dip>\n        <SampleRate>")
      print(io, S.fs[i])
      write(io, "</SampleRate>\n        <ClockDrift>")
      print(io, get(S.misc[i], "ClockDrift", 0.0))
      write(io, "</ClockDrift>\n        <Sensor>\n          <Description>")
      write(io, get(S.misc[i], "SensorDescription", "Unknown"))
      write(io, "</Description>\n        </Sensor>\n        <Response>\n")
      msr_to_xml(io, resp, S.gain[i], S.units[i])
      write(io, "        </Response>\n      </Channel>\n")

      done[j] = true
    end
    write(io, "    </Station>\n  </Network>\n")
  end
  write(io, "</FDSNStationXML>\n")

  return nothing
end

"""
    write_sxml(fname::String, S::GphysData[, chans=Cha])
    write_sxml(fname::String, C::GphysChannel)

Write station XML from the fields of `S` or `C` to file `fname`.

Use keyword `chans=Cha` to restrict station XML write to `Cha`. This keyword
can accept an Integer, UnitRange, or Array{Int64,1} as its argument.
"""
function write_sxml(fname::String, S::GphysData;
  chans::ChanSpec=Int64[])

  chans = mkchans(chans, S, keepempty=true)
  io = open(fname, "w")
  mk_xml!(io, S, chans)
  close(io)

  for i in chans
    fwrite_note!(S, i, "write_sxml", fname, string(", chans=", i))
  end
  return nothing
end
write_sxml(fname::String, C::GphysChannel) = write_sxml(fname, SeisData(C), chans=1)


text:
module OPMC
#
# Monte-Carlo type methods for stochastic dynamical systems.
#

export matchPoints, poissonDisk, selectPoints, resample

using Random

using ..OPccmt

# ===================================================================
function matchPoints(xv, xz)
#
# Representing a collection of particles as a chaos expansion
# requires mapping the coordinates from the physical domain to a
# reference domain.  Given the coordinates of the physical points,
# and a set of reference coordinates, map the physical to the
# reference coordinates.
#
# The algorithm used is "farthest nearest neighbor".  It is a greedy
# sequential algorithm assigning closest point pairs; but it does so
# in an order starting from the point whose minimum distance among
# all candidate matching points is greatest.  This prevents the case
# where all the close points are assigned first, leaving the last
# couple matches between points on opposite sides of the domain.
#
# Here we will call the points in physical space "v" and in reference
# space "z".
#
# Version:        Changes:
# --------        -------------
# 08.07.2021      Original code.
#
# Version:        Verification:
# --------        -------------
# 08.07.2021      Checked a 2D example problem.
#
# Inputs:
# -------
# xv              : A D-by-N matrix of points in physical space,
#                   where D is the number of dimensions and N the
#                   number of points.
# xz              : The matching points in the reference space, also
#                   a D-by-N matrix.
#
# Outputs:
# --------
# k               : A 2-by-N matrix of indices.  The first row is
#                   the index of a v point, and the second of a z
#                   point, such that each column identifies a pair.

   D = size(xv,1)  # Number of dimensions.
   N = size(xv,2)  # Number of points.

   # Compute the coordinate-wise cumulative distributions.  Future
   # improvement: this can likely be accelerated by the use of a
   # 1D sorting algorithm.
   Pv = zeros(D,N)
   Pz = zeros(D,N)
   for ip in 1:N
      for id in 1:D
         ind = (xv[id,:] .< xv[id,ip])
         Pv[id,ip] = sum(ind)
         ind = (xz[id,:] .< xz[id,ip])   # Future improvement: can compute
         Pz[id,ip] = sum(ind)            # this once upfront.
      end
   end

   # Build the matrix of distances between Pv and Pz points.
   Svz  = zeros(N,N)
   for ir in 1:N
      for ic in 1:N
         for id in 1:D
            Svz[ir,ic] = Svz[ir,ic] + (Pv[id,ir] - Pz[id,ic])^2
         end
         Svz[ir,ic] = sqrt(Svz[ir,ic])
      end
   end

   rng = range(1,stop=N)
   indr = ones(Bool,N)
   indc = ones(Bool,N)
   k = zeros(Int,2,N)
   for ip in 1:N
      nn = N+1-ip
      minr = zeros(nn)
      minc = zeros(nn)
      rngr = rng[indr]
      rngc = rng[indc]
      for jp in 1:nn
         minr[jp] = minimum(Svz[rngr[jp],rngc])
         minc[jp] = minimum(Svz[rngr,rngc[jp]])
      end
      irmax = argmax(minr)
      icmax = argmax(minc)
      rmax = rngr[irmax]
      cmax = rngc[icmax]

      if (minr[irmax] > minc[icmax])
         imin = argmin(Svz[rmax,rngc])
         cmin = rngc[imin]
         k[1,ip] = rmax
         k[2,ip] = cmin
         indr[rmax] = false
         indc[cmin] = false
      else
         imin = argmin(Svz[rngr,cmax])
         rmin = rngr[imin]
         k[1,ip] = rmin
         k[2,ip] = cmax
         indr[rmin] = false
         indc[cmax] = false
      end
   end

   return k

end  # function matchPoints.

# ===================================================================
function poissonDisk(D, N, k; rseed=-1)
#
# Fills space with a set of points that are randomly located, but a
# minimum distance from all other points in the domain.  This is
# known as "Poisson disk sampling", or "blue noise sampling".
#
# The domain is the D-dimensional unit cube.  The coordinates can be
# easily scaled afterwards.
#
# Bridson R.  Fast Poisson disk sampling in arbitrary dimensions.
# Proceedings of ACM SIGGRAPH 2007, p 22.
#
# Roberts M.  http://extremelearning.com.au/an-improved-version-of-
# bridsons-algorithm-n-for-poisson-disc-sampling/
#
# Version:        Changes:
# --------        -------------
# 26.07.2021      Original code.
#
# Version:        Verification:
# --------        -------------
# 26.07.2021      Plotted 2D cases, showing proper execution.  Ran
#                 3D and 4D cases.
#
# Inputs:
# -------
# D               : The dimension of the space.
# N               : The number of cells along each dimension.
# k               : The maximum number of trials before rejection.
#
# Outputs:
# --------
# x               : Coordinates of the points.

   if (rseed > 0)
      Random.seed!(rseed)
   end

   del = sqrt(eps())
   sqD = sqrt(D)

   # Initialize a background grid.
   M = N^D               # Number of cells.
   dx = 1.0/N
   r  = dx*sqrt(D)

   xb = zeros(D,3)
   xb[:,1] .= 0.5*dx
   xb[:,2] .= dx
   xb[:,3] .= 1.0 - 0.5*dx

   ic = range(1,stop=M)
   xc = xFromCell(ic,xb)

   # Dynamically allocate the number of points.
   Narray = 1

   act = zeros(Bool,Narray)     # Active point flag.
   xs  = zeros(D,Narray)        # Point coordinates.
   cx  = zeros(Int,Narray)      # Point cell association.
   ind = range(1, stop=Narray)  # [1:Narray].

   # Pick a starting coordinate at random.
   xs[:,1] = rand(D)
   jc, wc = nearestCells(xs[:,1],xb)
   cx[1] = jc[argmax(wc)]
   act[1] = true

   np = 1

   flg = true
   while (flg)

      # Pick a random one of the active points.
      iact = ind[act]
      nact = length(iact)

      if (nact == 0)
         flg = false
      else

         jp   = round(Int,(-0.5 + eps()) + nact*(rand() - eps()) + 1)
         ip   = iact[jp]

         # Sample k points about the hypersphere in sequence, looking for
         # one that is at least a distance r from existing points.
         kflg = true
         kk = 0
         while (kflg)

            kk = kk + 1

            vec = -1.0 .+ 2.0*rand(D)
            lv = 0.0
            for jj in 1:D
               lv = lv + vec[jj]^2
            end
            n = vec/sqrt(lv)

            xtry = xs[:,ip] + (r + del) .* n

            # Wrap around the domain.
            for jj in 1:D
               if (xtry[jj] < 0.0)
                  xtry[jj] = xtry[jj] + 1.0
               elseif (xtry[jj] > 1.0)
                  xtry[jj] = xtry[jj] - 1.0
               end
            end

            # Compute the minimum distance to other points.
            jc, wc = nearestCells(xtry,xb)
            isur = surroundingCells(xtry,xb)
            Ncell = length(isur)
            dmin = 2.0*r
            for icell in 1:Ncell
               jcx = (cx .== isur[icell])
               if (sum(jcx) >= 1)  # There is a point associated with this cell.
                  cind = argmax(jcx)
                  dist = 0.0
                  for jj in 1:D
                     dist = dist + (xs[jj,cind] - xtry[jj])^2
                  end
                  dist = sqrt(dist)
                  if (dist < dmin)
                     dmin = dist
                  end
               end
            end

            if (dmin >= r)

               # Resize the arrays if needed.
               if (Narray == np)
                  Narray = 2*np   # Double on each resizing.

                  temp = deepcopy(act)
                  act = zeros(Bool,Narray)
                  act[1:np] = temp

                  temp = deepcopy(xs)
                  xs = zeros(D,Narray)
                  xs[:,1:np] = temp

                  temp = deepcopy(cx)
                  cx = zeros(Int,Narray)
                  cx[1:np] = temp

                  ind = range(1, stop=Narray)

               end

               # Success!  Log the point.
               np = np + 1
               xs[:,np] = xtry
               cx[np] = jc[argmax(wc)]
               act[np] = true
               kflg = false

            elseif (kk == k)

               act[ip] = false  # No available points around ip.
               kflg = false

            end

         end  # while (kflg)

#print("nact = ",nact,"  jp = ",jp,"  ip = ",ip,"  np = ",np,"\n")

      end  # if (nact != 0)

   end  # while (flg)

   return xs[:,1:np]

end  # function poissonDisk.

# ===================================================================
function selectPoints(xin, N)
#
# It is not possible to specify an exact number of points generated
# by the poissonDisk function.  The best that can be done is to
# generate a few too many points, and then take a subset of the
# space containing the desired number of points, rescaling to fill
# the unit cube.
#
# This works by performing a bisecting search along the diagonal
# from (0,0,...,0) to (1,1,...,1) until the proper number of points
# is encapsulated.  Then the coordinates are scaled up to the unit
# cube.
#
# Version:        Changes:
# --------        -------------
# 28.07.2021      Original code.
#
# Version:        Verification:
# --------        -------------
# 28.07.2021      Tested a 2D case.
#
# Inputs:
# -------
# xin             : Point coordinates defined on the unit cube.
#                   Array of dimension D-by-M.
# N               : Desired number of points on the unit cube.  Must
#                   be less than M.
#
# Outputs:
# --------
# x               : Coordinates of the N points normalized to fill
#                   the unit cube.

   D = size(xin,1)
   M = size(xin,2)

   # Declare variables that would otherwise be undefined outside the
   # "for" loop.  [Unintuitive Julia name scoping conventions...]
   ind = ones(Bool,M)
   a = 0.0

   if (N >= M)

      return xin

   else

      aLB = 0.0
      aUB = 1.0

      flg = true
      while (flg)

         # Find the bisecting point.
         a = 0.5*(aLB + aUB)

         # Find the number of points below the bisecting point.
         ind = ones(Bool,M)
         for jj in 1:D
            ind = (ind .& (xin[jj,:] .<= a))
         end
         nb = sum(ind)

         if (nb == N)
            flg = false
         elseif (nb > N)
            # Move down to eliminate points.
            aUB = a
         else  # nb < N
            # Move up to capture more points.
            aLB = a
         end

      end  # while (flg)

      return (xin[:,ind] * (1.0/a))

   end  # if (N >= M)

end  # selectPoints

# ===================================================================
function resample(xi, pri; Msam = 0)
#
# Based on the "systematic resampling" algorithm of Arulampalam
# et al. A tutorial on particle filters for online nonlinear/
# non-Gaussian Bayesian tracking. IEEE Transactions on Signal
# Processing 50: 174-188.
#
# Version:        Changes:
# --------        -------------
# 16.08.2021      Original code.
#
# Version:        Verification:
# --------        -------------
# 16.08.2021      Checked using some straightforward cases.
#
# Inputs:
# -------
# xi              : D-by-N array of points.
# pri             : Input probabilities.
# Msam            : (Optional) Number of points to sample.
#
# Outputs:
# --------
# xo              : D-by-M array of resampled points.
# pro             : Output probabilities.

   D = size(xi,1)
   N = size(xi,2)

   if (Msam == 0)
      M = N
   else
      M = Msam
   end

   C = zeros(N)
   C[1] = pri[1]
   for ii in 2:N-1
      C[ii] = C[ii-1] + pri[ii]
   end
   C[N] = 1.0 + eps()

   Mm1 = 1.0/M

   u1 = rand()*Mm1

   xo = zeros(D,M)
   pro = zeros(M)
   io = zeros(Int,M)
   ii = 1
   for jj in 1:M
      u = u1 + (jj-1)*Mm1
      while (u > C[ii])
         ii = ii + 1
      end
      xo[:,jj] = xi[:,ii]
      pro[jj] = Mm1
      io[jj] = ii
   end

   return xo, pro, io

end  # function resample

end  # module OPMC.

max_stars_repo_name: SINTEF-Energy-Wind/OutliersPrimer

max_issues_repo_name: SINTEF-Energy-Wind/OutliersPrimer

max_forks_repo_name: SINTEF-Energy-Wind/OutliersPrimer


text:
"""
    struct DeclarationName <: Any
Represent a declaration name.

Note that, the underlying pointer is NOT a *pointer* to a `clang::DeclarationName` object.
Instead, it's the opaque pointer representation of the `clang::DeclarationName` itself.
"""
struct DeclarationName
    ptr::CXDeclarationName
end

Base.unsafe_convert(::Type{CXDeclarationName}, x::DeclarationName) = x.ptr
Base.cconvert(::Type{CXDeclarationName}, x::DeclarationName) = x


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     comment_magics: false
#     formats: ipynb,jl:light
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.4.2
#   kernelspec:
#     display_name: Julia 1.4.0
#     language: julia
#     name: julia-1.4
# ---

# ## Initialize a particle group

using GEMPIC, Plots, Test

# +
function test_sampling( sampling_type :: Symbol,
                        symmetric     :: Bool,
                        pg            :: ParticleGroup{D,V},
                        df             ) where {D, V}

   mean  = zeros(3)
   sigma = zeros(3)
   xi    = zeros(3)

   n_particles = pg.n_particles

   sampling = ParticleSampler{D,V}( sampling_type, symmetric, n_particles )

   GEMPIC.sample!( pg, sampling, df, mesh )

   for i_part = 1:n_particles
       xi = GEMPIC.get_x(pg, i_part)
       vi = GEMPIC.get_v(pg, i_part)
       mean[1] += xi[1]
       mean[2] += vi[1]
       mean[3] += vi[2]
   end

   mean = mean/n_particles

   for i_part = 1:n_particles
       xi = GEMPIC.get_x(pg, i_part)
       vi = GEMPIC.get_v(pg, i_part)
       sigma[1] += (xi[1] - mean[1])^2
       sigma[2] += (vi[1] - mean[2])^2
       sigma[3] += (vi[2] - mean[3])^2
   end

   sigma = sigma/(n_particles-1)

   mean, sigma

end
# -

n_particles = 100000
xmin        = 1.0 :: Float64
xmax        = 4π + 1.0
Lx          = xmax - xmin
nx          = 64

mesh = Mesh( xmax, xmin, nx)

pg = ParticleGroup{1,2}(n_particles)

params = ( k = [[0.5]],
           α = [0.01],
           σ = [[0.1, 2.0]],
           μ = [[0.0, 0.0]]
)

df1 = CosSumGaussian{1,2}(params...)

mean_ref  = [Lx*0.5+xmin, 0.0, 0.0]
sigma_ref = [Lx^2/12.0,   params.σ[1][1]^2, params.σ[1][2]^2 ]

@info "Sobol non symmetric"
mean, sigma = test_sampling( :sobol, false, pg, df1 )
@show abs.(mean  .- mean_ref)
@show abs.(sigma .- sigma_ref)
@test maximum(abs.(mean .- mean_ref)) ≈ 0.0 atol = 1e2/sqrt(n_particles)

@info "Sobol symmetric"
mean, sigma = test_sampling( :sobol, true,  pg, df1 )
@show abs.(mean  .- mean_ref)
@show abs.(sigma .- sigma_ref)
@test maximum(abs.(mean .- mean_ref)) ≈ 0.0 atol = 1e-12

@info "Random non symmetric"
mean, sigma = test_sampling( :random, false, pg, df1)
@show abs.(mean  .- mean_ref)
@show abs.(sigma .- sigma_ref)
@test maximum(abs.(mean .- mean_ref)) ≈ 0.0 atol = 1e2/sqrt(n_particles)

@info "Random symmetric"
mean, sigma = test_sampling( :random, true,  pg, df1)
@show abs.(mean  .- mean_ref)
@show abs.(sigma .- sigma_ref)
@test maximum(abs.(mean .- mean_ref)) ≈ 0.0 atol = 1e-12

xp = view(pg.array,1,:)
wp = view(pg.array,4,:);

histogram(xp, weights=wp, normalize=true, bins=200)

# Expected mean:
# 2π+1, 0, 0
# Expected variance:
# 16/12 π^2 , 0.01, 4

params = (
    k = [[0.5]],
    α = [0.01],
    σ = [[0.1, 2.0], [2.0, 2.0]],
    μ = [[0.0, 0.0], [1.0, 1.0]],
    δ = [0.7, 0.3]
)

df2 = CosSumGaussian{1,2}(params...)

mean_ref = [Lx*0.5+xmin, 0.3, 0.3]
sigma_ref[1] = Lx^2/12.0
for j=1:2
    sigma_ref[j+1] = - (df2.params.δ[1] * df2.params.μ[j][1]
                     + (df2.params.δ[2]
                     -  df2.params.δ[1]) * df2.params.μ[j][2])^2
    for k=1:2
        sigma_ref[j+1] += df2.params.δ[k] * (df2.params.σ[j][k]^2 + df2.params.μ[j][k]^2)
    end
end

@info "Sobol non symmetric"
mean, sigma =  test_sampling( :sobol, false, pg, df2)
@show abs.(mean  .- mean_ref)
@show abs.(sigma .- sigma_ref)
@test maximum(abs.(mean .- mean_ref)) ≈ 0.0 atol = 1e2/sqrt(n_particles)

@info "Sobol symmetric"
mean, sigma =  test_sampling( :sobol, true,  pg, df2)
@show abs.(mean  .- mean_ref)
@show abs.(sigma .- sigma_ref)
@test maximum(abs.(mean .- mean_ref)) ≈ 0.0 atol = 1e2/sqrt(n_particles)


text:
include("../../matplotlib_settings.jl")

set_config(:tab_compensated_summation, false)

using ChargedParticleDynamics.GuidingCenter4d: plot_integral_error, plot_poincare_surface

const coll_title = "Guiding Center 4D Fast 2nd Poincaré Integral Invariant"

XMIN = 1.5
XMAX = 2.0
YMIN = -0.25
YMAX = +0.25

include("../guiding_center_4d_poincare_invariant_2nd.jl")


text:
function dist_matrix(G::SimpleDigraph)
    VV = vlist(G)
    n = length(VV)
    A = zeros(Int, n, n)
    for u = 1:n
        d = dist(G, VV[u])
        for v = 1:n
            A[u, v] = d[VV[v]]
        end
    end
    return A
end

function diam(G::SimpleDigraph)
    A = dist_matrix(G)
    if minimum(A) < 0
        return -1
    end
    return maximum(A)
end


text:
using Optim
using Devectorize

"""
Linear Regression Cost function with Regularisation

Parameters:
	params      = vector containing rolled up X and Θ
		X       = array of features for each movie
		Θ       = array of weights for each user
	Y           = array of ratings for each movie
	R           = array of the existance of ratings for each movie
	numUsers    = number of users
	numMovies   = number of movies
	numFeatures = number of features
	λ           = regularisation constant

Returns:
	J = linear regression cost of the given weights with the given regularisation constant
"""
function cost(params, Y::Array{Int}, numUsers, numMovies, numFeatures, λ)
	# Unfold params -> X, Θ
	X, Θ = unroll(params, numUsers, numMovies, numFeatures)

	# Calculate cost
	J = 0
 	@simd for i = 1:size(Y)[2]
		@inbounds J += (dot(X[Y[2, i], :][:], Θ[Y[1, i], :][:]) - Y[3, i])^2
	end
	J = 1/2 * J
	# add regularisation
	J += λ/2 * (sum(Θ.^2) + sum(X.^2))

	return J
end

"""
Linear Regression Gradient with Regularisation

Parameters:
	params      = vector containing rolled up X and Θ
		X       = array of features for each movie
		Θ       = array of weights for each user
	gradient    = vector containing rolled up Xgrad and ΘGrad from previous iteration
	Y           = array of ratings for each movie
	R           = array of the existance of ratings for each movie
	numUsers    = number of users
	numMovies   = number of movies
	numFeatures = number of features
	λ           = regularisation constant

Modifies:
	gradient = XGrad and ΘGrad are updated based on the given weights and regularisation costant
"""
function grad!(params, gradient, Y, numUsers, numMovies, numFeatures, λ)
	@printf "hi"
	# Unfold params -> X, Θ
	X, Θ = unroll(params, numUsers, numMovies, numFeatures)

	XGrad = zeros(size(X))
	ΘGrad = zeros(size(Θ))

	for n = 1:size(Y)[2]
		i = Y[2, n]
		j = Y[1, n]
		@devec x = X[i, :]
		@devec þ = Θ[j, :]
		temp = dot(x[:], þ[:]) - Y[3, n]
		@devec XGrad[i, :] += temp .* þ + λ.*x
		@devec ΘGrad[j, :] += temp .* x + λ.*þ
	end

	gradient[:] = rollup(XGrad, ΘGrad)
end

function costGrad!(params, gradient, Y, numUsers, numMovies, numFeatures, λ)
	# Unfold params -> X, Θ
	X, Θ = unroll(params, numUsers, numMovies, numFeatures)

	XGrad = zeros(size(X))
	ΘGrad = zeros(size(Θ))

	# Calculate cost
	J = 0
	count = Int(floor(size(Y)[2] * 0.1))
	start = rand(1:size(Y)[2] - count)
	for n = start:start+count
	#for n = 1:size(Y)[2]
		i = Y[2, n]
		j = Y[1, n]
		@devec x = X[i, :]
		@devec þ = Θ[j, :]
		temp = dot(x[:], þ[:]) - Y[3, n]
		J += temp^2
		@devec XGrad[i, :] += temp .* þ + λ.*x
		@devec ΘGrad[j, :] += temp .* x + λ.*þ
	end
	J = 1/2 * J
	# add regularisation
	J += λ/2 * (sum(Θ.^2) + sum(X.^2))

	gradient[:] = rollup(XGrad, ΘGrad)

	return J
end

"""
Unroll parameter vector into X and Θ

Parameters:
	params = vector containing rolled up X and Θ
	numUsers    = number of users
	numMovies   = number of movies
	numFeatures = number of features

Returns:
	X = array of features for each movie
	Θ = array of weights for each user
"""
function unroll(params, numUsers, numMovies, numFeatures)
	X = reshape(params[1 : numMovies * numFeatures], numMovies, numFeatures)
	Θ = reshape(params[numMovies * numFeatures + 1 : end], numUsers, numFeatures)
	return X, Θ
end

"""
Roll up X and Θ into single vector

Parameters:
	X = array of features for each movie
	Θ = array of weights for each user

Returns:
	params = vector containing rolled up X and Θ
"""
function rollup(X, Θ)
	return params = vcat(X..., Θ...)
end

"""
Entry point

Params:
	dataPath = path to() CSVs
"""
function main()
	dataPath = ""
	if ARGS[1] == "small"
		dataPath = "../ml-latest-small/fixed/"
	elseif ARGS[1] == "full"
		dataPath = "../ml-latest/fixed/"
	end

	# Load data from CSV - using wc to speed up loading by preallocating
	# N.B. requires cat and wc to be installed - a shell should not be nessecary though
	lines = parse(readall(pipeline(`cat $[dataPath]ratings.csv`, `wc -l`)))
	@time ratings = readdlm(dataPath * "ratings.csv", ',', header=true, dims=(lines, 4))[1]

	lines = parse(readall(pipeline(`cat $[dataPath]train.csv`, `wc -l`)))
	@time train = readdlm(dataPath * "train.csv", ',', header=true, dims=(lines, 4))[1]

	lines = parse(readall(pipeline(`cat $[dataPath]cross.csv`, `wc -l`)))
	@time cross = readdlm(dataPath * "cross.csv", ',', header=true, dims=(lines, 4))[1]

	lines = parse(readall(pipeline(`cat $[dataPath]movies.csv`, `wc -l`)))
	@time movies = readdlm(dataPath * "movies.csv", ',', header=true, dims=(lines, 3))
	moviesHeader = movies[2]
	movies = movies[1]

	@printf "Finished loading, begin processing\n"

	train[:, 3] *= 2
	train = Array{Int}(train[:, 1:3]')
	cross[:, 3] *= 2
	cross = Array{Int}(cross[:, 1:3]')

	# Count users, movies and features
	numUsers = length(union(ratings[:,1]))
	numMovies = size(movies)[1]

	numFeatures = 18
	λ = 0.05
	X = rand(numMovies, 18)
	Θ = rand(numUsers, 18)

	if numFeatures == 18
		for i = 1:numMovies
			X[i, 1] = contains(movies[i, 3], "Action") ? 2 : 0
			X[i, 2] = contains(movies[i, 3], "Adventure") ? 2 : 0
			X[i, 3] = contains(movies[i, 3], "Animation") ? 2 : 0
			X[i, 4] = contains(movies[i, 3], "Children") ? 2 : 0
			X[i, 5] = contains(movies[i, 3], "Comedy") ? 2 : 0
			X[i, 6] = contains(movies[i, 3], "Crime") ? 2 : 0
			X[i, 7] = contains(movies[i, 3], "Documentary") ? 2 : 0
			X[i, 8] = contains(movies[i, 3], "Drama") ? 2 : 0
			X[i, 9] = contains(movies[i, 3], "Fantasy") ? 2 : 0
			X[i, 10] = contains(movies[i, 3], "Film-Noir") ? 2 : 0
			X[i, 11] = contains(movies[i, 3], "Horror") ? 2 : 0
			X[i, 12] = contains(movies[i, 3], "Musical") ? 2 : 0
			X[i, 13] = contains(movies[i, 3], "Mystery") ? 2 : 0
			X[i, 14] = contains(movies[i, 3], "Romance") ? 2 : 0
			X[i, 15] = contains(movies[i, 3], "Sci-Fi") ? 2 : 0
			X[i, 16] = contains(movies[i, 3], "Thriller") ? 2 : 0
			X[i, 17] = contains(movies[i, 3], "War") ? 2 : 0
			X[i, 18] = contains(movies[i, 3], "Western") ? 2 : 0
		end
	end
	scoreboard = Array{Any}(0,3)
	for i = 0.1:0.10:1.00
		tempTrain = train[:, 1:Int(floor(i*size(train)[2]))]

		function costTemp(x)
			cost(x, train, numUsers, numMovies, numFeatures, λ)
		end

		function gradTemp!(x, y)
			grad!(x, y, train, numUsers, numMovies, numFeatures, λ)
		end

		function costGradTemp!(x, y)
			costGrad!(x, y, tempTrain, numUsers, numMovies, numFeatures, λ)
		end

		tic()
		@time @show results = optimize(DifferentiableFunction(costTemp, gradTemp!, costGradTemp!), rollup(X,Θ), iterations=10)
		time = toc()

		@show score = cost(results.minimum, tempTrain, numUsers, numMovies, numFeatures, 0) / size(tempTrain)[2]
		@show score = cost(results.minimum, cross, numUsers, numMovies, numFeatures, 0) / size(cross)[2]
		scoreboard = [scoreboard; λ score time]
	end

	#writecsv("lambba-scoreboard.csv", scoreboard)

	@printf "Done :D\n"
end

main()