text:
function solve_boundary_condition(ω::T, k_eff::Complex{T}, eigvectors::Array{Complex{T}}, source::AbstractSource{Acoustic{T,3}}, material::Material{3,Sphere{T,3}}, ::WithoutSymmetry{3};
        basis_order::Int = 2,
        basis_field_order::Int = 4,
        # source_basis_field_order::Int = basis_field_order,
        source_basis_field_order::Int = Int(round(sqrt(size(eigvectors)[end]))) - 1,
        kws...
    ) where T
    # source_basis_field_order is often chosen so that there is the same number of source coefficients a_n as the number of unknowns α_n
    # Before was: source_basis_field_order = min(basis_field_order,Int(round(sqrt(size(eigvectors)[end])))) - 1

    scale_number_density = one(T) - one(T) / material.numberofparticles

    k = ω / source.medium.c

    species = material.species
    S = length(species)
    rs = outer_radius.(species)

    # dim 1 is the (n,n1) indices, dim 2 is the species, dim 3 are the different eigenvectors
    # eigvectors = reshape(eigvectors,(:,S,size(eigvectors)[end]))

    R = outer_radius(material.shape)

    # Set the maximmum order the basis functions
    L = basis_order
    L1 = basis_field_order
    # L1 = Int(sqrt(size(eigvectors,1) / (L+1)^2) - 1)
    Linc = source_basis_field_order
    L2 = 2L + L1

    # the kernel use to wieight the species and the field's basis order.
    Ns = [
        (R - rs[j]) * kernelN3D(l1,k*(R - rs[j]), k_eff*(R - rs[j])) * scale_number_density * number_density(species[j])
    for l1 = 0:L1, j in eachindex(species)] ./ (k^T(2) - k_eff^T(2))

    l1s = [l1 for l = 0:L for m = -l:l for l1 = 0:L1 for m1 = -l1:l1];

    vecs = [
        eigvectors[i] * Ns[l1s[i[1]]+1,i[2]]
    for i in CartesianIndices(eigvectors)];

    # sum over species
    vecs = sum(vecs, dims=2)
    vecs = reshape(vecs, size(eigvectors)[[1,3]])

    # The extinction_matrix
    function gaunt2(dl,dm,l1,m1,l,m,l2,m2)::Complex{T}
        minl3 = max(abs(dl-l),abs(l2-l1),abs(dm-m))
        maxl3 = min(abs(dl+l),abs(l2+l1))
        return if minl3 <= maxl3
            - sum(l3 ->
                gaunt_coefficient(dl,dm,l,m,l3,dm-m) *
                gaunt_coefficient(l2,m2,l1,m1,l3,dm-m)
            , minl3:maxl3)
        else
            zero(Complex{T})
        end
    end

    # in this form: extinction_matrix[(n,n2),(dn,n1)] ==  gaunt2(dl,dm,l1,m1,l,m,l2,m2)
    extinction_matrix = [
            gaunt2(dl,dm,l1,m1,l,m,l2,m2)
        for dl = 0:L for dm = -dl:dl
        for l1 = 0:L1 for m1 = -l1:l1
    for l = 0:L for m = -l:l
    for l2 = 0:L2 for m2 = -l2:l2]

    len = (L1+1)^2 * (L+1)^2
    extinction_matrix = reshape(extinction_matrix, (:,len))

    source_coefficients = regular_spherical_coefficients(source)(Linc,zeros(3),ω)

    forcing = [
        - sum(
            [gaunt_coefficient(dl,dm,l,m,l2,m2) for dl = 0:Linc for dm = -dl:dl] .*
            source_coefficients
        )
    for l = 0:L for m = -l:l
    for l2 = 0:L2 for m2 = -l2:l2]

    α = (extinction_matrix * vecs) \ forcing

    err = norm(forcing - extinction_matrix * vecs * α) / norm(forcing)

    if err > sqrt(eps(T))
        @warn "Extinction equation (like a boundary condition) was solved with an error: $err for the effective wavenumber: $k_eff"
    end

    return α

end

function solve_boundary_condition(ω::T, k_eff::Complex{T}, eigvectors::Array{Complex{T}}, source::AbstractSource{Acoustic{T,3}}, material::Material{3,Sphere{T,3}}, ::AbstractAzimuthalSymmetry{3};
        basis_order::Int = 2,
        basis_field_order::Int = 4,
        # source_basis_field_order::Int = basis_field_order,
        source_basis_field_order::Int = size(eigvectors)[end] - 1,
        kws...
    ) where T
    # source_basis_field_order is often chosen so that there is the same number of source coefficients a_n as the number of unknowns α_n. Before was: min(basis_field_order, size(eigvectors)[end]) - 1

    k = ω / source.medium.c

    scale_number_density = one(T) - one(T) / material.numberofparticles
    species = material.species
    S = length(species)
    rs = outer_radius.(species)

    # dim 1 is the (n,n1) indices, dim 2 is the species, dim 3 are the different eigenvectors
    # eigvectors = reshape(eigvectors,(:,S,size(eigvectors)[end]))

    R = outer_radius(material.shape)

    L = basis_order
    L1 = basis_field_order
    # L1 = Int((3 * size(eigvectors,1) - 3 - (2 - L) * L * (2 + L)) / (3 + 3*L*(2 + L)))
    Linc = source_basis_field_order

    # Linc = 2
    L2 = 2L + L1

    # the kernel use to wieight the species and the field's basis order.
    Ns = [
        (R - rs[j]) * kernelN3D(l1, k*(R - rs[j]), k_eff*(R - rs[j])) * scale_number_density * number_density(species[j])
    for l1 = 0:L1, j in eachindex(species)] ./ (k^T(2) - k_eff^T(2))

    l1s = [l1 for l = 0:L for m = -l:l for l1 = abs(m):L1];

    vecs = [
        eigvectors[i] * Ns[l1s[i[1]]+1,i[2]]
    for i in CartesianIndices(eigvectors)];

    # sum over species
    vecs = sum(vecs, dims=2)
    vecs = reshape(vecs, size(vecs)[[1,3]])

    # The extinction_matrix
    function gaunt2(dl,dm,l1,l,m,l2)::Complex{T}
        minl3 = max(abs(dl-l),abs(l2-l1),abs(dm-m))
        maxl3 = min(abs(dl+l),abs(l2+l1))
        return if minl3 <= maxl3
            - sum(l3 ->
                gaunt_coefficient(dl,dm,l,m,l3,dm-m) *
                gaunt_coefficient(l2,-m,l1,-dm,l3,dm-m)
            , minl3:maxl3)
        else
            zero(Complex{T})
        end
    end

    # in this form: extinction_matrix[(n,n2),(dn,n1)] ==  gaunt2(dl,dm,l1,m1,l,m,l2,m2)
    extinction_matrix = [
            gaunt2(dl,dm,l1,l,m,l2)
        for dl = 0:L for dm = -dl:dl
        for l1 = abs(dm):L1
    for l = 0:L for m = -l:l
    for l2 = abs(m):L2]

    extinction_matrix = reshape(extinction_matrix, (:,size(vecs,1)))

    # we expect the source direction to be aligned with the z-axis
    direction = SVector(zero(T),zero(T),one(T))

    # source = plane_source(psource.medium, psource.position, direction, psource.amplitude[1])
    # source.coefficients(Linc,zeros(3),ω)

    source_coefficients = regular_spherical_coefficients(source)(Linc,zeros(3),ω)

    forcing = [
        - sum(
            [gaunt_coefficient(dl,dm,l,m,l2,-m) for dl = 0:Linc for dm = -dl:dl] .*
            source_coefficients
        )
    for l = 0:L for m = -l:l
    for l2 = abs(m):L2]

    # If the forcing was exactly a plane-source in the z-axis direction:
    # forcing2 = - [
    #     (m != 0 || m2 != 0) ? 0.0im : 4pi * (T(1)*im)^(l2+l) * sqrt((2l2+1)*(2l+1))
    # for l = 0:L for m = -l:l
    # for l2 = 0:L2 for m2 = -l2:l2]
    #
    # findall( abs.(forcing) .> 1e-4)
    # norm(forcing - forcing2) / norm(forcing2)

    α = (extinction_matrix * vecs) \ forcing

    err = norm(forcing - extinction_matrix * vecs * α) / norm(forcing)

    if err > sqrt(eps(T))
        @warn "Extinction equation (like a boundary condition) was solved with an error: $err for the effective wavenumber: $k_eff"
    end

    return α

end

function solve_boundary_condition(ω::T, k_eff::Complex{T}, eigvectors::Array{Complex{T}}, source::AbstractSource{Acoustic{T,3}}, material::Material{3,Sphere{T,3}}, ::RadialSymmetry{3};
        basis_order::Int = 2,
        kws...
    ) where T

    k = ω / source.medium.c

    scale_number_density = one(T) - one(T) / material.numberofparticles
    species = material.species
    S = length(species)
    rs = outer_radius.(species)

    R = outer_radius(material.shape)

    # the kernel use to weight the species and the field's basis order.
    F = sum(
        T(2 *(i[1] - 1) + 1) * (-one(T))^(i[1]-1) * eigvectors[i] * (R - rs[i[2]]) *
        kernelN3D(i[1] - 1, k*(R - rs[i[2]]), k_eff*(R - rs[i[2]])) * scale_number_density * number_density(species[i[2]])
    for i in CartesianIndices(eigvectors)) / (k^T(2) - k_eff^T(2))

    # We expect there to only one component of a regular spherical wave expansion
    source_coefficients = regular_spherical_coefficients(source)(0,zeros(3),ω)

    forcing = source_coefficients / sqrt(T(4pi))

    α = F \ forcing

    return α
end

max_stars_repo_name: arturgower/EffectiveWaves


text:
include("logformat.jl")
convert(::Type{HTTP.Header}, pa::Pair{String,String}) = Pair(SubString(pa[1]), SubString(pa[2]))
sethd(r::HTTP.Response, pa::Pair) = sethd(r, convert(HTTP.Header, pa))
sethd(r::HTTP.Response, pa::HTTP.Header) = HTTP.setheader(r, pa)

const NEWPORT = 8091

@info "Start server which accepts websocket upgrades including with subprotocol " *
      "'xml' and immediately closes, following protocol."
addsubproto("xml")
serverws =  WebSockets.ServerWS(
    (r::HTTP.Request) -> HTTP.Response(200, "OK"),
    (r::HTTP.Request, ws::WebSocket) -> nothing)
@async WebSockets.serve(serverws, "127.0.0.1", NEWPORT, true)
# tas = @async WebSockets.serve(serverws, "127.0.0.1", NEWPORT, true)
# while !istaskstarted(tas);yield();end

sleep(2)

@info "Open client without subprotocol."
sleep(1)
URL = "ws://127.0.0.1:$NEWPORT"
res = WebSockets.open((_)->nothing, URL);
@test res.status == 101

@info "Open client with approved subprotocol."
sleep(1)
URL = "ws://127.0.0.1:$NEWPORT"
res = WebSockets.open((_)->nothing, URL, subprotocol = "xml");
@test res.status == 101

@info "Open with unknown subprotocol."
sleep(1)
res = WebSockets.open((_) -> nothing, URL, subprotocol = "unapproved");
@test res.status == 400

@info "Try to open a websocket with uknown port. Takes a few seconds."
sleep(1)
global caughterr = WebSockets.WebSocketClosedError("")
try
WebSockets.open((_)->nothing, "ws://127.0.0.1:8099");
catch err
    global caughterr = err
end
@test typeof(caughterr) <: WebSocketClosedError
@test caughterr.message == " while open ws|client: connect: connection refused (ECONNREFUSED)"

@info "Try open with uknown scheme."
sleep(1)
global caughterr = ArgumentError("")
try
WebSockets.open((_)->nothing, "ww://127.0.0.1:8099");
catch err
    global caughterr = err
end
@test typeof(caughterr) <: ArgumentError
@test caughterr.msg == " bad argument url: Scheme not ws or wss. Input scheme: ww"

global caughterr = ArgumentError("")
try
WebSockets.open((_)->nothing, "ws://127.0.0.1:8099/svg/#");
catch err
    global caughterr = err
end
@test typeof(caughterr) <: ArgumentError
@test caughterr.msg == " replace '#' with %23 in url: ws://127.0.0.1:8099/svg/#"

@info "start a client websocket that irritates by closing the TCP stream
 connection without a websocket closing handshake. This
 throws an error in the server coroutine, which can't be
 captured on the server side."
sleep(1)
WebSockets.open("ws://127.0.0.1:$(NEWPORT)") do ws
    close(ws.socket)
end

@info "Check that the server is still running regardless."
sleep(1)
res = WebSockets.open((_)->nothing, URL);
@test res.status == 101

@info "Open with a ws client handler that throws a domain error."
sleep(1)
@test_throws DomainError WebSockets.open((_)->sqrt(-2), URL);

@info "Stop the server in morse code."
sleep(1)
put!(serverws.in, "...-.-")
sleep(1)
@info "Emulate a correct first accept response from server, with BufferStream socket."
sleep(1)
req = HTTP.Request()
req.method = "GET"
key = base64encode(rand(UInt8, 16))
resp = HTTP.Response(101)
resp.request = req
sethd(resp,   "Sec-WebSocket-Version" => "13")
sethd(resp, "Upgrade" => "websocket")
sethd(resp, "Sec-WebSocket-Accept" => WebSockets.generate_websocket_key(key))
sethd(resp,   "Connection" => "Upgrade")
servsock = BufferStream()
s = HTTP.Stream(resp, HTTP.Transaction(HTTP.Connection(servsock)))
write(servsock, resp)
function dummywsh(dws::WebSockets.WebSocket{BufferStream})
    close(dws.socket)
    close(dws)
end
@test WebSockets._openstream(dummywsh, s, key) == WebSockets.CLOSED

@info "Emulate an incorrect first accept response from server."
sleep(1)
sethd(resp, "Sec-WebSocket-Accept" => WebSockets.generate_websocket_key(base64encode(rand(UInt8, 16))))
write(servsock, resp)
@test_throws WebSockets.WebSocketError WebSockets._openstream(dummywsh, s, key)
nothing

ext: jl


text: signaschar(x::Real) = signbit(x) ? '-' : '+'


text:
"""
  semiparametricAUC(model_formula, treatment_group, data)

  This function is used to fit semiparametric AUC regression model specified by
  giving a formula object of response and covariates and a separate argument of treatment
  group. It will convert variables other than response into factors, estimate model parameters,
  and display results.

  This function takes a formula object `model_formula` with response and covariates such as {response ~ x1 + x2},
  group argument `treatment_group`, which is treatment group for which a comparision is to be made, and
  a data argument `data` which is a DataFrame object that contains variables needed for the analysis.
"""
semiparametricAUC(model_formula = y ~ x1, treatment_group = :group, data = fasd)

function semiparametricAUC(; model_formula::DataFrames.Formula = throw(ArgumentError("Argument model_formula is missing")),
  treatment_group::Symbol = throw(ArgumentError("Argument treatment_group is missing")),
    data::DataFrames.DataFrame = throw(ArgumentError("Argument data is missing")))
  # fasd = DataFrames.readtable(joinpath(Pkg.dir("SemiparametricAUC"), "data/fasd.csv"))
  # model_formula = y ~ x1 + x2
  # treatment_group = :group
  # data = fasd

  if (isa(model_formula, Formula))
    input_covariates = DataFrames.Terms(model_formula).terms
    n1 = length(input_covariates)
    input_response = Symbol(DataFrames.Terms(model_formula).eterms[1])
  else error("Please put response and input as DataFrames.Formula object. For example, model_formula = response ~ x1 + x2")
  end

  if (!isa(treatment_group, Symbol))
    error("The parameter treatment_group should be Symbol object. For e.g. :x1")
  end

  if (!isa(data, DataFrames.DataFrame))
    error("The parameter data should be DataFrames.DataFrame object.")
  end
  input_treatment = Symbol(treatment_group)
  group_covariates = vcat(input_covariates, input_treatment)

  # if (sum([isa(data[:,i], DataFrames.PooledDataArray) for i in group_covariates]) == 0)
  #   error("Please put response and input as formula. For example, response ~ x1 + x2")
  # else
  #     println("Great")
  # end

  print_with_color(:red,"Data are being analyzed. Please, be patient.\n\n")
  # split by factors
  # TO-DO: make sure that oerder of the variables are aligned with order of coefnames(mf)
  ds = data[vcat(input_response,Symbol.(group_covariates))]
  grouped_d = DataFrames.groupby(ds, group_covariates)

  half_data = 0.5
  set1 =  Dict()
  for i in 1:2:length(grouped_d)
    set1[i] = grouped_d[i]
  end
  set1_sorted = DataStructures.SortedDict(set1)

  set2 = Dict()
  for i in 2:2:length(grouped_d)
    set2[i] = grouped_d[i]
  end
  set2_sorted = DataStructures.SortedDict(set2)

  # TO-DO: make sure that order of the variables are aligned with order of coefnames(mf)
  logitauchat_matrix = collect(calculate_auc(ya = set1_sorted[i][:,input_response], yb = set2_sorted[i+1][:,input_response])
    for i in 1:2:length(grouped_d))

  dff = DataFrames.DataFrame([y[i] for y in logitauchat_matrix, i in 1:length(logitauchat_matrix[1])])
  var_logitauchat = dff[1]
  gamma1 = dff[2]

  # Change to adjust more than 2 d array
  # get levels
  dict_levels = Dict()
  for i in 1:length(input_covariates)
    dict_levels[i] = DataArrays.levels(ds[input_covariates[i]])
  end

  # for expand.grid
  matrix_x = collect(Iterators.product(values(sort(dict_levels))...))

  df_from_tuple = DataFrames.DataFrame([y[i] for y in matrix_x, i in 1:length(matrix_x[1])])
  df_from_tuple[input_response] = gamma1
  df_from_tuple[:var_logit] = var_logitauchat

  function convert_to_factor(x)
      return(DataFrames.pool(x))
  end

  for i in input_covariates
    df_from_tuple[i] = convert_to_factor(df_from_tuple[i])
  end

  # model.matrix using DataFrames (ModelMatrix)
  mf = ModelFrame(DataFrames.Terms(model_formula), df_from_tuple)
  mm = ModelMatrix(mf)

  Z = mm.m
  tau  =  diagm([1/i for i in var_logitauchat])
  ztauz = inv(Z' * tau * Z)
  var_betas = diag(ztauz)
  std_error = sqrt(var_betas)
  betas = ztauz * Z' * tau * gamma1

  threshold = Distributions.quantile(Distributions.Normal(), 0.975)

  lo = betas - threshold*std_error
  up = betas + threshold*std_error
  ci = hcat(betas,lo,up)

  return(SemiparametricAUC.coefs_table(mf = mf, lo = lo, up = up, betass = betas, std_errors = std_error))
end


text:
input_string = join(ARGS, " ")
println(input_string[end:-1:1])


text:
const DOC_PLOT_PARAMS = """
- **`title`** : Text to display on the top of the plot.

- **`xlabel`** : Text to display on the x axis of the plot

- **`ylabel`** : Text to display on the y axis of the plot

- **`labels`** : Boolean. Can be used to hide the labels by
  setting `labels = false`.

- **`border`** : The style of the bounding box of the plot.
  Supports `:corners`, `:solid`, `:bold`, `:dashed`, `:dotted`,
  `:ascii`, and `:none`.

- **`margin`** : Number of empty characters to the left of the
  whole plot.

- **`padding`** : Space of the left and right of the plot between
  the labels and the canvas.

- **`color`** : Color of the drawing.
  Can be any of `:green`, `:blue`, `:red`, `:yellow`, `:cyan`,
  `:magenta`, `:white`, `:normal`

- **`width`** : Number of characters per row that should be used
  for plotting.
"""

transform_name(::typeof(identity), basename = "") = basename
function transform_name(f, basename = "")
    fname = string(nameof(f))
    fname = occursin("#", fname) ? "custom" : fname
    string(basename, " [", fname, "]")
end

roundable(num::Number) = isinteger(num) & (typemin(Int) <= num < typemax(Int))

ceil_neg_log10(x) = roundable(-log10(x)) ? ceil(Integer, -log10(x)) : floor(Integer, -log10(x))
round_neg_log10(x) = roundable(-log10(x)) ? round(Integer, -log10(x), RoundNearestTiesUp) : floor(Integer, -log10(x))
round_up_tick(x,m) = x == 0. ? 0. : (x > 0 ? ceil(x, digits=ceil_neg_log10(m)) : -floor(-x, digits=ceil_neg_log10(m)))
round_down_tick(x,m) = x == 0. ? 0. : (x > 0 ? floor(x, digits=ceil_neg_log10(m)) : -ceil(-x, digits=ceil_neg_log10(m)))
round_up_subtick(x,m) = x == 0. ? 0. : (x > 0 ? ceil(x, digits=ceil_neg_log10(m)+1) : -floor(-x, digits=ceil_neg_log10(m)+1))
round_down_subtick(x,m) = x == 0. ? 0. : (x > 0 ? floor(x, digits=ceil_neg_log10(m)+1) : -ceil(-x, digits=ceil_neg_log10(m)+1))
float_round_log10(x::F,m) where {F<:AbstractFloat} = x == 0. ? F(0) : (x > 0 ? round(x, digits=ceil_neg_log10(m)+1)::F : -round(-x, digits=ceil_neg_log10(m)+1)::F)
float_round_log10(x) = x > 0 ? float_round_log10(x,x) : float_round_log10(x,-x)

function plotting_range(xmin, xmax)
    diff = xmax - xmin
    xmax = round_up_tick(xmax, diff)
    xmin = round_down_tick(xmin, diff)
    Float64(xmin), Float64(xmax)
end

function plotting_range_narrow(xmin, xmax)
    diff = xmax - xmin
    xmax = round_up_subtick(xmax, diff)
    xmin = round_down_subtick(xmin, diff)
    Float64(xmin), Float64(xmax)
end

function extend_limits(vec, limits)
    mi, ma = map(Float64, extrema(limits))
    if mi == 0. && ma == 0.
        mi, ma = map(Float64, extrema(vec))
    end
    diff = ma - mi
    if diff == 0
        ma = mi + 1
        mi = mi - 1
    end
    (limits == (0.,0.) || limits == [0.,0.]) ? plotting_range_narrow(mi, ma) : (mi, ma)
end

const bordermap = Dict{Symbol,Dict{Symbol,String}}()
const border_solid  = Dict{Symbol,String}()
const border_corners = Dict{Symbol,String}()
const border_barplot = Dict{Symbol,String}()
const border_bold   = Dict{Symbol,String}()
const border_none   = Dict{Symbol,String}()
const border_dashed = Dict{Symbol,String}()
const border_dotted = Dict{Symbol,String}()
const border_ascii  = Dict{Symbol,String}()
border_solid[:tl] = "┌"
border_solid[:tr] = "┐"
border_solid[:bl] = "└"
border_solid[:br] = "┘"
border_solid[:t] = "─"
border_solid[:l] = "│"
border_solid[:b] = "─"
border_solid[:r] = "│"
border_corners[:tl] = "┌"
border_corners[:tr] = "┐"
border_corners[:bl] = "└"
border_corners[:br] = "┘"
border_corners[:t] = " "
border_corners[:l] = " "
border_corners[:b] = " "
border_corners[:r] = " "
border_barplot[:tl] = "┌"
border_barplot[:tr] = "┐"
border_barplot[:bl] = "└"
border_barplot[:br] = "┘"
border_barplot[:t] = " "
border_barplot[:l] = "┤"
border_barplot[:b] = " "
border_barplot[:r] = " "
border_bold[:tl] = "┏"
border_bold[:tr] = "┓"
border_bold[:bl] = "┗"
border_bold[:br] = "┛"
border_bold[:t] = "━"
border_bold[:l] = "┃"
border_bold[:b] = "━"
border_bold[:r] = "┃"
border_none[:tl] = " "
border_none[:tr] = " "
border_none[:bl] = " "
border_none[:br] = " "
border_none[:t] = " "
border_none[:l] = " "
border_none[:b] = " "
border_none[:r] = " "
border_dashed[:tl] = "┌"
border_dashed[:tr] = "┐"
border_dashed[:bl] = "└"
border_dashed[:br] = "┘"
border_dashed[:t] = "╌"
border_dashed[:l] = "┊"
border_dashed[:b] = "╌"
border_dashed[:r] = "┊"
border_dotted[:tl] = "⡤"
border_dotted[:tr] = "⢤"
border_dotted[:bl] = "⠓"
border_dotted[:br] = "⠚"
border_dotted[:t] = "⠤"
border_dotted[:l] = "⡇"
border_dotted[:b] = "⠒"
border_dotted[:r] = "⢸"
border_ascii[:tl] = "+"
border_ascii[:tr] = "+"
border_ascii[:bl] = "+"
border_ascii[:br] = "+"
border_ascii[:t] = "-"
border_ascii[:l] = "|"
border_ascii[:b] = "-"
border_ascii[:r] = "|"
bordermap[:solid]  = border_solid
bordermap[:corners] = border_corners
bordermap[:barplot] = border_barplot
bordermap[:bold]   = border_bold
bordermap[:none]   = border_none
bordermap[:dashed] = border_dashed
bordermap[:dotted] = border_dotted
bordermap[:ascii]  = border_ascii

const color_cycle = [:green, :blue, :red, :magenta, :yellow, :cyan]
const color_encode = Dict{Symbol,UInt8}()
const color_decode = Dict{UInt8,Symbol}()
color_encode[:normal]  = 0b000
color_encode[:blue]    = 0b001
color_encode[:red]     = 0b010
color_encode[:magenta] = 0b011
color_encode[:green]   = 0b100
color_encode[:cyan]    = 0b101
color_encode[:yellow]  = 0b110
for k in keys(color_encode)
    v = color_encode[k]
    color_decode[v] = k
end
color_encode[:white] = 0b111
color_decode[0b111]  = :white

function print_color(color::UInt8, io::IO, args...)
    col = color in keys(color_decode) ? color_decode[color] : Int(color)
    str = string(args...)
    printstyled(io, str; color = col)
end


text:
using Pkg
function is_dqmc_env_activated()
    project_file = Base.active_project()
    project = Pkg.Types.read_project(project_file)
    return !isnothing(project.name) && lowercase(project.name) == "qmc"
end

if !is_dqmc_env_activated()
  println("Activating DQMC environment.")
  haskey(ENV, "JULIA_DQMC") || error("DQMC environment not loaded and JULIA_DQMC env variable not set!")
  Pkg.activate(ENV["JULIA_DQMC"])
  is_dqmc_env_activated() || error("Invalid JULIA_DQMC env variable value.")
end

# #######################################################
#                       Framework
# #######################################################
# -------------------------------------------------------
#                       Includes
# -------------------------------------------------------
global const TIMING = true
include("dqmc_framework.jl")
using Parameters
using ProgressMeter
using TimerOutputs
using BinningAnalysis
using RecursiveArrayTools

const OBSERVABLES = Set((:chi_dyn,
                         :chi_dyn_symm,
                         :chi_dyn_components,
                         :binder,
                         :etpc,
                         :zfpc,
                         :etcdc,
                         :zfcdc,
                         :zfccc_xx,
                         :zfccc_yy,
                         :sfdensity,
                         :greens,
                         :tdgfs_gt0,
                         :tdgfs_g0t
))

const PLUSMINUS_OBSERVABLES = Set((:etpc,
                                   :zfpc,
                                   :etcdc,
                                   :zfcdc
))

const CALC_IN_ONE_GO = Set((:chi_dyn,       # always in one go
                            :chi_dyn_symm,  # always in one go
                            :chi_dyn_components, # always in one go
                            :binder         # always in one go
))

# -------------------------------------------------------
#                  Define Macros
# -------------------------------------------------------
macro addlightobs(name, zero)
  nc = esc(:(num_confs))
  o_inonego = esc(:(obs_inonego))
  o_insteps = esc(:(obs_insteps))
  mp = esc(:(mp))
  n = "$name"
  z = esc(:($zero))
  return quote
    if Symbol($n) in $(mp).todo_insteps
      $(o_insteps) = add($(o_insteps), $name = create_or_load_lightobs($n, $z, name = $n, alloc = $nc))
    else
      $(o_inonego) = add($(o_inonego), $name = LightObservable($z, name = $n, alloc = $nc))
    end
  end
end

macro addobs(name, T)
  nc = esc(:(num_confs))
  o_inonego = esc(:(obs_inonego))
  o_insteps = esc(:(obs_insteps))
  mp = esc(:(mp))
  n = "$name"
  t = esc(:($T))
  return quote
    if Symbol($n) in $(mp).todo_insteps
      $(o_insteps) = add($(o_insteps), $name = create_or_load_obs($n, $t, name = $n, alloc = $nc))
    else
      $(o_inonego) = add($(o_inonego), $name = Observable($t, name = $n, alloc = $nc))
    end
  end
end

# -------------------------------------------------------
#                  Define MeasParams
# -------------------------------------------------------
@with_kw mutable struct MeasParams
  requested::Set{Symbol} = Set{Symbol}()
  todo::Set{Symbol} = Set{Symbol}()
  todo_inonego::Set{Symbol} = Set{Symbol}()
  todo_insteps::Set{Symbol} = Set{Symbol}()

  # input & output files
  meas_infile::String = ""
  outfile::String = ""
  objfile::String = "" # for temporarily storing full observable objetcs
  p::Params = Params() # DQMC Params
  dqmc_outfile::String = ""

  # parameters
  overwrite::Bool = false
  num_threads::Int = 1
  walltimelimit::Dates.DateTime = Dates.DateTime("2099", "YYYY")
  save_after::Int = 10
  to::TimerOutput = TimerOutput()
  ignore_first_n::Int = 0
  confs_iterator_start::Int = 1 # this is the S in "for i in S:length(confs)"
  safe_mult::Union{Nothing, Int} = nothing # nothing = use what has been used for dqmc
end

function measxml2MeasParams(fname)
  # meas.xml: direct mapping of xml fields to kwargs
  mpdict = Dict(Symbol(lowercase(k))=>lowercase(v) for (k,v) in xml2dict(fname, false))
  mp = MeasParams()

  for (k, v) in mpdict
    if k in OBSERVABLES
      if parse(Bool, v)
        if k in PLUSMINUS_OBSERVABLES
          push!(mp.requested, Symbol(string(k)*"_minus"))
          push!(mp.requested, Symbol(string(k)*"_plus"))
        else
          push!(mp.requested, k)
        end
      end

    elseif k == :overwrite
      mp.overwrite = parse(Bool, v)
    elseif k == :num_threads
      mp.num_threads = parse(Int64, v)
    elseif k == :save_after
      mp.save_after = parse(Int64, v)
    elseif k == :ignore_first_n
      mp.ignore_first_n = parse(Int64, v)
      mp.confs_iterator_start = mp.ignore_first_n + 1
    elseif k == :safe_mult
      mp.safe_mult = parse(Int64, v)
    end
  end

  mp.meas_infile = fname

  mp.dqmc_outfile = replace(mp.meas_infile, ".meas.xml" => ".out.h5")
  if !isfile(mp.dqmc_outfile)
      mp.dqmc_outfile = replace(mp.dqmc_outfile, ".out.h5" => ".out.h5.running")
      isfile(mp.dqmc_outfile) || error("Couldn't find .out.h5[.running] file.")
  end
  mp.outfile = replace(mp.dqmc_outfile, ".out.h5" => ".meas.h5")
  mp.objfile = replace(mp.dqmc_outfile, ".out.h5" => ".meas.obj.h5")

  dqmc_infile = fname[1:end-9]*".in.xml"
  if isfile(dqmc_infile)
    println("Loading DQMC parameters (.in.xml)")
    xml2parameters!(mp.p, dqmc_infile)
  else
    println("Loading DQMC parameters (.out.h5)")
    try
      hdf52parameters!(mp.p, mp.dqmc_outfile)
    catch er
      println("Error while loading parameters from .out.h5 file")
      throw(er)
    end
  end

  if @isdefined start_time
    set_walltimelimit!(mp, start_time)
  else
    set_walltimelimit!(mp, now())
  end

  return mp
end

@inline function need_to_setup_mc(mp::MeasParams)
  symb = (:etpc_plus, :zfpc_plus, :etcdc_plus, :zfcdc_plus, :zfccc_xx, :zfccc_yy, :sfdensity, :greens, :tdgfs_gt0, :tdgfs_g0t)
  return any(in.(symb, Ref(mp.todo)))
end
@inline function need_etgf(mp::MeasParams)
  symb = (:etpc_plus, :etcdc_plus, :zfcdc_plus, :zfccc_xx, :zfccc_yy, :sfdensity, :greens)
  return any(in.(symb, Ref(mp.todo)))
end
@inline function need_to_meas_tdgfs(mp::MeasParams)
  symb = (:zfpc_plus, :zfcdc_plus, :zfccc_xx, :zfccc_yy, :sfdensity, :tdgfs_gt0, :tdgfs_g0t)
  return any(in.(symb, Ref(mp.todo)))
end
@inline function need_to_meas_zfccc(mp::MeasParams)
  symb = (:zfccc_xx, :zfccc_yy, :sfdensity)
  return any(in.(symb, Ref(mp.todo)))
end

@inline function need_etgf(ol::NamedTuple{K,V}) where {K,V}
  symb = (:etpc_plus, :etcdc_plus, :zfcdc_plus, :zfccc_xx, :zfccc_yy, :sfdensity, :greens)
  return any(Base.sym_in.(symb, Ref(K)))
end
@inline function need_to_meas_tdgfs(ol::NamedTuple{K,V}) where {K,V}
  symb = (:zfpc_plus, :zfcdc_plus, :zfccc_xx, :zfccc_yy, :sfdensity, :tdgfs_gt0, :tdgfs_g0t)
  return any(Base.sym_in.(symb, Ref(K)))
end
@inline function need_to_meas_zfccc(ol::NamedTuple{K,V}) where {K,V}
  return Base.sym_in(:sfdensity, K) || Base.sym_in(:zfccc_xx, K) || Base.sym_in(:zfccc_yy, K)
end

# -------------------------------------------------------
#              Initialize stack for meas
# -------------------------------------------------------
function initialize_stack_for_measurements(mc::AbstractDQMC, mp::MeasParams)
    _initialize_stack(mc)

    # allocate for measurement run based on todos
    :etpc_plus in mp.todo  && allocate_etpc!(mc)
    :zfpc_plus in mp.todo  && allocate_zfpc!(mc)
    :etcdc_plus in mp.todo && allocate_etcdc!(mc)
    :zfcdc_plus in mp.todo && allocate_zfcdc!(mc)

    need_to_meas_zfccc(mp) && allocate_zfccc!(mc)
    need_to_meas_tdgfs(mp) && allocate_tdgfs!(mc)

    nothing
end

"""
Groups todos according to whether they should be calculated
 * in one go
 * step by step with intermediate saves
"""
function create_todolist!(mp::MeasParams)
  mp.todo = Set{Symbol}()

  if !endswith(mp.outfile, ".running") && isfile(mp.outfile*".running")
    # job has finished since last measurement run
    mv(mp.outfile*".running", mp.outfile)
  end

  if isfile(mp.outfile) && !mp.overwrite
    # Start/continue unfinished measurements but don't overwrite anything
    println("Measurement file found.");
    println("Checking what is left todo...")

    local confs
    try
      confs = loadobs_frommemory(mp.dqmc_outfile, "obs/configurations")
    catch err
      println("Couldn't read configuration data. Probably no configurations yet? Exiting.")
      exit()
    end
    allobs = listobs(mp.outfile)
    N = length(confs)
    counts = Dict{Symbol, Int64}() # # measurement results in .meas.h5
    objcounts = Dict{Symbol, Int64}() # length of binner obj in .meas.obj.h5

    jldopen(mp.outfile, "r") do f
      for o in mp.requested
        o == :binder && continue # special case binder cumulant
        p = joinpath("obs/", string(o))

        if !HDF5.has(f.plain, p)
          push!(mp.todo, o)
        else
          count = read(f[joinpath(p, "count")])
          if count != N
            if strip_plusminus(o) in CALC_IN_ONE_GO
              push!(mp.todo, o)
            else # CALC_IN_STEPS
              push!(mp.todo, o)
              counts[o] = count
            end
          end
        end
      end # observable loop
    end # open file

    if isfile(mp.objfile)
      jldopen(mp.objfile, "r") do fobj
        for o in mp.requested
          o == :binder && continue # special case binder cumulant
          if !(strip_plusminus(o) in CALC_IN_ONE_GO)
            # CALC_IN_STEPS
            p = string(o)
            if HDF5.has(fobj.plain, p)
              objcounts[o] = length(read(fobj[p]))
            end
          end
        end
      end # open file
    end

    # Special case: binder
    if :binder in mp.requested && :chi_dyn in mp.todo
      push!(mp.todo, :binder)
    end

    if length(mp.todo) > 0
      c = isempty(counts) ? 0 : first(values(counts))
      objc = isempty(objcounts) ? 0 : first(values(objcounts))

      if length(unique(values(objcounts))) > 1 || length(unique(values(counts))) > 1
        @show counts
        @show objcounts
        @error "Found unfinished measurements of different length (not yet supported). Aborting."
      end

      println("Found ", c, " in-steps-measurement results in outfile.")
      println("Found ", objc, " in-steps-measurements in binners in objfile.")

      if c>0 && c!=objc
        @show counts
        @show objcounts
        @error "Measurement results present + incompatible binner objects (unsupported)."
      end

      mp.confs_iterator_start = mp.ignore_first_n + objc + 1
      @show mp.confs_iterator_start
    else
      mp.confs_iterator_start = -1
    end

  else # there is no meas outfile yet or we are in overwrite mode
    mp.todo = copy(mp.requested)
    isfile(mp.objfile) && rm(mp.objfile)
  end

  mp.todo
end

function group_todos!(mp)
  mp.todo_inonego

  for o in mp.todo
    if strip_plusminus(o) in CALC_IN_ONE_GO
      push!(mp.todo_inonego, o)
    else
      push!(mp.todo_insteps, o)
    end
  end

  nothing
end

strip_plusminus(str::AbstractString) = replace(replace(str, "_plus" => ""), "_minus" => "")
strip_plusminus(s::Symbol) = Symbol(strip_plusminus(string(s)))

# -------------------------------------------------------
#                        Main
# -------------------------------------------------------
function measure(mp::MeasParams; debug=false)

  create_todolist!(mp)

  if isempty(mp.todo)
    println("Nothing on the todo list. Exiting.")
    exit()
  end

  group_todos!(mp)

  @show mp.todo_inonego
  @show mp.todo_insteps
  flush(stdout)

  # Overwrite mode?
  if isfile(mp.outfile) && mp.overwrite
    if mp.overwrite && !isempty(mp.todo_insteps)
      println("\nOVERWRITE MODE: Ignored to avoid bug #23.")
    else
      println("\nOVERWRITE MODE")
      isfile(mp.objfile) && rm(mp.objfile)
      jldopen(mp.outfile, "r+") do f
        for o in mp.todo
          p = joinpath("obs/", string(o))
          if HDF5.has(f.plain, p)
            println("Clearing $p")
            o_delete(f, p)
          end
        end
      end
      println()
    end
  end

  # -------------- Load DQMC params/results -----------------
  println("\nLoading DQMC results...")

  greens = nothing;
  mc = nothing;
  local confs

  confs = loadobs_frommemory(mp.dqmc_outfile, "obs/configurations")

  if need_etgf(mp)
    if "greens" in listobs(mp.dqmc_outfile)
      greens = loadobs_frommemory(mp.dqmc_outfile, "obs/greens")
    else
      # error("Couldn't find ETGF in $(mp.dqmc_outfile) although need_etgf == true.")
      println("Couldn't find ETGF in $(mp.dqmc_outfile). Will calculate it on-the-fly.")
      # TODO: Calculate ETGF if it isn't found in dqmc output.
    end
  end

  # prepare dqmc framework if necessary
  if need_to_setup_mc(mp)
    if !isnothing(mp.safe_mult)
      mp.p.safe_mult = mp.safe_mult
      @show mp.p.safe_mult
    end
    mc = DQMC(mp.p)
    initialize_stack_for_measurements(mc, mp)
  end

  # ------------------- Create list of observables to measure --------------------
  println("\nPreparing observables......................")
  num_confs = length(confs)
  nsweeps = num_confs * mp.p.write_every_nth
  obs_inonego = NamedTuple()
  obs_insteps = NamedTuple()

  # chi_dyn
  if :chi_dyn_symm in mp.todo
    @addobs chi_dyn_symm Array{Float64, 3}
  end

  if :chi_dyn in mp.todo
    @addobs chi_dyn Array{Float64, 3}
  end

  if :chi_dyn_components in mp.todo
    @addobs chi_dyn_components Array{Float64, 4}
  end

  # binder
  if :binder in mp.todo
    @addobs m2s Float64
    @addobs m4s Float64
    obs_inonego = add(obs_inonego, binder = Observable(Float64, name="binder"))
  end

  # etpc
  if :etpc_plus in mp.todo
    zero_etpc = zeros(Float64, mp.p.L, mp.p.L)
    @addlightobs etpc_plus zero_etpc
    @addlightobs etpc_minus zero_etpc
  end

  # zfpc
  if :zfpc_plus in mp.todo
    zero_zfpc = zeros(Float64, mp.p.L, mp.p.L)
    @addlightobs zfpc_plus zero_zfpc
    @addlightobs zfpc_minus zero_zfpc
  end

  # etcdc
  if :etcdc_plus in mp.todo
    zero_etcdc = zeros(ComplexF64, mp.p.L, mp.p.L)
    @addlightobs etcdc_plus zero_etcdc
    @addlightobs etcdc_minus zero_etcdc
  end

  # zfcdc
  if :zfcdc_plus in mp.todo
    zero_zfcdc = zeros(ComplexF64, mp.p.L, mp.p.L)
    @addlightobs zfcdc_plus zero_zfcdc
    @addlightobs zfcdc_minus zero_zfcdc
  end

  # zfccc
  if :zfccc in mp.todo
    zero_zfccc = zeros(ComplexF64, mp.p.L, mp.p.L)
    @addlightobs zfccc_xx zero_zfccc
    @addlightobs zfccc_yy zero_zfccc
  end

  # sfdensity
  if :sfdensity in mp.todo
    @addlightobs sfdensity Float64
  end

  # greens
  if :greens in mp.todo
    @addobs greens Matrix{geltype(mc)}
  end

  # tdgfs
  if :tdgfs_g0t in mp.todo
    Nflv = mc.p.flv * mc.l.sites
    tdgf_size = (Nflv, Nflv, mc.p.slices)
    zero_tdgf = zeros(geltype(mc), tdgf_size)
    @addlightobs tdgfs_gt0 zero_tdgf
    @addlightobs tdgfs_g0t zero_tdgf
  end

  print("........................... Done. ");
  println("\n")
  @show keys(obs_inonego)
  @show keys(obs_insteps)
  println("\n")

  # --------------------- Measure -------------------------
  if length(obs_inonego) > 0 && !debug
    reset_timer!(mp.to)
    measure_inonego(mp, obs_inonego, confs, greens, mc)
    export_results(mp, obs_inonego, nsweeps)
  end

  if length(obs_insteps) > 0 && !debug
    reset_timer!(mp.to)
    measure_insteps(mp, obs_insteps, confs, greens, mc)
    export_results(mp, obs_insteps, nsweeps)
  end

  debug && (return obs_inonego, obs_insteps)

  nothing
end

# -------------------------------------------------------
#                 MEASUREMENT LOOPS
# -------------------------------------------------------
function measure_inonego(mp::MeasParams, obs::NamedTuple{K,V}, confs, greens, mc) where {K,V}
    length(obs) > 0 && println("---------------------- Measuring (in-one-go) ----------------------");
    flush(stdout)

    @mytimeit mp.to "measure loop" begin
      @inbounds @views @showprogress for i in 1:length(confs)
          # println(i); flush(stdout);
          @mytimeit mp.to "load conf" conf = confs[i]

          @mytimeit mp.to "bosonic" measure_bosonic(mp, obs, conf, i)

          @mytimeit mp.to "fermionic" measure_fermionic(mp, obs, conf, greens, mc, i)
      end
    end

    if TIMING
      println(); display(mp.to); println(); flush(stdout);
    end

    if :binder in keys(obs)
        # binder postprocessing
        m2ev2 = mean(obs[:m2s])^2
        m4ev = mean(obs[:m4s])

        push!(obs[:binder], m4ev/m2ev2)
    end
end

function measure_insteps(mp::MeasParams, obs::NamedTuple{K,V}, confs, greens, mc) where {K,V}
    length(obs) > 0 && println("---------------------- Measuring (in-steps) ----------------------");
    flush(stdout)

    @inbounds @views @showprogress for i in mp.confs_iterator_start:length(confs)
        @mytimeit mp.to "iteration" begin
          println(i); flush(stdout);
          # sleep(2)
          @mytimeit mp.to "load conf" conf = confs[i]

          @mytimeit mp.to "bosonic" measure_bosonic(mp, obs, conf, i)

          @mytimeit mp.to "fermionic" measure_fermionic(mp, obs, conf, greens, mc, i)

          @mytimeit mp.to "intermediate save" if mod1(i, mp.save_after) == mp.save_after
            save_obs_objects(mp, obs)

            # Estimate whether we'll make it to another save before hitting WTL. If not, flush and restart.
            titer = mp.to["iteration"]
            iterdur = TimerOutputs.time(titer) *10.0^(-9)/TimerOutputs.ncalls(titer)
            secs_to_save = mp.save_after * iterdur
            secs_to_save *= 1.1 # add 10 percent buffer
            next_save_date = now() + Millisecond(ceil(Int, secs_to_save * 1000))

            if next_save_date >= mp.walltimelimit
              println("Approaching wall-time limit. Won't make it to next save.")
              println("Exporting intermediate results.")
              export_results(mp, obs, nothing)
              println("Safely exiting. (i = $(i)). Current date: $(Dates.format(now(), "d.u yyyy HH:MM")).")
              exit(42)
            end
          end # if save
        end # time one iteration
    end

    # @mytimeit mp.to "final save" save_obs_objects(mp, obs)
    # if isfile(mp.objfile)
    #   println("Deleting obj file: ", mp.objfile)
    #   rm(mp.objfile)
    # end

    if TIMING
      println(); display(mp.to); println(); flush(stdout);
    end
end

function save_obs_objects(mp, obs)
  println("Intermediate save..."); flush(stdout)
  jldopen(mp.objfile*"_tmp", "w") do f # always overwrite => "w"
    for (o, obj) in pairs(obs)
      f[string(o)] = obj
    end
  end
  mv(mp.objfile*"_tmp", mp.objfile; force=true)
  # OPT: Be smarter for Observable objects (bosonic observables)
  nothing
end

# -------------------------------------------------------
#                      MEASUREMENTS
# -------------------------------------------------------
function measure_bosonic(mp, obs, conf, i)
      # chi_dyn
      @mytimeit mp.to "chi" if :chi_dyn in keys(obs)
        chi = measure_chi_dynamic(conf)
        push!(obs[:chi_dyn], chi)

        if :chi_dyn_symm in keys(obs)
          chi = (permutedims(chi, [2,1,3]) + chi)/2 # C4 is basically flipping qx and qy (which only go from 0 to pi since we perform a real fft.)
          push!(obs[:chi_dyn_symm], chi)
        end

        if :chi_dyn_components in keys(obs)
          chi_components = measure_chi_dynamic_components(conf)
          push!(obs[:chi_dyn_components], chi_components)
        end
      end

      # binder
      @mytimeit mp.to "binder" if :binder in keys(obs)
        m = mean(conf, dims=(2,3))
        push!(obs[:m2s], dot(m, m))
        push!(obs[:m4s], last(obs[:m2s])^2)
      end
      nothing
end

function measure_fermionic(mp, obs::NamedTuple{K,V}, conf, greens, mc, i) where {K,V}
    if !isnothing(mc)
      mc.p.hsfield = conf
    end

    if need_etgf(obs)
      if !isnothing(greens)
        @mytimeit mp.to "load etgf" etgf = greens[i] # MCO.jl: load from disk
      else
        @mytimeit mp.to "calculate Greens" etgf = measure_greens(mc)
      end
    end

    # greens
    @mytimeit mp.to "greens" if :greens in K
        push!(obs[:greens], etgf)
    end

    # etpc
    @mytimeit mp.to "etpc" if :etpc_plus in K
        measure_etpc!(mc, etgf)
        push!(obs[:etpc_plus], mc.s.meas.etpc_plus)
        push!(obs[:etpc_minus], mc.s.meas.etpc_minus)
    end

    # etcdc
    @mytimeit mp.to "etcdc" if :etcdc_plus in K
        measure_etcdc!(mc, etgf)
        push!(obs[:etcdc_plus], mc.s.meas.etcdc_plus)
        push!(obs[:etcdc_minus], mc.s.meas.etcdc_minus)
    end

    @mytimeit mp.to "calculate TDGFs" if need_to_meas_tdgfs(obs)
      measure_tdgfs!(mc)
      GC.gc()
      Gt0 = mc.s.meas.Gt0
      G0t = mc.s.meas.G0t
    end

    @mytimeit mp.to "zfccc" if need_to_meas_zfccc(obs)
      measure_zfccc!(mc, etgf, Gt0, G0t)
      zfccc_xx = mc.s.meas.zfccc_xx
      zfccc_yy = mc.s.meas.zfccc_yy
    end

    # zfpc
    @mytimeit mp.to "zfpc" if :zfpc_plus in K
        measure_zfpc!(mc, Gt0)
        push!(obs[:zfpc_plus], mc.s.meas.zfpc_plus)
        push!(obs[:zfpc_minus], mc.s.meas.zfpc_minus)
    end

    # zfcdc
    @mytimeit mp.to "zfcdc" if :zfcdc_plus in K
        measure_zfcdc!(mc, etgf, Gt0, G0t)
        push!(obs[:zfcdc_plus], mc.s.meas.zfcdc_plus)
        push!(obs[:zfcdc_minus], mc.s.meas.zfcdc_minus)
    end

    # sfdensity
    @mytimeit mp.to "sfdensity" if :sfdensity in K
        push!(obs[:sfdensity], measure_sfdensity(mc, zfccc_xx, zfccc_yy))
    end

    # tdgfs
    if :tdgfs_gt0 in keys(obs)
        push!(obs[:tdgfs_gt0], VectorOfArray(Gt0))
        push!(obs[:tdgfs_g0t], VectorOfArray(G0t))
    end

    # zfccc
    if :zfccc_xx in K || :zfccc_yy in K
        push!(obs[:zfccc_xx], zfccc_xx)
        push!(obs[:zfccc_yy], zfccc_yy)
    end

    nothing
end

# -------------------------------------------------------
#                      Export
# -------------------------------------------------------
function export_results(mp, obs, nsweeps)
    println("Calculating errors and exporting..."); flush(stdout)

    # bosonic
    :chi_dyn in keys(obs) && export_result(obs[:chi_dyn], mp.outfile, "obs/chi_dyn"; timeseries=true)
    :chi_dyn_symm in keys(obs) && export_result(obs[:chi_dyn_symm], mp.outfile, "obs/chi_dyn_symm"; timeseries=true)
    :chi_dyn_components in keys(obs) && export_result(obs[:chi_dyn_components], mp.outfile, "obs/chi_dyn_components"; timeseries=true)
    :binder in keys(obs) && export_result(obs[:binder], mp.outfile, "obs/binder", error=false) # jackknife for error

    # fermionic
    :greens in keys(obs) && export_result(obs[:greens], mp.outfile, "obs/greens"; timeseries=true)
    :etpc_plus in keys(obs) && export_result(obs[:etpc_plus], mp.outfile, "obs/etpc_plus")
    :etpc_minus in keys(obs) && export_result(obs[:etpc_minus], mp.outfile, "obs/etpc_minus")
    :zfpc_plus in keys(obs) && export_result(obs[:zfpc_plus], mp.outfile, "obs/zfpc_plus")
    :zfpc_minus in keys(obs) && export_result(obs[:zfpc_minus], mp.outfile, "obs/zfpc_minus")
    :etcdc_plus in keys(obs) && export_result(obs[:etcdc_plus], mp.outfile, "obs/etcdc_plus")
    :etcdc_minus in keys(obs) && export_result(obs[:etcdc_minus], mp.outfile, "obs/etcdc_minus")
    :zfcdc_plus in keys(obs) && export_result(obs[:zfcdc_plus], mp.outfile, "obs/zfcdc_plus")
    :zfcdc_minus in keys(obs) && export_result(obs[:zfcdc_minus], mp.outfile, "obs/zfcdc_minus")
    :zfccc_xx in keys(obs) && export_result(obs[:zfccc_xx], mp.outfile, "obs/zfccc_xx")
    :zfccc_yy in keys(obs) && export_result(obs[:zfccc_yy], mp.outfile, "obs/zfccc_yy")
    :sfdensity in keys(obs) && export_result(obs[:sfdensity], mp.outfile, "obs/sfdensity")
    :tdgfs_gt0 in keys(obs) && export_result(obs[:tdgfs_gt0], mp.outfile, "obs/tdgfs_gt0")
    :tdgfs_g0t in keys(obs) && export_result(obs[:tdgfs_g0t], mp.outfile, "obs/tdgfs_g0t")

    # meta data
    if !isnothing(nsweeps)
      h5open(mp.outfile, "r+") do fout
          HDF5.has(fout, "nsweeps") && HDF5.o_delete(fout, "nsweeps")
          HDF5.has(fout, "write_every_nth") && HDF5.o_delete(fout, "write_every_nth")
          fout["nsweeps"] = nsweeps
          fout["write_every_nth"] = mp.p.write_every_nth
      end
    end
    println("Done. Exported to $(mp.outfile).")
    flush(stdout)
end

# -------------------------------------------------------
#                 Helpers and stuff
# -------------------------------------------------------

"""
Allows one to use keyword syntax to add new entries to a `NamedTuple` `obs`.
"""
add(obs::NamedTuple; kw...) = combine(obs, values(kw))
"""
Combine two `NamedTuple`s into a new one.
"""
combine(x::NamedTuple, y::NamedTuple) = NamedTuple{(keys(x)..., keys(y)...)}((x...,y...))

"""
Create or load LightObservable
"""
function create_or_load_lightobs(abbrev::String, args...; kwargs...)
  print(abbrev, ": looking for old LightObservable in ", mp.objfile ," ... ")
  isfile(mp.objfile*"_tmp") && rm(mp.objfile*"_tmp")
  if isfile(mp.objfile)
    x = nothing
    jldopen(mp.objfile, "r") do f
      if HDF5.has(f.plain, abbrev)
        x = read(f[abbrev])
        println("found and loaded!")
      end
    end
    !isnothing(x) && return x
  end

  # couldn't be loaded -> create new LightObservable
  # println("created new lightobs")
  println("not found.")
  return LightObservable(args...; kwargs...)
end

"""
Create or load Observable
"""
function create_or_load_obs(abbrev::String, args...; kwargs...)
  print(abbrev, ": looking for old Observable in ", mp.objfile ," ... ")
  isfile(mp.objfile*"_tmp") && rm(mp.objfile*"_tmp")
  if isfile(mp.objfile)
    x = nothing
    jldopen(mp.objfile, "r") do f
      if HDF5.has(f.plain, abbrev)
        x = read(f[abbrev])
        println("found and loaded!")
      end
    end
    !isnothing(x) && return x
  end

  # couldn't be loaded -> create new LightObservable
  println("not found.")
  return Observable(args...; kwargs...)
end

nothing

max_stars_repo_name: crstnbr/dqmc

max_issues_repo_name: crstnbr/dqmc

max_forks_repo_name: crstnbr/dqmc


text:
include("Include.jl")

# Load the ensemble for AhR -
time_AhR = [0.0 48.0]
(raw_data_AhR,scaled_data_AhR) = process_ensemble_data(time_AhR,38);
writedlm("./plot/scaled_data_AhR.dat",scaled_data_AhR);

# Load the ensemble for CD38 -
time_CD38 = [0.0 24.0 48.0 72.0]
(raw_data_CD38,scaled_data_CD38) = process_ensemble_data(time_CD38,41);
writedlm("./plot/scaled_data_CD38.dat",scaled_data_CD38);

# Load the ensemble for CD11b -
time_CD11b = [0.0 24.0 48.0 72.0]
(raw_data_CD11b,scaled_data_CD11b) = process_ensemble_data(time_CD11b,39);
writedlm("./plot/scaled_data_CD11b.dat",scaled_data_CD11b);

# Load the ensemble for P47Phox -
time_P47Phox = [0.0 48.0]
(raw_data_P47Phox,scaled_data_P47Phox) = process_ensemble_data(time_P47Phox,50);
writedlm("./plot/scaled_data_P47Phox.dat",scaled_data_P47Phox);

# Load the ensemble for EGR1 -
time_EGR1 = [0.0 48.0]
(raw_data_EGR1,scaled_data_EGR1) = process_ensemble_data(time_EGR1,44);
writedlm("./plot/scaled_data_EGR1.dat",scaled_data_EGR1);

# Load the ensemble for IRF1 -
time_IRF1 = [0.0 48.0]
(raw_data_IRF1,scaled_data_IRF1) = process_ensemble_data(time_IRF1,46);
writedlm("./plot/scaled_data_IRF1.dat",scaled_data_IRF1);

# Load the ensemble for GIF1 -
time_GIF1 = [0.0 48.0]
(raw_data_GIF1,scaled_data_GIF1) = process_ensemble_data(time_GIF1,45);
writedlm("./plot/scaled_data_GIF1.dat",scaled_data_GIF1);

# Load the ensemble for PU1 -
time_PU1 = [0.0 48.0]
(raw_data_PU1,scaled_data_PU1) = process_ensemble_data(time_PU1,52);
writedlm("./plot/scaled_data_PU1.dat",scaled_data_PU1);

# Load the ensemble for OCT4 -
time_PU1 = [0.0 48.0]
(raw_data_OCT4,scaled_data_OCT4) = process_ensemble_data(time_OCT4,48);
writedlm("./plot/scaled_data_OCT4.dat",scaled_data_OCT4);

# Load the ensemble for CEBPa -
time_CEBPa = [0.0 48.0]
(raw_data_PU1,scaled_data_CEBPa) = process_ensemble_data(time_CEBPa,42);
writedlm("./plot/scaled_data_CEBPa.dat",scaled_data_CEBPa);


text:
using CSV, GLM, Distributions, Plots; pyplot()

data = CSV.read("../data/weightHeight.csv")
df = sort(data, :Weight)[1:20,:]
model = lm(@formula(Height ~ Weight), df)
pred(x) = coef(model)'*[1, x]
tStat = coef(model)[2]/stderror(model)[2]
n = size(df)[1]
pVal = 2*ccdf(TDist(n-2),tStat)
println("Manual Pval: ", pVal)
println("Manual Confidence Interval: ",
    (coef(model)[2] - quantile(TDist(n-2), 0.975)*stderror(model)[2],
     coef(model)[2] + quantile(TDist(n-2), 0.975)*stderror(model)[2]))
println(model)

scatter(df.Weight, df.Height,c=:blue, msw=0)
xlims = [minimum(df.Weight), maximum(df.Weight)]
plot!(xlims, pred.(xlims),
    c=:red, legend=:none, xlabel = "Weight (kg)", ylabel = "Height (cm)")

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
using Tar, Inflate, SHA
filename="C:\\Users\\pi96doc\\Documents\\Programming\\Java\\View5D\\View5D_v2.3.2-jar.tar.gz"
println("git-tree-sha1 = \"", Tar.tree_hash(IOBuffer(inflate_gzip(filename))),"\"")
println("sha256 = \"", bytes2hex(open(sha256, filename)), "\"")


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 34ff0254-da58-11eb-3cd1-9de133b02ae1
begin
	import Pkg
	Pkg.add(["Images","TestImages",
			"TiledIteration",
			"PlutoUI", "Plots", "Combinatorics",
			"MosaicViews", "PerceptualColourMaps", "VideoIO"])
	using Images,TestImages,LinearAlgebra,TiledIteration,PlutoUI,Plots
	using MosaicViews, PerceptualColourMaps, FFTW, VideoIO, Combinatorics
end

# ╔═╡ d238625c-7fcd-4d34-8285-6f889240248d
begin
	vid = openvideo(VideoIO.testvideo("ladybird")) do f
			collect(f)
		end
end

# ╔═╡ 340f339b-9f00-4228-bd65-0a27f58fc996
size(vid)

# ╔═╡ 55b4bd75-e696-4a3c-997a-22ad4e59aa92
begin
	function crop_to_multiple(im, factor; center::Bool = true)
		newsize = div.(size(im), factor).*factor
		offsets = div.(mod.(size(im), factor), center ? 2 : 1)
		return im[range.(offsets.+1, newsize.+offsets, step=1)...]
	end
	vidarr = crop_to_multiple(cat(vid...,dims=3),40) # video as array
end

# ╔═╡ 63c1dc77-4c8f-4afb-a220-bc06238e1d92
begin
	function array2frames(Z::Array{T,3}) where T
		collect(eachslice(Z;dims=3))
	end
	array2frames(vidarr)
end

# ╔═╡ efd07cac-cddf-474f-9cd2-3589dc62b190
begin
	tile_sizes = [5,10,20,40]
	frame_sizes = [5,10,20,40]
end

# ╔═╡ d9a53542-9a5b-4735-bc78-08b86dd35707
@bind tilesz_str Select(string.(tile_sizes) .=> string.(tile_sizes), default="10")

# ╔═╡ 1c843884-070a-41ed-aaac-c46c15bff20d
tilesz = parse(Int, tilesz_str)

# ╔═╡ 5d6be6c2-f011-49cb-9366-64aa48b72a7a
@bind framesz_str Select(string.(frame_sizes) .=> string.(frame_sizes), default="10")

# ╔═╡ 9ee7409a-9bda-4827-9630-56c0d6be3f18
framesz = parse(Int, framesz_str)

# ╔═╡ 22987ff8-c802-41d3-9a70-135fd3daf629
begin
	function tile_reshape(im,tilesz,framesz)
	  hcat(vec(map(tile -> vec(im[tile...]),
			TileIterator(axes((im)), (tilesz,tilesz,framesz))))...)
	end
	X = tile_reshape(vidarr, tilesz, framesz)
end

# ╔═╡ 3674247e-224b-485d-b981-90c0161139be
#tU,tS,tV = svd(Float64.(Gray.(X))); # this is too slow, ugh...

# ╔═╡ 0d8c3e6f-05c8-4ce2-8289-aef937afd8aa
size(X)

# ╔═╡ be92a63e-b329-46d7-a66c-07e7981e683b
Q,R = qr(Float64.(Gray.(X))');

# ╔═╡ ae872b1c-75ea-4c39-a59b-fbf66fc717ce
tVQ, tS, tU = svd(R);

# ╔═╡ 9f140f91-be19-4f03-96fc-4c82cdb14ab2
tV = Q*tVQ;

# ╔═╡ 898ca787-59ba-4775-a57d-7d7c00b3faaf
begin
	# Need to be careful here as we need to normalize things across all images!
	# Here's the code for apply diverging colour map...
#=
	    minv = minimum(img)
    maxv = maximum(img)

    if refval < minv || refval > maxv
        @warn("Reference value is outside the range of image values")
    end

    dr = maximum([maxv - refval, refval - minv])
    rnge = [-dr dr] .+ refval

    return applycolourmap(img, cmap, rnge)
	=#
	function array2frames_cmap(Z::Array{T,3}) where T
		refval = 0
		minv,maxv = extrema(Z)
		dr = maximum([maxv - refval, refval - minv])
	    rnge = [-dr dr] .+ refval

		frames = array2frames(Z)

		map(A -> colorview(RGB,
			permutedims(applycolourmap(
                        copy(A), cmap("D2"), rnge),
					   (3,1,2))), frames)
	end
	array2frames_cmap(reshape(tU[:,1], (tilesz,tilesz,framesz)))
end

# ╔═╡ 77e39090-d27f-4f3b-a2f1-fb8b8aad21ed
begin
	map(i->
		(array2frames_cmap(reshape(tU[:,i], (tilesz,tilesz,framesz)))),
		1:10)
end

# ╔═╡ 9f631d3e-d0ae-4b6e-81e8-814f604f17e9
begin

	# way more complicated than strictly necessary,
	# but this is a general solution to find the ith index
	# in a multinomial expansion

	# https://math.stackexchange.com/questions/2576122/can-stars-and-bars-method-be-used-with-constraints-on-the-number-of-stars
function stars_and_bars_upper_constraints(n::Integer,upper::NTuple{N,Int}) where N
  k = length(upper)
  #@show k
  nall = binomial(n-1,k-1)
  #@show nall
  for v in upper
    if n-1-v <= 0
      continue
    end
    nall -= binomial(n-1-v,k-1)
  end
  return nall
end

""" Find the ith tuple with sum equal to s and no element larger than the upperbounds. """
function tuple_sum_unrank(s::Integer, upper::NTuple{N,Int}, i::Integer) where N
	# TODO, if sum of upper >= s, then there is a trick...
	# https://math.stackexchange.com/questions/1388764/integer-solutions-using-pie/1388847#1388847
	# see drhab answer...
	for ta = multiexponents(N, s-N) # sum to s-N (because )
		t = NTuple{N,Int}(ta)
		t  = t.+1 # adjust to sum to s so we add one to each index
		if all(t .<= upper) # check for upper-boudns
			i -=1
			if i==0
				return t
			end
		end
	end
	return i
end

##

function idct_id(offset,sz,fsz)
  A = zeros(sz,sz,fsz)
  # find the triangular index
	for s=3:sz+sz+fsz # these are the diagonals...
		# size of the current diagonal
		ndiag = stars_and_bars_upper_constraints(s, (sz,sz,fsz))
		if offset > ndiag
			offset -= ndiag
		else
			# ideally we'd like to unrank this... but ugh...
			# instead, we are just going to enumerate and throw out
			# invalid stuff...
			#@show ndiag, offset
			t = tuple_sum_unrank(s, (sz, sz, fsz), offset)
			#@show t, typeof(t)
			@assert typeof(t) <: Tuple
			A[t...] = 1
			break
		end
	end
	return A
end

	map(i->
		array2frames_cmap(idct(idct_id(i,tilesz,framesz))),
		1:10)
end

# ╔═╡ 3182eb0a-b5c8-4292-ab9f-201473769637
map(i->array2frames_cmap(dct(reshape(tU[:,i], (tilesz,tilesz,framesz)))), 1:10)

# ╔═╡ 40cd74a9-791b-44ab-81ee-940d24fe2e4b
begin
	# Ahh, I stole this from a github issue/pull request for ...
	# https://github.com/JuliaStats/MultivariateStats.jl/pull/130/files
	## Comparison to other implementations
	# The implementation of varimax in R row-normlises the input matrix before
	# application of the algorithm and rescales the rows afterwards.
	## Reference
	# Mikkel B. Stegmann, Karl Sjöstrand, Rasmus Larsen, "Sparse modeling of
	# landmark and texture variability using the orthomax criterion,"
	# Proc. SPIE 6144, Medical Imaging 2006: Image Processing, 61441G
	# (10 March 2006); doi: 10.1117/12.651293
	function orthomax(F::AbstractMatrix, γ, miniter, maxiter, ϵ)
		n, p = size(F)
		if n < 2
			return (F, Matrix{eltype(F)}(I, p, p))
		end

		# Warm up step
		# Do one step. If that first step did not lead away from the identity
		# matrix enough use a random orthogonal matrix as a starting point.
		M = svd(F' * (F .^ 3 - γ / n * F * Diagonal(vec(sum(F .^ 2, dims=1)))))
		R = M.U * M.Vt
		if norm(R - Matrix{eltype(R)}(I, p, p)) < ϵ
			R = qr(randn(p, p)).Q
		end

		# Main iterations
		d = 0
		lastchange = NaN
		converged = false
		for i in 1:maxiter
			dold  = d
			B = F * R
			M = svd(F' * (B .^ 3 - γ / n * B * Diagonal(vec(sum(B .^ 2, dims=1)))))
			R = M.U * M.Vt
			d = sum(M.S)
			lastchange = abs(d - dold) / d
			if lastchange < ϵ && i >= miniter
				converged = true
				break
			end
		end

		converged || throw(ConvergenceException(maxiter, lastchange, ϵ))

		(F * R, R)
	 end
	tU2 = orthomax(tU[:,1:32], 1.0, 20, 1000, 1e-12)[1]
end

# ╔═╡ f3fbac2b-3582-4404-9391-a2099c97e124
begin
	map(i->
		(array2frames_cmap(reshape(tU2[:,i], (tilesz,tilesz,framesz)))),
		1:size(tU2,2))
end

# ╔═╡ 7d617284-7779-470a-821a-296572f00578
begin
	rtile = 3
	function tiled_svd_approx(r::Int, tilesz::Int, framesz::Int; tU, tS, tV)
		Y = tU[:,1:r]*Diagonal(tS[1:r])*tV[:,1:r]'
		B = similar(vidarr, Float64)
		for (i,tile) in enumerate(TileIterator(axes(vidarr), (tilesz,tilesz,framesz)))
		  B[tile...] .= reshape(Y[:,i], (tilesz,tilesz,framesz))
		end
		return B, length(tU[:,1:r])+length(tV[:,1:r])
	end
	tile_svd_approx_gray, nparam_tile_gray = tiled_svd_approx(rtile,tilesz,framesz; tU, tS, tV)
	# save a video
	frames = map(M -> Gray{N0f8}.(clamp01nan.(M)), array2frames((tile_svd_approx_gray)))
end

# ╔═╡ 6d40196d-9f4a-4baa-b0e8-d8f43244ef7d
VideoIO.save("ladybird-grey-tile-$(tilesz)-frame-$(framesz)-rank-$rtile.mp4", frames)

# ╔═╡ 77cfd1a7-2326-47da-b710-9aa443aa62cc
frames

# ╔═╡ Cell order:
# ╠═34ff0254-da58-11eb-3cd1-9de133b02ae1
# ╠═d238625c-7fcd-4d34-8285-6f889240248d
# ╠═340f339b-9f00-4228-bd65-0a27f58fc996
# ╠═55b4bd75-e696-4a3c-997a-22ad4e59aa92
# ╠═63c1dc77-4c8f-4afb-a220-bc06238e1d92
# ╠═efd07cac-cddf-474f-9cd2-3589dc62b190
# ╠═d9a53542-9a5b-4735-bc78-08b86dd35707
# ╠═1c843884-070a-41ed-aaac-c46c15bff20d
# ╠═5d6be6c2-f011-49cb-9366-64aa48b72a7a
# ╠═9ee7409a-9bda-4827-9630-56c0d6be3f18
# ╠═22987ff8-c802-41d3-9a70-135fd3daf629
# ╠═3674247e-224b-485d-b981-90c0161139be
# ╠═0d8c3e6f-05c8-4ce2-8289-aef937afd8aa
# ╠═be92a63e-b329-46d7-a66c-07e7981e683b
# ╠═ae872b1c-75ea-4c39-a59b-fbf66fc717ce
# ╠═9f140f91-be19-4f03-96fc-4c82cdb14ab2
# ╠═898ca787-59ba-4775-a57d-7d7c00b3faaf
# ╠═77e39090-d27f-4f3b-a2f1-fb8b8aad21ed
# ╠═9f631d3e-d0ae-4b6e-81e8-814f604f17e9
# ╠═3182eb0a-b5c8-4292-ab9f-201473769637
# ╠═40cd74a9-791b-44ab-81ee-940d24fe2e4b
# ╠═f3fbac2b-3582-4404-9391-a2099c97e124
# ╠═7d617284-7779-470a-821a-296572f00578
# ╠═6d40196d-9f4a-4baa-b0e8-d8f43244ef7d
# ╠═77cfd1a7-2326-47da-b710-9aa443aa62cc

max_stars_repo_name: dgleich/better-svd-for-images

max_issues_repo_name: dgleich/better-svd-for-images

max_forks_repo_name: dgleich/better-svd-for-images


text:
# Bucket all tables with the same column counts together in anticipation of get index
# overloading. This will minimize the number of overloaded functions added to the system.
abstract AbstractDataTable{N, R<:NTuple{N}, C<:NTuple{N, AbstractVector}}

# Tuple unroller: Un-tuple-ize a vector of tuples into a tuple of vectors
@generated function Base.convert{R<:Tuple}(::Type{Tuple{[Vector{r} for r in R.parameters]...}}, vs::AbstractVector{R})
	x = :(t = (); for i = 1:length(vs) end; t;)
	x.args[1].args[2].args = [:(Vector{$(r)}(length(vs))) for r in R.parameters][:]
	x.args[2].args[1].args[2].args = [:(@inbounds t[$(r)][i] = vs[i][$(r)]) for r = 1:length(R.parameters)][:]
	x
end

# Vector unroller: index a tuple of vectors into a tuple
@generated function Base.getindex{N}(D::AbstractDataTable{N}, r::Int64)
	x = :(())
	prepend!(x.args, [:(Base.getindex(getcolumns(D)[$(c)], r)) for c = 1:N])
	x
end

# Unroll value setting
@generated function Base.setindex!{N, R, C}(D::AbstractDataTable{N, R, C}, v::R, r::Int64)
	x = quote v end
	prepend!(x.args, [:(Base.setindex!(getcolumns(D)[$(c)], v[$(c)], r)) for c = 1:N])
	x
end

# Unroll vector emptying
@generate function Base.empty!{N}(D::AbstractDataTable{N})
	x = quote D end
	prepend!(x.args, [:(Base.empty!(getcolumns(D)[$(c)])) for c = 1:N])
	x
end

# All concrete types must implement getcolumns that returns NTuple{N, AbstractVector}
# getcolumns()

# Courtesy method to check that the table has been constructed properly
Base.isvalid{N, R, C}(D::AbstractDataTable{N, R, C}) =
	length(getcolumns(D)) == N &&
	length(R.parameters) == N &&
	length(C.parameters) == N &&
	minimum([length(getcolumns(D)[c]) for c = 1:N]) == maximum([length(getcolumns(D)[c]) for c = 1:N]) &&
	typeof(getcolumns(D)) == C &&
	all([(R.parameters[i] == C.parameters[i].parameters[1])::Bool for i = 1:N]) &&
	isvalid(Val{DataView}, D)

# Iterable methods
Base.start(::AbstractDataTable) = 1
Base.next(D::AbstractDataTable, r::Int64) = (D[r], r + 1)
Base.done(D::AbstractDataTable, r::Int64) = r > length(D)
Base.eltype{N, R, C}(::Type{AbstractDataTable{N, R, C}}) = R
Base.length{N}(D::AbstractDataTable{N}) = minimum([length(getcolumns(D)[c]) for c = 1:N])

# Collection methods
Base.isempty(D::AbstractDataTable) = length(D) < 1
Base.endof(D::AbstractDataTable) = length(D)

# Indexable methods
# A single column returns a vector; type unstable
# A single row and column pair returns an element; type unstable
Base.getindex{N, R, C}(D::AbstractDataTable{N, R, C}, c::Symbol) = getcolumns(D)[getcolumnindex(D, c)]::C.parameters[c]
Base.getindex{N, R, C}(D::AbstractDataTable{N, R, C}, c::AbstractString) = D[Symbol(c)]::C.parameters[c]
Base.getindex{N, R, C}(D::AbstractDataTable{N, R, C}, r::Int64, c::Int64) = getcolumns(D)[c][r]::R.parameters[c]
Base.getindex{N, R, C}(D::AbstractDataTable{N, R, C}, r::Int64, c::Symbol) = D[r, getcolumnindex(D, c)]::R.parameters[c]
Base.getindex{N, R, C}(D::AbstractDataTable{N, R, C}, r::Int64, c::AbstractString) = D[r, Symbol(c)]::R.parameters[c]

# Array methods
Base.ndims(::AbstractDataTable) = 2
Base.size{N}(D::AbstractDataTable{N}) = (length(D), N)
Base.linearindexing(::Type{AbstractDataTable}) = Base.LinearFast()


text:
function z2gk(A::Float64, B::Float64, g::Float64, k::Float64, z::Float64, c=0.8)
    temp = exp(-g*z)
    A + B*(1.0+c*(1.0-temp)/(1.0+temp))*(1.0+z^2)^k*z
end

function rgk(pars)
    (A,B,g,k) = pars
    z = randn()
    z2gk(A, B, g, k, z)
end

##Simulates U(0,1) order statistics specified in "orderstats" from n total sims
##orderstats should be in ascending order
##See Ripley "Stochastic Simulation" pg 98
function unif_os(orderstats::Array{Int,1}, n::Int)
    p = size(orderstats)[1]
    w = Array{Float64}(p+1)
    w[1] = rand(Gamma(orderstats[1]))
    for i in 2:p
        w[i] = rand(Gamma(orderstats[i] - orderstats[i-1]))
    end
    w[p+1] = rand(Gamma(n + 1 - orderstats[p]))
    wsums = cumsum(w)
    wsums[1:p] / wsums[p+1]
end

##Efficiently simulates g&k order statistics specified in "orderstats" from n total sims
##orderstats should be in ascending order
function rgk_os(pars::Array{Float64,1}, orderstats::Array{Int,1}, n::Int)
    (A,B,g,k) = pars
    u = unif_os(orderstats, n)
    z = quantile.(Normal(), u)
    map((x)->z2gk(A,B,g,k,x), z)
end


text:
include("solver_preamble.jl")
using COSMO
factory = with_optimizer(COSMO.Optimizer)
config = MOI.Test.TestConfig(atol=1e-4, rtol=1e-4)
@testset "Linear" begin
    Tests.linear_test(factory, config)
end
@testset "SOC" begin
    Tests.soc_test(factory, config)
end
@testset "SDP" begin
    Tests.sd_test(factory, config)
end


text:
function revers(side :: HorizonSide)
    side = HorizonSide((Int(side)+2)%4)
    return side
end

function movements_to_side!(r, side :: HorizonSide)
    while isborder(r,side) == false
        move!(r,side)
    end
end

function movements_find_temp!(r,side::HorizonSide, side_part :: HorizonSide)

    temp = 0
    count = 0

    while isborder(r,side_part) == false
        while isborder(r,side) ==  false
            if ismarker(r)
                temp = temp + temperature(r)
                count = count + 1
            end
            move!(r,side)
        end
        move!(r, side_part)
        side = revers(side)

        if ismarker(r)
            temp = temp + temperature(r)
            count = count + 1
        end

    end

    while isborder(r,side) ==  false
        if ismarker(r)
            temp = temp + temperature(r)
            count = count + 1
        end
        move!(r,side)
    end

    if ismarker(r)
        temp = temp + temperature(r)
        count = count + 1
    end

    return temp, count

end

function mid_temp(r)

    mid_temp = 0
    temp = 0
    count = 0

    temp, count = movements_find_temp!(r,Ost,Nord)

    movements_to_side!(r,Sud) #вернул его домой для порядка
    movements_to_side!(r,West)

    println(count)

    mid_temp = temp / count

    print(mid_temp)
end

mid_temp(r)
show(r)

max_stars_repo_name: AntonDudDm/mirea-progs

max_issues_repo_name: AntonDudDm/mirea-progs

max_forks_repo_name: AntonDudDm/mirea-progs


text:
workspace()

using IntervalArithmetic
using EAGOBranchBound

X = [Interval(1,5)]
model = BnBModel(X)
solver = BnBSolver()
set_to_default!(solver)
solver.opt = [x -> x[1]^2]
solver.Lower_Prob = (X,k,p,opt,temp) -> [opt[1](X).lo,mid.(X),true,[]]
solver.Upper_Prob = (X,k,p,opt,temp) -> [opt[1](X).hi,mid.(X),true,[]]
a = solver.Lower_Prob(X,1,1,solver.opt,[])
b = solver.Upper_Prob(X,1,1,solver.opt,[])
solveBnB!(solver,model)

solver.iter_lim = true
solver.max_iter = 10
solver.Lower_Prob = (X,k,p,opt,temp) -> [opt[1](X).lo,mid.(X),false,[]]
model2 = BnBModel(X)
solveBnB!(solver,model2)
c1 = getobjval(model2)

solver.Lower_Prob = (X,k,p,opt,temp) -> [opt[1](X).lo,mid.(X),true,[]]
solver.Upper_Prob = (X,k,p,opt,temp) -> [opt[1](X).hi,mid.(X),false,[]]
model3 = BnBModel(X)
solveBnB!(solver,model3)
c2 = getobjval(model3)

solver.Lower_Prob = (X,k,p,opt,temp) -> [opt[1](X).lo,mid.(X),true,[]]
solver.Upper_Prob = (X,k,p,opt,temp) -> [((k==2) ? opt[1](X).hi : Inf),mid.(X),false,[]]
model4 = BnBModel(X)
solveBnB!(solver,model4)
c3 = getobjval(model4)
println("c3: $c3")
#@test c3 == Inf

solver.Lower_Prob = (X,k,p,opt,temp) -> [opt[1](X).lo,mid.(X),true,[]]
solver.Upper_Prob = (X,k,p,opt,temp) -> [opt[1](X).hi,mid.(X),true,[]]
solver.Repeat_Chk = (S,M,X1,X2) -> (M.lbcnt == 3)
model5 = BnBModel(X)
solveBnB!(solver,model5)
c5 = getobjval(model5)
println("c5: $c5")
#@test c4 == Inf


text:
# A Simple Computational Algebra System

_mul(::Val{X}, ::Val{Y}) where {X, Y} = Val(X*Y)
_mul(x::Number, y::Number) = x * y
_neg(::Val{X}) where X = Val(-X)
_neg(x::Number) = -x

# negate
Base.:(-)(x::AbstractBlock) = Scale(Val(-1), x)
Base.:(-)(x::Scale{Val{-1}}) = content(x)
Base.:(-)(x::Scale) = Scale(_neg(x.alpha), content(x))

# scaler multiply block
Base.:(*)(α::T, x::AbstractBlock) where {T<:Union{Val, Number}} = Scale(α, x)
Base.:(*)(α::T, x::Scale{T2}) where {T<:Val,T2<:Val} = Scale(_mul(α, x.alpha), content(x))
Base.:(*)(α::T, x::Scale{T2}) where {T<:Number,T2<:Number} = Scale(_mul(α, x.alpha), content(x))
Base.:(*)(x::AbstractBlock, α::Union{Val, Number}) = α * x

# block multiply block
Base.:(*)(x::Scale, y::Scale) = (_mul(x.alpha, y.alpha)) * (content(x) * content(y))
Base.:(*)(x::Scale, y::AbstractBlock) = x.alpha * chain(y, content(x))
Base.:(*)(y::AbstractBlock, x::Scale) = x.alpha * chain(content(x), y)
Base.:(*)(xs::AbstractBlock...) = chain(Iterators.reverse(xs)...)

# add
Base.:(+)(xs::AbstractBlock...) = Add(xs...)

# div
Base.:(/)(A::AbstractBlock, x::Number) = (1 / x) * A

# reduce
Base.prod(blocks::AbstractVector{<:AbstractBlock{D}}) where D =
    chain(Iterators.reverse(blocks)...)
Base.sum(blocks::AbstractVector{<:AbstractBlock{D}}) where D = +(blocks...)

#sub
Base.:(-)(lhs::AbstractBlock, rhs::AbstractBlock) = Add(lhs, -rhs)

# pow
Base.:(^)(x::AbstractBlock, n::Int) = chain((copy(x) for k = 1:n)...)


text:
# basic definitions
@enum Sourcetype line=1 plane=2

@enum Boundarytype rigid=1 free=2 absorbing=3

struct Displacementfield
    v::Array{Float64, 3}
    zeta::Array{Float64, 2}
end

struct Boundarydisplacements
    vt::Array{Float64, 2}
    vb::Array{Float64, 2}
    vl::Array{Float64, 2}
    vr::Array{Float64, 2}
end

struct Boundaries
    top::Boundarytype
    bottom::Boundarytype
    left::Boundarytype
    right::Boundarytype
end

struct Params
    delta::Float64
    #LX,LZ: number of grid points in x and z direction (note: boundaries are at points with indices 2, lx-1,lz-1)
    lx::Int32
    lz::Int32
    #vp, eps, kappa: background velocity, std.dev. of random fluctuations, parameter for density fluctuations (3F10.0)
    vp::Float64 #//[vp]=km/s
    eps::Float64 # 0.05
    kappa::Float64 # 0.0;
    #source signal (KUEPPERS,1958):
    nqx::Int32 # NQX: x index of source, last config: 2000-pt lattice 1000
    nqz::Int32 # NQZ: z index of source, last config: 2000-pt lattice 100
    iquell::Sourcetype #IQUELL:=1: line source,  =2: plane wave starting at z=0  (F10.0,4I10)
    dx::Float64 # 0.12;//[dx]=km
    nsig::Int32 #2; //NSIG: number of exrema
    nmax::Int32 #16000; ?
end

function vmax()
end

function deltat(p)
    p.dx/sqrt(2.0)/vmax()
end

# Einganspuls (Kuepperwavelet)
function signal(t::Float64, delta::Float64=6.0, nsig::Int=2)
    ampmax = zeros(10)
    ampmax[1]=1.333
    ampmax[2]=1.299
    ampmax[3]=1.6
    ampmax[4]=1.585
    ampmax[5]=1.714
    ampmax[6]=1.706
    ampmax[7]=1.778
    ampmax[8]=1.773
    ampmax[9]=1.818
    ampmax[10]=1.815
    fm = (1.0 * nsig + 2.0) / (1.0 * nsig)
    del = nsig * pi / delta
    signal = (t < 0 || t > delta) ? 0.0 : (sin(del * t)-((sin(fm* del *t))/fm)) / ampmax[nsig]
end

function timeloop(p::Params, bd::Boundarydisplacements, d::Displacementfield)
    for nt in 1:p.nmax
        saveboundarydisplacements(p, bd, d)
    end
end

function setinitialcondtions(field::Displacementfield, p::Params)
    if p.iquell == 2
        for j in 2:p.lz-1
            arg = -abs(p.nqz)/p.nqz*(j-p.nqz)*p.dx/p.vp
            s1 = signal(arg-deltat(p), p.delta, p.nsig)
            s2 = signal(arg, p.delta, p.nsig)
            for i in 1:p.lx)
                field.v[i][j][1] = s1
                field.v[i][j][2] = s2
            end
        end
    end
end

function saveboundarydisplacements!(p::Params, bd::Boundarydisplacements, d::Displacementfield)
    for i in 2:p.lx-1
        bd.vt[i][1] = d.v[i][2][1]
        bd.vt[i][2] = d.v[i][3][1]
        bd.vb[i][1] = d.v[i][p.lz-1][1]
        bd.vb[i][2] = d.v[i][p.lz-2][1]
    end
    for j in 2:p.lz-1
        bd.vl[j][1] = d.v[2][j][1]
        bd.vl[j][2] = d.v[3][j][1]
        bd.vr[j][1] = d.v[p.lx-1][j][1]
        bd.vr[j][2] = d.v[p.lx-2][j][1]
    end
end

function calculatedisplacement!(p::Params, d::Displacementfield)
    for i in 2:p.lx-1
        for j in 2:p.lz-1
            r5 = 1. / (1. + p.kappa * fij)
            f1 = 1. + r5 * (1. + p.kappa * d.zeta[i + 1][j])
            f2 = 1. + r5 * (1. + p.kappa * d.zeta[i - 1][j])
            f3 = 1. + r5 * (1. + p.kappa * d.zeta[i][j+1])
            f4 = 1. + r5 * (1. + p.kappa * d.zeta[i][j-1])
            v[i][j][1] = -d.v[i][j][1] + fact1 * (1. + d.zeta[i][j]) * (1. + d.zeta[i][j]) *
                ((d.v[i+1][j][2] - d.v[i][j][2]) / f1 - (d.v[i][j][2] - d.v[i-1][j][2]) / f2 +
                 (d.v[i][j+1][2] - d.v[i][j][2]) / f3 - (d.v[i][j][2] - d.v[i][j-1][2]) / f4) + 2. * d.v[i][j][2]
        end
    end
end

function applyboundaryconditions!()
end


text:
###################

using CryptoUtils: sqrt_mod_prime

# Field element to integer. Does not make sense for binary curves

<|(::Type{AffinePoint{EQ, F}}, args::Tuple{F, F}) where {EQ <: EllipticCurve, F <: Field} = AffinePoint{EQ, F}(args...)

function <|(::Type{P}, bytes::Tuple{Vector{UInt8}, Vector{UInt8}}) where P <: AbstractPoint
    xb, yb = bytes
    F = field(P)
    P <| (F <| xb, F <| yb)
end

function <|(::Type{P}, compressed_point::Tuple{Vector{UInt8}, UInt8}) where P <: AbstractPoint

    xb, ỹb = compressed_point

    F = field(P)

    x = F <| xb
    ỹ = Bool(ỹb)

    return P <| (x, ỹ)
end

#

function <|(::Type{AffinePoint{EQ, F}}, compressed_point::Tuple{BigInt, Bool}) where {EQ <: Weierstrass, F <: PrimeField}

    P = AffinePoint{EQ, F}

    x, ỹ = compressed_point

    p = modulus(F)

    α = mod(x^3 + a(EQ) * x + b(EQ), p)

    β = sqrt_mod_prime(α, p)

    # Rightmost bit of β. Perhaps this works:
    # Alternativelly I could do circshift
    if mod(β, 2) == ỹ
        y = β
    else
        y = p - β
    end

    return P <: (x, y)
end

function int2octet(x::Integer)

    hex = string(x, base=16)
    if mod(length(hex), 2) != 0
        hex = string("0", hex)
    end

    return hex2bytes(hex)
end

octet2int(x::Vector{UInt8}) = parse(BigInt, bytes2hex(x), base=16)

<|(::Type{Vector{UInt8}}, x::PrimeField) = int2octet(value(x))
<|(::Type{F}, x::Vector{UInt8}) where F <: PrimeField = F <| octet2int(x)

function octet2bits(x::Vector{UInt8})
    bv = BitVector(u << -i % Bool for u in x for i in 7:-1:0)
    return bv
end

octet2bits(x::Vector{UInt8}, N::Int) = octet2bits(x)[end - N + 1:end]

function bits2uint8(x::BitVector)

    s = UInt8(0)

    for (i, j) in enumerate(reverse(x))

        if j == 1
            s += UInt8(2)^(i - 1)
        end
    end

    return s
end

function bits2octet(_x::BitVector)

    x = copy(_x)

    if mod(length(x), 8) != 0

        padding = BitVector(0 for i in 1:(8 - mod(length(x), 8)))
        prepend!(x, padding)

    end

    # For now assuming that x is in the length of octets

    N = div(length(x), 8, RoundUp)

    b = reshape(x, 8, N)

    bytes = UInt8[bits2uint8(b[:, i]) for i in 1:N]

    return bytes
end

<|(::Type{Vector{UInt8}}, x::BinaryField) = bits2octet(tobits(x))
<|(::Type{F}, x::Vector{UInt8}) where F <: BinaryField = F <| octet2bits(x, bitlength(F))

function <|(::Type{P}, po::Vector{UInt8}) where P <: AbstractPoint

    pc = po[1]

    if pc == 0
        return zero(P)
    elseif pc == 2
        ỹ = false
        x = po[2:end]
        return P <| (x, ỹ)
    elseif pc == 3
        ỹ = true
        x = po[2:end]
        return P <| (x, ỹ)
    elseif pc == 4 || pc == 6 || pc == 7
        l = div(length(po) - 1, 2)
        x = po[2:(l+1)]
        y = po[(l+2):(2l+1)]
        # Can write assertion for the rightmost bit of y here

        #@infiltrate
        return P <| (x, y)
    else
        error("Wrong PC: $pc")
    end

end

function <|(::Type{Vector{UInt8}}, p::AbstractPoint; option::Symbol=:uncompressed)

    if option == :uncompressed
        x = Vector{UInt8} <| gx(p)
        y = Vector{UInt8} <| gy(p)

        return UInt8[4, x..., y...]
    else
        error("To be implemented")
    end
end

function <|(x::Tuple{DataType, Symbol}, y) where T
    return <|(x[1], y; option=x[2])
end


text:
using Requests, OAuth, Compat

type TwStream
  twstream::Requests.ResponseStream

  function TwStream(auth::TwAuth, options::Dict)
    new(init_stream(auth, "POST", "https://stream.twitter.com/1.1/statuses/filter.json", options))
  end

  TwStream(options::Dict) = TwStream(JuliaTwitter.auth, options)
end

#Initialize a Requests.ResponseStream
#Getting 401 errors? make sure to use post_streaming (typ = "POST")
function init_stream(auth::TwAuth, typ::ASCIIString, endpoint::AbstractString, options::Dict)
  #Validate typ
  typ = uppercase(typ)
  if typ != "POST" && typ != "GET"
    throw(ArgumentError("$typ is an invalid type. Use GET or POST"))
  end

  #Build query string
  query_str = Requests.format_query_str(options)

  #Build oauth_header
  oauth_header_val = oauth_header(typ, endpoint, options, auth.consumer_key, auth.consumer_secret, auth.access_token, auth.access_secret)

  #Make request
  headers = @compat(
     Dict{AbstractString,AbstractString}(
         "Content-Type" => "application/x-www-form-urlencoded",
         "Authorization" => oauth_header_val,
         "Accept" => "*/*"
     )
  )

  try
    #init stream
    #Requests.ResponseStream{MbedTLS.SSLContext}
    if typ == "POST"
      return Requests.post_streaming(URI("$(endpoint)?$query_str"); headers = headers)
    else
      return Requests.get_streaming(URI("$(endpoint)?$query_str"); headers = headers)
    end

  catch stream_error
    println("[ERROR] Cannot connect to Twitter")
    println("        See message below.")
    throw(stream_error)
  end

end

max_stars_repo_name: nick-paul/JuliaTwitter

max_issues_repo_name: nick-paul/JuliaTwitter

max_forks_repo_name: nick-paul/JuliaTwitter


text:
module SKPreprocessors

using PyCall

# standard included modules
using DataFrames
using Random
using ..AbsTypes
using ..Utils

import ..AbsTypes: fit, fit!, transform, transform!
export fit, fit!, transform, transform!
export SKPreprocessor, skpreprocessors

const preprocessor_dict = Dict{String,PyObject}()
const DEC  = PyNULL()
const FS   = PyNULL()
const IMP  = PyNULL()
const PREP = PyNULL()

function __init__()
   copy!(DEC , pyimport_conda("sklearn.decomposition","scikit-learn"))
   copy!(FS  , pyimport_conda("sklearn.feature_selection","scikit-learn"))
   copy!(IMP , pyimport_conda("sklearn.impute","scikit-learn"))
   copy!(PREP, pyimport_conda("sklearn.preprocessing","scikit-learn"))

   # Available scikit-learn learners.
   preprocessor_dict["DictionaryLearning"]          = DEC
   preprocessor_dict["FactorAnalysis"]              = DEC
   preprocessor_dict["FastICA"]                     = DEC
   preprocessor_dict["IncrementalPCA"]              = DEC
   preprocessor_dict["KernelPCA"]                   = DEC
   preprocessor_dict["LatentDirichletAllocation"]   = DEC
   preprocessor_dict["MiniBatchDictionaryLearning"] = DEC
   preprocessor_dict["MiniBatchSparsePCA"]          = DEC
   preprocessor_dict["NMF"]                         = DEC
   preprocessor_dict["PCA"]                         = DEC
   preprocessor_dict["SparsePCA"]                   = DEC
   preprocessor_dict["SparseCoder"]                 = DEC
   preprocessor_dict["TruncatedSVD"]                = DEC
   preprocessor_dict["dict_learning"]               = DEC
   preprocessor_dict["dict_learning_online"]        = DEC
   preprocessor_dict["fastica"]                     = DEC
   preprocessor_dict["non_negative_factorization"]  = DEC
   preprocessor_dict["sparse_encode"]               = DEC
   preprocessor_dict["GenericUnivariateSelect"]     = FS
   preprocessor_dict["SelectPercentile"]            = FS
   preprocessor_dict["SelectKBest"]                 = FS
   preprocessor_dict["SelectFpr"]                   = FS
   preprocessor_dict["SelectFdr"]                   = FS
   preprocessor_dict["SelectFromModel"]             = FS
   preprocessor_dict["SelectFwe"]                   = FS
   preprocessor_dict["RFE"]                         = FS
   preprocessor_dict["RFECV"]                       = FS
   preprocessor_dict["VarianceThreshold"]           = FS
   preprocessor_dict["chi2"]                        = FS
   preprocessor_dict["f_classif"]                   = FS
   preprocessor_dict["f_regression"]                = FS
   preprocessor_dict["mutual_info_classif"]         = FS
   preprocessor_dict["mutual_info_regression"]      = FS
   preprocessor_dict["SimpleImputer"]               = IMP
   preprocessor_dict["MissingIndicator"]            = IMP
   preprocessor_dict["Binarizer"]                   = PREP
   preprocessor_dict["FunctionTransformer"]         = PREP
   preprocessor_dict["KBinsDiscretizer"]            = PREP
   preprocessor_dict["KernelCenterer"]              = PREP
   preprocessor_dict["LabelBinarizer"]              = PREP
   preprocessor_dict["LabelEncoder"]                = PREP
   preprocessor_dict["MultiLabelBinarizer"]         = PREP
   preprocessor_dict["MaxAbsScaler"]                = PREP
   preprocessor_dict["MinMaxScaler"]                = PREP
   preprocessor_dict["Normalizer"]                  = PREP
   preprocessor_dict["OneHotEncoder"]               = PREP
   preprocessor_dict["OrdinalEncoder"]              = PREP
   preprocessor_dict["PolynomialFeatures"]          = PREP
   preprocessor_dict["PowerTransformer"]            = PREP
   preprocessor_dict["QuantileTransformer"]         = PREP
   preprocessor_dict["RobustScaler"]                = PREP
   preprocessor_dict["StandardScaler"]              = PREP
   #"IterativeImputer" => IMP.IterativeImputer,
   #"KNNImputer" => IMP.KNNImputer,
   #"add_dummy_feature" => PREP.add_dummy_feature,
   #"binarize" => PREP.binarize,
   #"label_binarize" => PREP.label_binarize,
   #"maxabs_scale" => PREP.maxabs_scale,
   #"minmax_scale" => PREP.minmax_scale,
   #"normalize" => PREP.normalize,
   #"quantile_transform" => PREP.quantile_transform,
   #"robust_scale" => PREP.robust_scale,
   #"scale" => PREP.scale,
   #"power_transform" => PREP.power_transform
end

"""
    SKPreprocessor(preprocessor::String,args::Dict=Dict())

A wrapper for Scikitlearn preprocessor functions.
Invoking `skpreprocessors()` will list the acceptable
and supported functions. Please check Scikitlearn
documentation for arguments to pass.

Implements `fit!` and `transform!`.
"""
mutable struct SKPreprocessor <: Transformer
   name::String
   model::Dict{Symbol,Any}

   function SKPreprocessor(args=Dict())
      default_args=Dict(
         :name => "skprep",
         :preprocessor => "PCA",
         :autocomponent=>false,
         :impl_args => Dict()
      )
      cargs = nested_dict_merge(default_args, args)
      cargs[:name] = cargs[:name]*"_"*randstring(3)
      prep = cargs[:preprocessor]
      if !(prep in keys(preprocessor_dict))
         println("$prep is not supported.")
         println()
         skpreprocessors()
         error("Argument keyword error")
      end
      new(cargs[:name],cargs)
   end
end

function SKPreprocessor(prep::String,args::Dict)
  SKPreprocessor(Dict(:preprocessor => prep,:name=>prep,args...))
end

function SKPreprocessor(prep::String; args...)
   SKPreprocessor(Dict(:preprocessor => prep,:name=>prep,:impl_args=>Dict(pairs(args))))
end

function (skp::SKPreprocessor)(;objargs...)
   skp.model[:impl_args] = Dict(pairs(objargs))
   prepname = skp.model[:preprocessor]
   skobj = getproperty(preprocessor_dict[prepname],prepname)
   newskobj = skobj(;objargs...)
   skp.model[:skpreprocessor] = newskobj
   return skp
end

function skpreprocessors()
  processors = keys(preprocessor_dict) |> collect |> x-> sort(x,lt=(x,y)->lowercase(x)<lowercase(y))
  println("syntax: SKPreprocessor(name::String, args::Dict=Dict())")
  println("where *name* can be one of:")
  println()
  [print(processor," ") for processor in processors]
  println()
  println()
  println("and *args* are the corresponding preprocessor's initial parameters.")
  println("Note: Please consult Scikitlearn's online help for more details about the preprocessor's arguments.")
end

function fit!(skp::SKPreprocessor, x::DataFrame, y::Vector=[])::Nothing
   features = x |> Array
   impl_args = copy(skp.model[:impl_args])
   autocomp = skp.model[:autocomponent]
   if autocomp == true
      cols = ncol(x)
      ncomponents = 1
      if cols > 0
         ncomponents = round(sqrt(cols),digits=0) |> Integer
         push!(impl_args,:n_components => ncomponents)
      end
   end
   preprocessor = skp.model[:preprocessor]
   py_preprocessor = getproperty(preprocessor_dict[preprocessor],preprocessor)

   # Train model
   preproc = py_preprocessor(;impl_args...)
   preproc.fit(features)
   skp.model[:skpreprocessor] = preproc
   skp.model[:impl_args] = impl_args
   return nothing
end

function fit(skp::SKPreprocessor, x::DataFrame, y::Vector=[])::SKPreprocessor
   fit!(skp,x,y)
   return deepcopy(skp)
end

function transform!(skp::SKPreprocessor, x::DataFrame)::DataFrame
   features = deepcopy(x) |> Array
   model=skp.model[:skpreprocessor]
   return collect(model.transform(features)) |> x->DataFrame(x,:auto)
end

transform(skp::SKPreprocessor, x::DataFrame)::DataFrame = transform!(skp,x)

end


text:
delays = 0:50
runs = 10
runs2 = 10
T_steps = 50

println("\nTesting MCDTS complete tree, Lorenz63 univariate:")
@time begin
@testset "MCDTS single rollout on univariate data" begin

    # L
    Random.seed!(1234)
    pecuzal = MCDTS.PecuzalOptim()
    tree = mcdts_embedding(Dataset(data[:,1]), pecuzal, w1, delays, runs)
    best_node = MCDTS.best_embedding(tree)
    @test best_node.τs == [0, 9, 42, 20]
    @test -0.94 < L(best_node) < -0.929

    # L with tws
    Random.seed!(1234)
    tws = 2:4:delays[end]
    optmodel2 = MCDTS.MCDTSOptimGoal(MCDTS.L_statistic(0,3,tws), MCDTS.Continuity_function())
    tree2 = mcdts_embedding(Dataset(data[:,1]), optmodel2, w1, delays, runs)
    best_node2 = MCDTS.best_embedding(tree2)
    @test best_node2.τs == best_node.τs
    @test L(best_node2) > L(best_node)

    # L with tws and less fiducials for computation
    Random.seed!(1234)
    tws = 2:2:delays[end]
    samplesize = 0.5
    optmodel3 = MCDTS.MCDTSOptimGoal(MCDTS.L_statistic(0,3,tws,samplesize), MCDTS.Continuity_function())
    tree3 = mcdts_embedding(Dataset(data[:,1]), optmodel3, w1, delays, runs)
    best_node3 = MCDTS.best_embedding(tree3)
    @test best_node3.τs == best_node.τs
    @test L(best_node) - 0.1 < L(best_node3) < L(best_node) + 0.1

    # L with tws and threshold
    Random.seed!(1234)
    optmodel4 = MCDTS.MCDTSOptimGoal(MCDTS.L_statistic(-0.5,3,tws), MCDTS.Continuity_function())
    tree4 = mcdts_embedding(Dataset(data[:,1]), optmodel4, w1, delays, runs)
    best_node4 = MCDTS.best_embedding(tree4)
    @test length(best_node4.τs) < length(best_node2.τs)
    @test best_node4.τs == [0, 9]
    @test L(best_node4) > L(best_node2)

    # FNN with threshold
    Random.seed!(1234)
    optmodel4 = MCDTS.MCDTSOptimGoal(MCDTS.FNN_statistic(0.05), MCDTS.Continuity_function())
    tree = mcdts_embedding(Dataset(data[:,1]), optmodel4, w1, delays, runs)
    best_node = MCDTS.best_embedding(tree)
    @test best_node.τs == [0, 27, 21, 6]

    # FNN with threshold and less fid-points
    Random.seed!(1234)
    optmodel4 = MCDTS.MCDTSOptimGoal(MCDTS.FNN_statistic(0.05,2,0.8), MCDTS.Continuity_function())
    tree = mcdts_embedding(Dataset(data[:,1]), optmodel4, w1, delays, runs)
    best_node_less = MCDTS.best_embedding(tree)
    @test best_node_less.τs == [0, 27, 21, 6]

    L_YY = MCDTS.compute_delta_L(data[:,1], best_node.τs, delays[end];  w = w1)
    L_YY2 = MCDTS.compute_delta_L(data[:,1], best_node_less.τs, delays[end];  w = w1)
    @test L_YY == L_YY2

end

println("\nTesting MCDTS complete tree, Lorenz63 multivariate:")
@testset "MCDTS single rollout on multivariate data" begin

    Random.seed!(1234)
    tws = 2:4:delays[end]
    optmodel = MCDTS.MCDTSOptimGoal(MCDTS.L_statistic(-0.05,3,tws), MCDTS.Continuity_function())
    tree = MCDTS.mc_delay(data, optmodel, w, delays, runs2)
    best_node = MCDTS.best_embedding(tree)

    Random.seed!(1234)
    optmodel2 = MCDTS.MCDTSOptimGoal(MCDTS.FNN_statistic(0.05), MCDTS.Continuity_function())
    tree2 = MCDTS.mc_delay(data, optmodel2, w, delays, runs2)
    best_node2 = MCDTS.best_embedding(tree2)
    L_YY = MCDTS.compute_delta_L(data, best_node2.τs, best_node2.ts, delays[end];  w = w, tws = 2:4:delays[end])

    @test best_node.τs == [0, 9, 42, 20]
    @test best_node.ts == [1, 1, 1, 1]
    @test best_node2.τs == [0, 22, 16, 7]
    @test best_node2.ts == [1, 2, 2, 1]
    @test -0.92 < L(best_node) < -0.91
    @test -0.43 < L_YY < -0.42

    # less fid points
    Random.seed!(1234)
    optmodel2 = MCDTS.MCDTSOptimGoal(MCDTS.FNN_statistic(0.05,2,0.2), MCDTS.Continuity_function())
    tree3 = MCDTS.mc_delay(data, optmodel2, w, delays, runs2)
    best_node3 = MCDTS.best_embedding(tree3)

    @test best_node3.τs == [0, 48, 39, 0, 2]
    @test best_node3.ts == [1, 3, 2, 3, 2]

end
end
true