text:
# Packages used in this chapter
# To check for installed packages in V1.0 we now need to use the PKG API

using Pkg
ipks = Pkg.installed()
pkgs = ["BenchmarkTools","UnicodePlots","PyPlot"]

for p in pkgs
    !haskey(ipks,p) && Pkg.add(p)
end

using BenchmarkTools, UnicodePlots, PyPlot

# Define a function to compute the sum of the squares of (x,y)

sumsq(x,y) = x*x + y*y;

# Set the number of trials and initialise the counter
# Using zero(Integer) will ensure the counter is an integer

N = 10^8;
K = zero(Integer);

using Printf

# Sum the number of pairs which lie within the inscribed circle

for i = 1:N
  if sumsq(rand(), rand()) < 1.0
    global K += 1
  end
end

# ... and output the value for PI

@Printf.printf "Estimate of PI for %d trials is %8.5f\n" N 4.0*(K / N);

# We can time this using the @time macro
# Note that his takes some time, will be different in the REPL

@time for i = 1:N
  if sumsq(rand(), rand()) < 1.0
    global K += 1
  end
end

# Now need the global K to access top level varables in the REPL
# The error message is less than informative
# i.e. UndefVarError: K not defined

# Try v0.7 to get a warning and understand the problem
# ┌ Warning: Deprecated syntax `implicit assignment to global variable `K``.
# │ Use `global K` instead.
# └ @ none:0

# This ONLY applies to the top level loops
# So wrapping in let/end, begin/end block is OK, since K is in local scope

# Has is by using using a function and passing K,N as parameters
#
# Here is a contrived example which demonstrates this

while (true)
    N = 1000000
    K = zero(Integer)
    for i = 1:N
        if sumsq(rand(), rand()) < 1.0
          K += 1
        end
    end
    println(4.0*(K/N))
    break
end

# We can rewrite as a generator expressions
# This does not require the variable 'K'
#
# Comprehensions can also be written without the enclosing square brackets,
# producing an object known as a generator.

((rand()^2 + rand()^2) < 1.0 for i = 1:N)

# The generator can be "evaluated" as the argument to an aggregate function

count((rand()^2 + rand()^2) < 1.0 for i = 1:N) * 4 / N

# And this runs much more quickly with a much better memory allocation

@time count((rand()^2 + rand()^2) < 1.0 for i = 1:N) * 4 / N

# JULIA_HOME is not longer defined in Version 1
# See: https://docs.julialang.org/en/v1/manual/environment-variables/
#
# JULIA_BINDIR sets Sys.BINDIR

Sys.BINDIR

# Rather than just using the @time macro we can get a more accuracy timing using the BenchmarkTools package

using BenchmarkTools

# Define the function to sum the series
# It converges very slowly

function basel(N::Integer)
  @assert N > 0
  s = 0.0
  for i = 1:N
    s += 1.0/float(i)^2
  end
  return s
end

basel(10^8)   # Evaluate it over 10^8 terms

π*π/6.0    # Compare this to the analytic solution

# Use the BenchmarkTools macro @benchmark to get more unbiased timings

@benchmark basel(10^8) samples=10

using UnicodePlots

# Generate an array of the numbers from 1 to 100
# The ordinate value is create using a list comprehensive

x = collect(1:100);
y = [x[i]*sin(0.3*x[i])*exp(-0.03*x[i]) for i = 1:length(x)];

myPlot = lineplot(x, y, title = "My Plot", name = "chp-01")

# Alternatively this can be done using a map() construction

t = collect(0.0:0.1:10.0);
y =  map(x -> x*sin(3.0*x)*exp(-0.3*x), t);
myPlot = lineplot(t, y, title = "Second Plot", name = "chp-01")

using PyPlot
plt = PyPlot   # Will get function name clashes with UnicodePlots.

S0  = 100;      # Spot price
K   = 102;      # Strike price
r   = 0.05;     # Risk free rate
q   = 0.0;      # Dividend yield
v   = 0.2;      # Volatility
tma = 0.25;     # Time to maturity
T   = 90;       # Number of time steps

dt  = tma/T;
N = T + 1;
x = collect(0:T);

# Qualify the function calls since UnicodePlot has exports same functions into Main.

plt.title("Asian Option trajectories");
plt.xlabel("Time");
plt.ylabel("Stock price");

for k = 1:5
  S = zeros(Float64,N)
  S[1] = S0;
  dW = randn(N)*sqrt(dt);
  [ S[t] = S[t-1] * (1 + (r - q - 0.5*v*v)*dt + v*dW[t] + 0.5*v*v*dW[t]*dW[t]) for t=2:N ]
  plt.plot(x,S)
end

# Notice that the reference to the PyPlot object is in scope inside this loop
# The same is true when passing an array into a loop
# Confusing YES?
#
K = [0]; N=10^6;
for i = 1:N
  ((rand()^2 + rand()^2) .<  1.0) && (K[1] += 1)
end
mypi = 4*(K[1]/N)


text:
const DOC_ROOT_PATH   = "public"
const CONFIG_PATH     = "config"
const ENV_PATH        = joinpath(CONFIG_PATH, "env")
const APP_PATH        = "app"
const RESOURCES_PATH  = joinpath(APP_PATH, "resources")
const TEST_PATH       = "test"
const TEST_PATH_UNIT  = joinpath(TEST_PATH, "unit")
const LIB_PATH        = "lib"
const HELPERS_PATH    = joinpath(APP_PATH, "helpers")
const LOG_PATH        = "log"
const LAYOUTS_PATH    = joinpath(APP_PATH, "layouts")
const TASKS_PATH      = "tasks"
const BUILD_PATH      = "build"

const GENIE_CONTROLLER_FILE_POSTFIX     = "Controller.jl"
const GENIE_CHANNEL_FILE_POSTFIX        = "Channel.jl"

const ROUTES_FILE_NAME  = joinpath(CONFIG_PATH, "routes.jl")
const CHANNELS_FILE_NAME  = joinpath(CONFIG_PATH, "channels.jl")

const PARAMS_REQUEST_KEY    = :REQUEST
const PARAMS_RESPONSE_KEY   = :RESPONSE
const PARAMS_SESSION_KEY    = :SESSION
const PARAMS_FLASH_KEY      = :FLASH
const PARAMS_POST_KEY       = :POST
const PARAMS_GET_KEY        = :GET
const PARAMS_WS_CLIENT      = :WS_CLIENT
const PARAMS_JSON_PAYLOAD   = :JSON_PAYLOAD
const PARAMS_RAW_PAYLOAD    = :RAW_PAYLOAD
const PARAMS_FILES          = :FILES

const TEST_FILE_IDENTIFIER = "_test.jl"

# Used to store log info during app bootstrap, when the logger itself is not available.
# The queue is automatically emptied by the logger upon load.
const GENIE_LOG_QUEUE = Vector{Tuple{String,Symbol}}()

const VIEWS_FOLDER = "views"
const LAYOUTS_FOLDER = "layouts"

ext: jl


text:
using Documenter, SYNJULIA

makedocs(
    modules = [SYNJULIA],
    format = Documenter.HTML(; prettyurls = get(ENV, "CI", nothing) == "true"),
    authors = "sang",
    sitename = "SYNJULIA.jl",
    pages = Any["index.md"]
    # strict = true,
    # clean = true,
    # checkdocs = :exports,
)

deploydocs(
    repo = "github.com/SANG/SYNJULIA.jl.git",
    push_preview = true
)


text:
# The following function allow to permute the frames
export meanderingFramePermutation, fullFramePermutation

# TODO the following requires a test
function meanderingFramePermutation(f::MPIFile)
  N = tuple(calibSize(f)...)

  perm = Array{Int}(undef, N)
  for i in CartesianIndices(N)
    idx = [i[k] for k=1:length(i)]
    for d=2:3
      if isodd(sum(idx[d:3])-length(idx[d:3]))
        idx[d-1] = N[d-1] + 1 - idx[d-1]
      end
    end
    perm[i] = (LinearIndices(N))[idx...]
  end
  return vec(perm)
end

function fullFramePermutation(f::MPIFile, meandering::Bool)
  perm1=cat(measFGFrameIdx(f),measBGFrameIdx(f),dims=1)
  if meandering
    perm2=cat(meanderingFramePermutation(f),
               (length(perm1)-acqNumBGFrames(f)+1):length(perm1),dims=1)
    permJoint = perm1[perm2]
    return permJoint
  else
    return perm1
  end
end


text:
using Dates, FileIO, SeisIO, SeisNoise, PlotlyJS, Statistics, Printf
using IterTools
using PyCall
using PyPlot

@pyimport scipy.cluster.hierarchy as hi
@pyimport scipy.spatial.distance as di
@pyimport scipy.signal as sig

include("stacking.jl")

# set up plot
# layout = Layout(width=1200, height=800,
#                 xaxis_title="Lags [s]",
#                 yaxis_title="Days",
#                 font=attr(size=14),
#                 showlegend=false)
#
# p = PlotlyJS.plot([NaN], layout)
#
# distance_threshold = 0.8

"""
    clustered_selective_stacking(data::CorrData, reference::CorrData; threshold::Float64=0.0, slice::Union{Bool, Float64, Array{Float64,1}}=false, metric::String="cc", win_len::Float64=10.0, win_step::Float64=5.0)

Stack the windows in a CorrData object that exceed a correlation-coefficient threshold with respect to a reference.

# Arguments
- `data::CorrData,`    : Input CorrData matrix used to define the reference and each window
- `reference::CorrData`    : Reference cross-correlation function to use in computing the correlation coefficient
- `threshold::Float64,`    : Value of correlation-coefficient below which we zero windows
- `slice::Union{Bool, Float64, Array{Float64,1}}`    : whether to slice the cross-correlations before computing convergence. If so, a single Float64 specifies the
                                                       lag at which to slice to keep just the ballistic wave and an array of Float64 specifies the start and end lag
                                                       to keep only the coda.

# Output
- `stackedData::CorrData,`    : Slectively stacked data
- `cList::Array{Float64,1}`    : Correlation coefficient of each window with respect to the reference
"""
function clustered_selective_stacking(data::CorrData, reference::CorrData; threshold::Float64=0.0, slice::Union{Bool, Float64, Array{Float64,1}}=false, metric::String="cc", coh_win_len::Float64=10.0, coh_win_step::Float64=5.0, filter::Union{Bool, Array{Float64,1}}=false)
    # slice data if given a time window
    # TODO: find a better way to pass arguments to this function that only apply to coh, or only to cc. e.g., win_len has no meaning if metric=="cc"
    #if typeof(slice) != Bool
    ref = copy(reference)
    d_origin = copy(data)

    # time vector
    tvec = -reference.maxlag:1/reference.fs:reference.maxlag
    # find all times within [-slice,slice] time
    if slice ==false
        slice = 0.0
    end

    if typeof(slice) != Bool
        if typeof(slice)==Float64
            # find all times within [-slice,slice] time
            slice_inds = findall(x->(x<=-slice || slice <= x), tvec)
        elseif typeof(slice)==Array{Float64,1}
            # convert startlag/endlag[s] to startlag/windowlength[samples]
            slice_inds = findall(x->(-slice[2]<=x || x<=-slice[1] || slice[1]<=x || x<=slice[2]), tvec)
        else
            println("Please choose an allowable slicing operation. Exiting.")
            exit()
        end
    else
        # default to full reference cross-correlation
        tvec = -reference.maxlag:1/reference.fs:reference.maxlag
        slice_inds = findall(x->(x<=slice ||slice<= x), tvec)
    end

    #-------------------------------------------#
    # Work flow
    # 1. clustered by N groups ( N = 2-15 according to computation power)
    # 2. stacking them within the cluster
    # 3. grid search for best combination of stacking by taking correlation between reference
    # 4. select a set of cc function
    #-------------------------------------------#
    plot=false
    wiener = false
    # 1. clustering cc functions
    # distance is computed by pdist correlation coefficient

    if metric == "cc"

        # before doing clustering, remove incoherent shorttime window ccs by correlation coefficient like canonical selective stacking
        cList = get_cc(d_origin, ref)
        good_fit = findall(x->(x>=threshold), cList)
        d = copy(d_origin)
        #update d::CorrData only with high correlation coefficient
        d.corr = d_origin.corr[:, good_fit]

        pdis = di.pdist(transpose(d.corr[slice_inds,:]), metric="correlation") #'euclidean' 'seuclidean' 'matching' correlation
        #pdis = []
        #numofcc = size(d.corr)[2]
        #println(numofcc)
        # compute pdist array [d(2,1), d(3,1), ..., d(n, 1), d(3,2), ..., d(n, n-1)]
        # @simd for j = 1:numofcc-1
        #     print(j)
        #     @simd for i = j+1:numofcc
        #         # u1 = d.corr[:,i]
        #         # u2 = d.corr[:,j]
        #         # dt = 1/C.fs
        #         # b = 10
        #         # direction = 1
        #         # (_, _, _, error) = dtwdt(u1, u2, dt, maxLag=round(Int64, 0.05*C.maxlag*C.fs), b=b, direction=direction)
        #         dtemp = 1 - cor(d.corr[:,i], d.corr[:,j]) #https://www.researchgate.net/publication/265604603_Minimum_Pearson_Distance_Detection_for_Multilevel_Channels_With_Gain_andor_Offset_Mismatch
        #         push!(pdis, dtemp)
        #     end
        # end

    else
        error("Please specify metric to calculate distance array.")
    end

    Z = hi.linkage(pdis, method="ward") #method="single"

    #normalize distance
    Z[:,3] = Z[:,3]./maximum(Z[:,3])
    #println(Z)
    # this is used for dendrogram plotting
    maxclustnum = 12 # this is equal to the num of grid searching windows

    if plot
        PyPlot.figure(figsize=(12,8))
        fig = hi.dendrogram(Z, p=maxclustnum, truncate_mode="lastp")
        PyPlot.show()
        str = d.id
        PyPlot.savefig("./fig/dendrogram_cc_$str.png", dpi=300)
        PyPlot.close()
    end
    #clustering
    #group = hi.fcluster(Z, distance_threshold, criterion="distance")
    group = hi.fcluster(Z, maxclustnum, criterion="maxclust")

    #println(group)

    #stack each clusters

    #---plotting parameters ---#

    plot_xlag = length(d.corr[:,1]) # xlag of cc function
    plot_maxgroupnum = maxclustnum # max num of plotting group
    plot_maxynum = 15 # max num of traces in y axis
    xlabel = "clusters"
    ylabel = "clustered cc functions"

    if plot;
        layout = Layout(width=1400, height=800,
                xaxis_title=xlabel,
                yaxis_title=ylabel,
                font =attr(size=18),
                showlegend=false,
                title="Clustered cross-correlation function",
                paper_bgcolor="rgba(0,0,0,0)",
                plot_bgcolor="rgba(0,0,0,0)")
        p = PlotlyJS.plot([NaN], layout);
    end

    Stackedcc = copy(d)
    Stackedcc.corr = zeros(length(d.corr[:,1]), maximum(group))

    numtrace = 0

    #sort by correlation coefficient
    corrcoef = []
    for i = 1:maximum(group)
        idx = findall(x -> (x==i), group)
        tempData = copy(d)
        tempData.corr = tempData.corr[:, idx]
        stackedtemp = stack(tempData, allstack=true)
        Stackedcc.corr[:,i] = stackedtemp.corr[:,1]
        push!(corrcoef, cor(stackedtemp.corr[slice_inds,1], ref.corr[slice_inds,1]))
        #plot cc fnctions
    end

    sortidx = sortperm(corrcoef, rev=true)
    #print("corrcoef is: ")
    #println(corrcoef)

    #plot with sorted corr
    if plot
        traceall = GenericTrace{Dict{Symbol,Any}}[]
        for (i, id) = enumerate(sortidx)
            if id <= plot_maxgroupnum
                #idx = findall(x -> (x==id), group)
                idx = findall(x -> (x==i), group)
                tempData = copy(d)
                tempData.corr = tempData.corr[:, idx]
                stackedtemp = stack(tempData, allstack=true)

                # plot each cc function
                lags = (collect(1:plot_xlag) .+ (i-1)*(1.5 * plot_xlag)) ./ tempData.fs
                for j = 1:minimum([length(idx), plot_maxynum])
                    plot_ccfunc = 0.5 * tempData.corr[:, j] ./ maximum(tempData.corr[:, j]) .+ float(plot_maxynum - j)
                    trace1 = PlotlyJS.scatter(;x=lags, y=plot_ccfunc, mode="lines", line_color="rgb(0,0,0)")
                    #println(trace1)
                    #println(trace1.fields)
                    #PlotlyJS.addtraces!(p, trace1) # this may cause bug
                    push!(traceall, trace1)
                end
                # plot stacked trace
                plot_stackedccfunc = 0.5 * stackedtemp.corr[:, 1] ./ maximum(stackedtemp.corr[:, 1]) .+ float(plot_maxynum)
                #trace2 = PlotlyJS.scatter(;x=lags, y=plot_stackedccfunc, mode="lines+text",
                # line_color="rgb(255,0,0)", name="stacked trace", textposition="top center", text=[@sprintf("%4.3f",corrcoef[id])])
                trace2 = PlotlyJS.scatter(;x=lags, y=plot_stackedccfunc, mode="lines+text",
                 line_color="rgb(255,0,0)", name="stacked trace")
                push!(traceall, trace2)
                #PlotlyJS.addtraces!(p, trace2) # this may cause bug
            end
        end
        p = PlotlyJS.plot(traceall, layout)
        display(p)
        str = d.id
        PlotlyJS.savefig(p, "./fig/clustered_ccfunctions_$str.png")

    end

    if plot
        #plot reference
        lags = collect(1:plot_xlag)./ ref.fs
        plot_ref = 0.5 * ref.corr[:, 1] ./ maximum(ref.corr[:, 1])
        trace1 = PlotlyJS.scatter(;x=lags, y=plot_ref, mode="lines", line_color="rgb(0,0,255)", name="reference")
        p = PlotlyJS.plot(trace1, layout)
        display(p)
        str = d.id
        PlotlyJS.savefig(p, "./fig/refference_$str.png")
    end

    # performing grid search to find best combination

    # init coeff is of all stack with correlation coefficient > threshold
    stackedalltemp = stack(d, allstack=true)

    init_coef = cor(stackedalltemp.corr[slice_inds,1], ref.corr[slice_inds,1])
    maxcrosscoef = init_coef
    bestgroups = collect(1:maximum(group))
    #println(corrcoef)
    totalcalcnum = length(IterTools.subsets(1:maximum(group)))

    for (i, subs) = enumerate(IterTools.subsets(1:maximum(group)))
        if isempty(subs); continue; end

        if mod((totalcalcnum - i), 1000) == 0; println(totalcalcnum - i); end

        tempData = copy(Stackedcc)
        tempData.corr = tempData.corr[:, subs]
        stackedtemp = stack(tempData, allstack=true)
        # evaluate
        trial_crosscoef = cor(stackedtemp.corr[slice_inds,1], ref.corr[slice_inds,1])
        if trial_crosscoef > maxcrosscoef
            # crosscoef is increased, so update info
            #print("coef update: ")
            #println([maxcrosscoef, trial_crosscoef])
            bestgroups = subs
            maxcrosscoef = trial_crosscoef
        end
    end

    println([init_coef, maxcrosscoef])

    # prepare stacked data
    #println(bestgroups)
    stackedData = copy(d)
    stackedData.corr = zeros(size(d.corr)[1],1)
    for mergegroupid = bestgroups
        #println(mergegroupid)
        idx = findall(x -> (x==mergegroupid), group)
        tempData = copy(d)
        tempData.corr = tempData.corr[:, idx]
        stackedtemp = stack(tempData, allstack=true)
        stackedData.corr += stackedtemp.corr
    end

    if wiener
        # apply wiener filter before clustering
        wienerwindowsize = [3,3]
        d.corr = sig.wiener(stackedData.corr, wienerwindowsize)
    end

    return stackedData, init_coef, maxcrosscoef

end

max_stars_repo_name: kura-okubo/SeisXcorrelation

max_issues_repo_name: kura-okubo/SeisXcorrelation

max_forks_repo_name: kura-okubo/SeisXcorrelation


text:
module SegmentMSTs
using ..Types
using ..Evaluate

using Base.Threads

export segment, segment!, sgm2error, sgm2ec
#
# """
# equality
# """
# function Base.is(sgm1::SegMST, sgm2::SegMST)
#   sgm1.segmentation==sgm2.segmentation && sgm1.segmentPairs==sgm2.segmentPairs && sgm1.segmentPairAffinities==sgm2.segmentPairAffinities
# end

function Base.size(sgm::SegMST)
    return size( sgm.segmentation )
end

function Base.ndims(sgm::SegMST)
    ndims(sgm.segmentation)
end

"""
merge supervoxels with high affinity
"""
function Base.merge!(sgm::SegMST, thd::AbstractFloat)
    # the dict of parent and child
    # key->child, value->parent
    pd = Dict{UInt32,UInt32}()
    idxlst = Vector{Int64}()
    for idx in 1:length(sgm.segmentPairAffinities)
        if sgm.segmentPairAffinities[idx] > thd
            # the first one is child, the second one is parent
            pd[sgm.segmentPairs[idx,1]] = sgm.segmentPairs[idx,2]
        end
    end

    # find the root id
    for (c,p) in pd
        # list of child node, for path compression
        clst = [c]
        # find the root
        while haskey(pd, p)
            push!(clst, p)
            p = pd[p]
        end
        # now p is the root id
        # path compression
        for c in clst
            pd[c] = p
        end
    end

    # set each segment id as root id
    gc_enable(false)
    @threads for i in eachindex(sgm.segmentation)
        sgm.segmentation[i] = get(pd, sgm.segmentation[i], sgm.segmentation[i])
    end
    gc_enable(true)

    # update the segmentPairsrogram
    sgm.segmentPairAffinities = sgm.segmentPairAffinities[idxlst]
    sgm.segmentPairs = sgm.segmentPairs[:, idxlst]

    return sgm.segmentation
end

function Base.merge(sgm::SegMST, thd::AbstractFloat)
  sgm2 = deepcopy(sgm)
  return merge!(sgm2, thd)
end

"""
segment the sgm, only return segmentation
"""
function segment(sgm::SegMST, thd::AbstractFloat)
  sgm2 = merge(sgm, thd)
  return sgm2.segmentation
end

function segment!(sgm::SegMST, thd::AbstractFloat)
  sgm = merge!(sgm, thd)
  return sgm.segmentation
end

"""
compute segmentation error using one threshold
"""
function sgm2error(sgm::SegMST, lbl::Segmentation, thd::AbstractFloat)
    @assert thd<=1 && thd>=0
    seg = segment(sgm, thd)
    return evaluate(seg, lbl)
end

"""
compute error curve based on a segmentation (including segmentPairsrogram) and groundtruth label
"""
function sgm2ec(sgm::SegMST, lbl::Segmentation, thds = 0:0.1:1)
    ec = ScoreCurve()
    for thd in thds
        e = sgm2error(sgm, lbl, thd)
        e[:thd] = thd
        append!(ec, e)
    end
    ec
end

end # end of module


text:
#utilities
#some of these are candidates for Finometrics

@inline function finiteormissing!(df::DataFrame, s::Symbol)::Nothing
  df[!,s] .= (f::MFloat64->coalesce(isfinite(f),false) ? f : missing).(df[!,s])

  return nothing
end

@inline function finiteormissing!(df::AbstractDataFrame, s::Symbol)::Nothing
  @mpar for r ∈ eachrow(df)
    r[s] = coalesce(isfinite(r[s]),false) ? r[s] : missing
  end

  return nothing
end

#@inline Base.isfinite(::Missing)=missing

#NOTE: The below adjusts for logged vs unlogged returns
#const RETAINED_FACT_COLS = (( RETURN_FIELD == :lret) ?
#  [RETURN_FIELD; FACTOR_FIELDS; :lrfr, :lmkt] :
#  [RETURN_FIELD; FACTOR_FIELDS; :rfr, :mkt]

#consider promiting to the below to finometrics
@inline obliterate!(df::DataFrame)::DataFrame = select!(df, Not(names(df)))

#checks if a column or multiple columns are a key
@inline function checkdfkey(df::AbstractDataFrame, keys::Union{Symbol, Vector{Symbol}})
  N::Int = size(df,1)
  Nunique::Int = length(groupby(df, keys))

  (N==Nunique) || error("$keys is not a key!!")

  return nothing
end

#converts gzip to lz4
#=function gz2lz4(file::String; filetype::String = "csv", location::String = "",
  gzextension::String = "gz", lz4extension::String = "lz4",
  deleteold::Bool = false)

  infile::String = "$location\\$file.$filetype.$gzextension"
  outfile::String = "$location\\$file.$filetype.$lz4extension"

  isfile(infile) || error("file not found. Full path: $location\\$file.$filetype.$gzextension")

  #Below doesnt work
  indata = open(infile)
  outdata = open(outfile, "w")
  s = GzipDecompressorStream(LZ4HCCompressorStream(outdata))
  write(s, indata)
  close(s)
end=#

#converts gzip to lz4
function gz2zstd(file::String; filetype::String = "csv", location::Union{String,Nothing} = nothing,
  gzextension::String = "gz", zstdextension::String = "zstd",
  deleteold::Bool = false, level=1)

  local prefix::String

  prefix = isnothing(location) ?  "" : "$location\\"
  infile::String = "$prefix$file.$filetype.$gzextension"
  outfile::String = "$prefix$file.$filetype.$zstdextension"

  isfile(infile) || error("file not found. Full path: $prefix$file.$filetype.$gzextension")

  #Below doesnt work
  indata = open(infile)
  outdata = open(outfile, "w")
  s = GzipDecompressorStream(ZstdCompressorStream(outdata, level=level))

  try
    write(s, indata)
  finally
    close(s)
  end
end

#=function archivefile(file::String; inpath::String=DATA_PATH, inext::String = "",
  archivepath::String = ARCHIVE_PATH)

  (length(inext) > 0) && (file = "$file.$inext")
  inp::String = "$inpath\\$file"
  outp::String = "$archivepath\\$file.gz"

  open(GzipCompressorStream, outp, "w") do sout
    open(inp) do sin
      write(sout, sin)
    end
  end

  return nothing
end

function unarchivefile(file::String; outpath::String=DATA_PATH, outext::String = "",
  archivepath::String = ARCHIVE_PATH)

  (length(outext) > 0) && (file = "$file.$outext")
  outp::String = "$outpath\\$file"
  inp::String = "$archivepath\\$file.gz"

  open(GzipDecompressorStream, inp) do sin
    open(outp, "w") do sout
      write(sout, sin)
    end
  end

end=#

function recompressfiles(files = FILES_TO_RECOMPRESS; datapath=DATA_PATH)::Nothing

  for f ∈ files
    gz2zstd(f, location=datapath)
  end

  return nothing
end

recompressfile(file; datapath=DATA_PATH)::Nothing = recompressfiles(
  [file]; datapath=datapath)

function uniquebyyear(df::DataFrame; Frequired::Union{Nothing, Vector{Symbol}} = nothing,
  Fidentifier::Symbol = :permno, Fgroup::Symbol = :fyear)

  for sdf::SubDataFrame ∈ groupby(df, Fgroup)
    local ssdf::SubDataFrame

    if !isnothing(Frequired)
      ssdf = view(sdf, completecases(sdf[!, Frequired]), :)
    else
      ssdf = sdf
    end

    N::Int = length(unique(ssdf[!, Fidentifier]))
    N==0 && continue
    println("$(ssdf[1, Fgroup]): $N unique records by $Fidentifier ",
      isnothing(Frequired) ? "" : "containing $Frequired")
  end
end

#some useful functions for conforming numbers
@inline floororx(::Missing, ::Float64 = DENOM_FLOOR) = missing
@inline floororx(x::Real, floor::Float64 = DENOM_FLOOR) = max(x, floor)

@inline boundunit(::Missing, ::Float64=UNIT_TOL) = missing
@inline boundunit(x::Real, tol::Float64=UNIT_TOL) = max(min(x, 1.0 - tol), tol)
@inline boundunit1(x::Union{Real,Missing}) = boundunit(x, 0.0)

@inline boundabsunit(::Missing, ::Float64=UNIT_TOL) = missing
@inline boundabsunit(x::Real, tol::Float64=UNIT_TOL) = max(min(x, 1.0 -tol), -1.0 + tol)

@inline winsorizelevel(::Missing, ::Real, ::Real) = missing
@inline winsorizelevel(x::T, low::T, high::T) where {T<:Real} = min(max(low, x), high)

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
#!/bin/bash
#=
export JULIA_PROJECT=$(dirname $(dirname $(readlink -e "${BASH_SOURCE[0]}")))
exec julia --color=yes --startup-file=no "${BASH_SOURCE[0]}" "$@"
=#

if length(ARGS) < 2
    error("usage: $(PROGRAM_FILE) OUTFILE INFILE1 [INFILE2 [...]]")
else
    using UVH5Splicer
    outfile = popfirst!(ARGS)
    uvh5_splice(outfile, ARGS)
end


text:
module spaces

import Base.copy, Base.getindex, Base.setindex!, Base.keys
import Graphs
using SparseArrays
using LinearAlgebra

export Graph, DiGraph, getdomain

mutable struct Graph
    adjacency
    n
    updates_since_dropzeros

    function Graph(adjacency :: AbstractMatrix{Bool})
        m, n = size(adjacency)

        if m != n
            error("adjacency matrix must be square")
        end

        adjacency = adjacency .| adjacency'
        adjacency[diagind(adjacency)] .= 0
        new(adjacency, n, 0)
    end
end

function Graph(graph :: Graphs.Graph)
    A = Graphs.LinAlg.symmetrize(Graphs.LinAlg.adjacency_matrix(graph))
    A = Bool.(A)
    DiGraph(A)
end

Base.keys(g :: Graph) = [(i, j) for i in 1:g.n for j in 1:i-1]

function Base.getindex(g :: Graph, p)
  i, j = p
  g.adjacency[i, j]
end

getdomain(g :: Graph, i) = [false, true]

function Base.setindex!(g :: Graph, v, p)
    if typeof(g.adjacency) <: AbstractSparseMatrix
        # if sparse, don't insert zeros unnecessarily
        if v == 0 && g[p] == 0
            return
        end

        # if sparse, occasionally drop zeros
        if g.updates_since_dropzeros > 500
            dropzeros!(g.adjacency)
            g.updates_since_dropzeros = 0
        end

        g.updates_since_dropzeros += 1
    end

    i, j = p
    g.adjacency[i, j] = v
    g.adjacency[j, i] = v

end

function Base.copy(g :: Graph)
  Graph(copy(g.adjacency))
end

mutable struct DiGraph
    adjacency
    n
    updates_since_dropzeros

    function DiGraph(adjacency :: AbstractMatrix{Bool})
        m, n = size(adjacency)

        if m != n
            error("adjacency matrix must be square")
        end

        adjacency[diagind(adjacency)] .= 0
        new(adjacency, n, 0)
    end
end

function DiGraph(graph :: Graphs.DiGraph)
    A = Graphs.LinAlg.adjacency_matrix(graph)
    A = Bool.(A)
    DiGraph(A)
end

Base.keys(g :: DiGraph) = [(i, j) for i in 1:g.n for j in 1:g.n if i != j]

function Base.getindex(g :: DiGraph, p)
  i, j = p
  g.adjacency[i, j]
end

getdomain(g :: DiGraph, i) = [false, true]

function Base.setindex!(g :: DiGraph, v, p)
    if typeof(g.adjacency) <: AbstractSparseMatrix
        # if sparse, don't insert zeros unnecessarily
        if v == 0 && g[p] == 0
            return
        end

        # if sparse, occasionally drop zeros
        if g.updates_since_dropzeros > 500
            dropzeros!(g.adjacency)
            g.updates_since_dropzeros = 0
        end

        g.updates_since_dropzeros += 1
    end

    i, j = p
    g.adjacency[i, j] = v

end

function Base.copy(g :: DiGraph)
  DiGraph(copy(g.adjacency))
end

end


text:
using HnswAnn
using Test

ldpath = "/home/jpboth/Rust/hnswlib-rs/target/release/"
setRustlibPath(ldpath)

using Logging
using Base.CoreLogging

logger = ConsoleLogger(stdout, CoreLogging.Debug)
global_logger(logger)

include("distf32.jl")
@testset "distf32" begin
    @test testdistl1()
    @test testdistptr()
end

include("dumpload.jl")
@testset "dumpreload" begin
    @test testdump()
    @test testreload()
end


text:
using ITensors
using ITensorTDVP
using Observers
using Printf

using ITensorTDVP: tdvp_solver, process_sweeps, TDVPOrder

function tdvp_nonuniform_timesteps(
  solver,
  PH,
  psi::MPS;
  time_steps,
  reverse_step=true,
  time_start=0.0,
  order=2,
  (step_observer!)=Observer(),
  kwargs...,
)
  nsweeps = length(time_steps)
  maxdim, mindim, cutoff, noise = process_sweeps(; nsweeps, kwargs...)
  tdvp_order = TDVPOrder(order, Base.Forward)
  current_time = time_start
  for sw in 1:nsweeps
    sw_time = @elapsed begin
      psi, PH, info = tdvp(
        tdvp_order,
        solver,
        PH,
        time_steps[sw],
        psi;
        kwargs...,
        current_time,
        reverse_step,
        sweep=sw,
        maxdim=maxdim[sw],
        mindim=mindim[sw],
        cutoff=cutoff[sw],
        noise=noise[sw],
      )
    end
    current_time += time_steps[sw]

    update!(step_observer!; psi, sweep=sw, outputlevel, current_time)

    if outputlevel ≥ 1
      print("After sweep ", sw, ":")
      print(" maxlinkdim=", maxlinkdim(psi))
      @printf(" maxerr=%.2E", info.maxtruncerr)
      print(" current_time=", round(current_time; digits=3))
      print(" time=", round(sw_time; digits=3))
      println()
      flush(stdout)
    end
  end
  return psi
end

function tdvp_nonuniform_timesteps(H, psi::MPS; kwargs...)
  return tdvp_nonuniform_timesteps(tdvp_solver(; kwargs...), H, psi; kwargs...)
end


text:
"""
   set_size!(o, t)

Internal function to set a (width, height) tuple `t` to an object `o`.
"""
set_size!(o::Union{Figure,Colorbar}, t::T2F) = (o.size = t;)

"""
   set_texlabels!(o, v)

Internal function to set a bool variable indicating whether the figure has
LaTex labels or not. If there is no texpreamble defined, a simple one is
instantiated with `amssymb`.
"""
function set_texlabels!(
            o::Figure,
            v::Bool
         )::Nothing

   o.texlabels = ifelse(v, v, ∅)
   isnothing(o.texpreamble) || return
   set_texpreamble!(o, t"\usepackage{amssymb}")
   return
end

"""
   set_texscale!(o, v)

Sets a string variable indicating the type of scaling (`"fixed"`, `"scale"`
or `"none"`)
   * fixed: keep LaTeX in std size as close as possible to the ambient fsize
   * scale: force LaTeX to use a size corresponding to the ambient fsize
   * none:  do not scale LaTeX.
"""
function set_texscale!(o::Figure, v::String)
   v ∈ GLE_TEXSCALE || throw(
      OptionValueError("texscale", v)
   )
   o.texscale = v
   return
end

"""
   set_texpreamble!(o, v)

Sets the tex preamble that should be used (e.g. with font packages or symbol
packages). Bear in mind that only `pdflatex` can be used and so not all font
packages can be used nor some commands like `fontspec` which are meant for
XeLaTeX or LuaLaTex.
"""
set_texpreamble!(o::Figure, v::String) = (o.texpreamble = v;)

"""
   set_transparency!(o, v)

Sets a bool variable corresponding to whether to use cairo (supports
transparency) or not. This is for instance useful if some of the fill colors
used are transparent.
"""
set_transparency!(o::Figure, v::Bool) = (o.transparency = v;)


text:
"""
This is a sample input file for the CloudScat code.  To use it run julia
and type

julia> include("variable_droplet_radius.jl"); run()

The output file will be based on this file's name so to run different simulations
you can just copy this file, change the configuration and run it as above with
the new file name.  You can also run from the command line as e.g.

> julia variable_droplet_radius.jl

- - -
2020 Alejandro Luque IAA-CSIC (aluque@iaa.es)

"""

using CloudScat

const co = CloudScat.constants

# We specify a variable composition through a new type that allows CloudScat
# to dispatch on it whenever it has to access e.g. the radius or density.
struct Composition end

# For the droplet size we will set up a parabolic profile with a maximum at zmax
# where the radius is rmax.
const zmax, rmax = 9 * co.kilo, 20 * co.micro

# The minimum radius is reached at zmin and has a value of rmin.
const zmin, rmin = 12 * co.kilo, 10 * co.micro

# In this case we consider a constant density of scatterers.  One may use
# a more complex method in CloudScat.density below.
const nscat = 100 * co.centi^-3

# Now we overload the CloudScat.radius and CloudScat.density functions to
# dispacth on the type that we have just defined.
function CloudScat.radius(::Composition, r::Point)
    rmax - (rmax - rmin) * (r[3] - zmax)^2 / (zmin - zmax)^2
end
CloudScat.density(::Composition, r::Point) = nscat

function run()
    params = Params(
        # NUmber of simulated photons
        N = 400000,

        # The wavelength of the scattered light
        λ = 337 * co.nano,

        # The source is defined as a straight line running from source_a
        # to source_b.  Within this line photons are distributed uniformly with
        # initially random, isotropic directions.  To use a point source set
        # source_b = source_a
        source_a = [0, 0, 9 * co.kilo],
        source_b = [0, 0, 9 * co.kilo])

    # World geometry
    # Describe the cloud geometry here.  It may consist in combinations of
    # geometrical figures (cylinders, spheres, cones, semi-planes).
    # See cloud_gemoetry.jl for a more detailed description
    #
    # A cylinder is defined as Cylinder(zbottom, ztop, xcenter, ycenter, radius)
    c1 = Cylinder(6 * co.kilo,  12 * co.kilo, 0, 0, 200 * co.kilo)
    cloud = c1

    # Geometry of the full domain.  Photons outside the domain are
    # immediately discarded.
    domain = Cylinder(7 * co.kilo, 60 * co.kilo, 0, 0, 200 * co.kilo)

    # We also have to provide a max. collision rate for the null collision
    # method.  The only requirement is that the actual rate never exceeds this
    # but the code is more efficient if the bound is tight.  One option
    # is to provide the highest values of droplet radius and density but
    # if these peak at different places we incur some inefficiency.
    # It works for this case however.
    composition = VariableNR(params.λ, Composition(), rmax, nscat)

    # Define the full simulation world
    world = World(cloud, domain, composition)

    # Observers record a time curve and an image of received light.  You can
    # set up as many observers as you wish but they have a significant impact
    # on performance.
    observers = [Observer(
        # Location of the observer
        position = [0, 0, 400 * co.kilo],

        # Sampling interval
        tsample = 1e-5,

        # Collected samples
        nsamples = 10000,

        # fov of the camera (half of the diagonal fov).
        fov = 40,

        # camera pixels
        pixels = 1024)]

    basename = splitext(splitdir(@__FILE__)[2])[1]

    # The saveto option specifies the output file.  Leave it as it is to base
    # the filename on the name of this file specifying the parameters.
    CloudScat.main(params, world, observers,
                   saveto="$(basename).h5")
end

# Check if the code has been 'included' or run from the shell.  In the latter
# case, run the simulation.
isinteractive() || run()


text:
using Plots

# Window a time series and plot a feature
@shorthands featurepath
@recipe function f(::Type{Val{:featurepath}}, x, y, z; features=catch24, windows=slidingWindow(length(x)÷20), zscore=false)
    # x is time, y is time series
    X, t = windows(y)
    t = x[t]
    F = features(X)
    if zscore
        F = standardise(F)
    end

    seriestype := :path
    framestyle --> :box
    linewidth --> 2.5
    size --> (1200, 800)
    legend --> :outerright
    xguide --> "Time"
    yguide --> "Feature"
    for f ∈ getnames(features)
        @series begin
            label := string(f)
            x := t[1:end-1].+diff(t)./2
            y := (size(F, 1)==1 ? F[:] : vec(F[f, :]))
        end
    end
end


text:
export symmetric, hermitian, conjugate, diagonal, upper_triangular, lower_triangular
export unit_upper_triangular, unit_lower_triangular

#=
struct UniversalWrapper{T,S,F,C,D} <:AbstractMatrix{T}
    parent::S
end
indexflip(A::UniversalWrapper, i, j) = (i, j)
indexflip(A::UniversalWrapper{<:Any,<:Any,Union{Transpose,Adjoint}}, i, j) = (j, i)

mapper(::Type{<:Real},::Val) = identity
mapper(::Type{<:Complex},::Val{0}) = identity
mapper(::Type{<:Complex},::Val{1}) = conj
mapper(A::UniversalWrapper{T,<:Any,<:Any,C,:N}, i, j) where {T,C} = mapper(T, Val(C))
mapper(A::UniversalWrapper{T,<:Any,<:Any,C,:O}, i, j) where {T,C} = aij -> i == j ? one(T) : mapper(T, Val(C))(aij)
mapper(A::UniversalWrapper{T,<:Any,<:Any,C,:Z}, i, j) where {T,C} = aij -> i == j ? zero(T) : mapper(T, Val(C))(aij)
mapper(A::UniversalWrapper{T,<:Any,<:Any,C,:R}, i, j) where {T,C} = aij -> i == j ? real(aij) : mapper(T, Val(C))(aij)

Base.size(A::UniversalWrapper) = indexflip(A, size(A.parent))
Base.getindex(A::UniversalWrapper, i, j) = mapper(A, i, j)(getindex(A.parent, indexflip(A, i, j)...))
=#

"""
    combine wrappers of the following kinds to achieve

  1. try to reduce number of involved wrappers
  2. standardized sequence of combination
"""
WRAPPER_TYPES = [
                 Conjugate, Transpose, Adjoint, Symmetric, Hermitian,
                 UpperTriangular, LowerTriangular, UnitUpperTriangular, UnitLowerTriangular,
                 Diagonal, Bidiagonal, Tridiagonal, SymTridiagonal, HermiteTridiagonal]

import LinearAlgebra: Transpose, Adjoint, Symmetric, Hermitian
import LinearAlgebra: transpose, adjoint, symmetric, hermitian

# missing types: HermiteTridiagonal, Conjugate to be defined
conjugate(B::Conjugate) = B.parent
conjugate(B::Transpose) = adjoint(B.parent)
conjugate(B::Adjoint) = transpose(B.parent)
conjugate(B::Diagonal) = Diagonal(conj.(B.diag))
conjugate(B::AbstractMatrix{<:Real}) = B
conjugate(B::AbstractMatrix) = Conjugate(B)

transpose(B::Conjugate) = adjoint(B.parent)
## transpose(B::Transpose) = B.parent
transpose(B::Adjoint) = conjugate(B.parent)
# transpose(B::Symmetric) = B
# transpose(B::Hermitian) = conjugate(B)
# transpose(B::Diagonal) = B
## transpose(B::AbstractMatrix) = Transpose(B)

adjoint(B::Conjugate) = transpose(B.parent)
adjoint(B::Transpose) = conjugate(B.parent)
# adjoint(B::Adjoint) = B.parent
# adjoint(B::Symmetric) = conjugate(B)
## adjoint(B::Hermitian) = B
# adjoint(B::Diagonal) = conjugate(B)
## adjoint(B::AbstractMatrix) = Adjoint(B)

symmetric(B::Conjugate, uplo::Symbol=:U) = conjugate(symmetric(B.parent, uplo))
symmetric(B::Transpose, uplo::Symbol=:U) = symmetric(B.parent, toggle(uplo))
symmetric(B::Adjoint, uplo::Symbol=:U) = conjugate(symmetric(B.parent, toggle(uplo)))
## symmetric(B::Symmetric, uplo::Symbol=:U) = B
function symmetric(B::Hermitian, uplo::Symbol=:U)
    base = isreal(diag(B.data)) ? B.data : hermtotria(B)
    uplo == up(B) ? symmetric(base, uplo) : conjugate(symmetric(base, up(B)))
end
symmetric(B::Diagonal, uplo::Symbol=:U) = B
symmetric(B::UpperTriangular, uplo::Symbol=:U) = uplo == :U ? symmetric(B.data, :U) : Diagonal(diag(B.data))
symmetric(B::LowerTriangular, uplo::Symbol=:L) = uplo == :L ? symmetric(B.data, :L) : Diagonal(diag(B.data))
symmetric(B::UnitUpperTriangular, uplo::Symbol=:U) = uplo == :U ? symmetric(sparsecsc(B), :U) : Diagonal(oneunit.(diag(B.data)))
symmetric(B::UnitLowerTriangular, uplo::Symbol=:L) = uplo == :L ? symmetric(sparsecsc(B), :L) : Diagonal(oneunit.(real(diag(B.data))))
symmetric(B::Hermitian{<:Real}, uplo::Symbol=:U) = Symmetric(B.data, uplo)
symmetric(B::AbstractMatrix) = Symmetric(B, :U)

hermitian(B::Conjugate, uplo::Symbol=:U) = conjugate(hermitian(B.parent, uplo))
hermitian(B::Transpose, uplo::Symbol=:U) = conjugate(hermitian(B.parent, toggle(uplo)))
hermitian(B::Adjoint, uplo::Symbol=:U) = hermitian(B.parent, toggle(uplo))
hermitian(B::Symmetric, uplo::Symbol=:U) = uplo == up(B) ? hermitian(B.data, uplo) : conjugate(hermitian(B.data, up(B)))
hermitian(B::Hermitian, uplo::Symbol=:U) = B
hermitian(B::Diagonal, uplo::Symbol=:U) = real(B)
hermitian(B::UpperTriangular, uplo::Symbol=:U) = uplo == :U ? hermitian(B.data, :U) : Diagonal(real(diag(B.data)))
hermitian(B::LowerTriangular, uplo::Symbol=:L) = uplo == :L ? hermitian(B.data, :L) : Diagonal(real(diag(B.data)))
hermitian(B::UnitUpperTriangular, uplo::Symbol=:U) = uplo == :U ? hermitian(sparsecsc(B), :U) : Diagonal(oneunit.(diag(B.data)))
hermitian(B::UnitLowerTriangular, uplo::Symbol=:L) = uplo == :L ? hermitian(sparsecsc(B), :L) : Diagonal(oneunit.(real(diag(B.data))))
hermitian(B::AbstractMatrix{<:Real}, uplo::Symbol=:U) = symmetric(B, uplo)
hermitian(B::AbstractMatrix) = Hermitian(B, :U)

upper_triangular(B::Conjugate) = conjugate(upper_triangular(B.parent))
upper_triangular(B::Transpose) = transpose(lower_triangular(B.parent))
upper_triangular(B::Adjoint) = adjoint(lower_triangular(B.parent))
upper_triangular(B::Symmetric) = up(B) == :U ? upper_triangular(B.data) : transpose(lower_triangular(B.data))
function upper_triangular(B::Hermitian)
    base = isreal(diag(B.data)) ? B.data : hermtotria(B)
    up(B) == :U ? upper_triangular(base) : adjoint(lower_triangular(base))
end
upper_triangular(B::Diagonal) = B
upper_triangular(B::UpperTriangular) = B
upper_triangular(B::LowerTriangular) = Diagonal(diag(B))
upper_triangular(B::UnitUpperTriangular) = B
upper_triangular(B::UnitLowerTriangular) = Diagonal(oneunit.(diag(B.data)))
upper_triangular(B::AbstractMatrix) = UpperTriangular(B)

lower_triangular(B::Conjugate) = conjugate(lower_triangular(B.parent))
lower_triangular(B::Transpose) = transpose(upper_triangular(B.parent))
lower_triangular(B::Adjoint) = Adjoint(upper_triangular(B.parent))
lower_triangular(B::Symmetric) = up(B) == :L ? lower_triangular(B.data) : transpose(upper_triangular(B.data))
function lower_triangular(B::Hermitian)
    base = isreal(diag(B.data)) ? B.data : hermtotria(B)
    up(B) == :L ? lower_triangular(base) : adjoint(upper_triangular(base))
end
lower_triangular(B::Diagonal) = B
lower_Triangular(B::LowerTriangular) = B
lower_triangular(B::UpperTriangular) = Diagonal(diag(B.data))
lower_triangular(B::UnitLowerTriangular) = B
lower_triangular(B::UnitUpperTriangular) = Diagonal(oneunit.(diag(B.data)))
lower_triangular(B::AbstractMatrix) = LowerTriangular(B)

unit_upper_triangular(B::Conjugate) = Conjugate(unit_upper_triangular(B.parent))
unit_upper_triangular(B::Transpose) = Transpose(UnitLowerTriangular(B.parent))
unit_upper_triangular(B::Adjoint) = Adjoint(UnitLowerTriangular(B.parent))
unit_upper_triangular(B::Symmetric) = up(B) == :U ? unit_upper_triangular(B.data) : transpose(unit_lower_triangular(B.data))
function unit_upper_triangular(B::Hermitian)
    base = isreal(diag(B.data)) ? B.data : hermtotria(B)
    up(B) == :U ? unit_upper_triangular(base) : adjoint(unit_lower_triangular(base))
end
unit_upper_triangular(B::Diagonal) = Diagonal(oneunit.(diag(B)))
unit_upper_triangular(B::UnitUpperTriangular) = B
unit_upper_triangular(B::UnitLowerTriangular) = Diagonal(oneunit.(diag(B.data)))
unit_upper_triangular(B::UpperTriangular) = unit_upper_triangular(B.data)
unit_upper_triangular(B::LowerTriangular) = Diagonal(oneunit.(diag(B.data)))
unit_upper_triangular(B::AbstractMatrix) = UnitUpperTriangular(B)

unit_lower_triangular(B::Conjugate) = Conjugate(unit_lower_triangular(B.parent))
unit_lower_triangular(B::Transpose) = Transpose(unit_upper_triangular(B.parent))
unit_lower_triangular(B::Adjoint) = Adjoint(unit_upper_triangular(B.parent))
unit_lower_triangular(B::Symmetric) = up(B) == :L ? unit_lower_triangular(B.data) : Transpose(unit_upper_triangular(B.data))
function unit_lower_triangular(B::Hermitian)
    base = isreal(diag(B.data)) ? B.data : hermtotria(B)
    up(B) == :L ? unit_lower_triangular(base) : adjoint(unit_upper_triangular(base))
end
unit_lower_triangular(B::Diagonal) = Diagonal(oneunit.(diag(B)))
unit_lower_triangular(B::UnitLowerTriangular) = B
unit_lower_triangular(B::UnitUpperTriangular) = Diagonal(oneunit.(diag(B.data)))
unit_lower_triangular(B::LowerTriangular) = unit_lower_triangular(B.data)
unit_lower_triangular(B::UpperTriangular) = Diagonal(oneunit.(diag(B.data)))
unit_lower_triangular(B::AbstractMatrix) = UnitLowerTriangular(B)

diagonal(B::Conjugate) = conjugate(Diagonal(diag(B.parent)))
diagonal(B::Transpose) = Diagonal(diag(B.parent))
diagonal(B::Adjoint) = conjugate(Diagonal(diag(B.parent)))
diagonal(B::Symmetric) = Diagonal(diag(B.data))
diagonal(B::Hermitian) = Diagonal(real(diag(B.data)))
diagonal(B::Diagonal) = B
diagonal(B::UpperTriangular) = Diagonal(diag(B.data))
diagonal(B::LowerTriangular) = Diagonal(diag(B.data))
diagonal(B::UnitUpperTriangular) = Diagonal(oneunit.(diag(B.data)))
diagonal(B::UnitLowerTriangular) = Diagonal(oneunit.(diag(B.data)))
diagonal(B::AbstractMatrix) = Diagonal(diag(B))

function hermtotria(B::Hermitian)
    A = sparsecsc((up(B) == :U ? UpperTriangular : LowerTriangular)(B.data))
    for i = 1:size(A, 1)
        A[i,i] = real(A[i,i])
    end
    A
end

#
# verification of type depth of combinations of wrappers.
# does the combiniation of wrappers define a finite group?
#
wrappertype(A) = begin wr = _wrt(A); isempty(wr) ? [identity] : wr end
_wr(A::T) where T = T.name.wrapper

_wrt(A::Union{Conjugate,Transpose,Adjoint}) = Any[ _wr(A); _wrt(A.parent)]
_wrt(A::Union{Symmetric,Hermitian}) = Any[ (_wr(A), up(A)); _wrt(A.data)]
_wrt(A::Union{UpperTriangular,LowerTriangular}) = Any[ _wr(A); _wrt(A.data)]
_wrt(A::Union{UnitUpperTriangular,UnitLowerTriangular}) = Any[ _wr(A); _wrt(A.data)]
_wrt(A::Union{Diagonal,Tridiagonal,SymTridiagonal}) = Any[ _wr(A) ]
_wrt(A::Bidiagonal) = Any[ (_wr(A), up(A)) ]
_wrt(A::Any) = []

apply(T) = A -> T(A)
apply(t::Tuple) = A -> t[1](A, t[2])
apply(a::Vector) = A -> (x = apply(a[end])(A); length(a) == 1 ? x : apply(a[1:end-1])(x))

WRAPPERS_BASE = [identity, Transpose, Adjoint, Conjugate,
                (Symmetric, :U), (Symmetric, :L), (Hermitian, :U), (Hermitian, :L),
                UpperTriangular, LowerTriangular,
                UnitUpperTriangular, UnitLowerTriangular, Diagonal
               ]

WRAPPERS_COMPLETE = [WRAPPERS_BASE...,
                     Any[Transpose, UpperTriangular], Any[Transpose, LowerTriangular],
                     Any[Transpose, UnitUpperTriangular], Any[Transpose, UnitLowerTriangular],
                     Any[Adjoint, UpperTriangular], Any[Adjoint, LowerTriangular],
                     Any[Adjoint, UnitUpperTriangular], Any[Adjoint, UnitLowerTriangular],
                     Any[Conjugate, UpperTriangular], Any[Conjugate, LowerTriangular],
                     Any[Conjugate, UnitUpperTriangular], Any[Conjugate, UnitLowerTriangular],
                     [Conjugate, (Symmetric, :U)], [Conjugate, (Symmetric, :L)],
                     [Conjugate, (Hermitian, :U)], [Conjugate, (Hermitian, :L)],
                     Any[Conjugate, Diagonal]
                    ]


text:
#-----------------------------------------------------------------# partitionplot
@userplot struct PartitionPlot
    args
end

@recipe function f(o::PartitionPlot; nparts=100, stat=Extrema(), by=nothing, dropmissing=false)
    t = o.args[1]
    sel = map(x -> lowerselection(t, x), o.args[2:end])
    # Create the (Indexed)Partition object
    o = if length(sel) == 1
        Partition(stat, nparts)
    elseif length(sel) == 2
        T = typeof(collect(rows(t)[1])[sel[1]])  # Type of first selection
        IndexedPartition(T, stat, nparts)
    else
        throw(ArgumentError("too many arguments for partitionplot"))
    end
    # Wrap the (Indexed)Partition in an FTSeries to possibly remove missing values
    s = FTSeries(o; filter = dropmissing ? !IndexedTables._ismissing : x -> true)
    if by === nothing
        reduce(s, t; select=sel)
    else
        grp = groupreduce(s, t, by; select=sel)
        grp = isa(t, DDataset) ? collect(grp) : grp
        for row in rows(grp)
            @series begin row[2] end
        end
    end
end

"""
    partitionplot(table, y;    stat=Extrema(), nparts=100, by=nothing, dropmissing=false)
    partitionplot(table, x, y; stat=Extrema(), nparts=100, by=nothing, dropmissing=false)

Plot a summary of variable `y` against `x` (`1:length(y)` if not specified).  Using `nparts`
approximately-equal sections along the x-axis, the data in `y` over each section is
summarized by `stat`.
"""
partitionplot


text:
export FDSNevq, FDSNevt

"""
    (H,R) = FDSNevq(ot)

Multi-server query for the events with the closest origin time to `ot`.
Returns an Array{SeisHdr,1} in H with event headers and an Array{SeisSrc,1}
in R in H with corresponding source process info.

### Keywords
| KW       | Default      | T [^1]    | Meaning                        |
|----------|:-----------  |:----------|:-------------------------------|
| evw      | [600., 600.] | Float64   | search window in seconds [^2]  |
| mag      | [6.0, 9.9]   | Float64   | search magitude range          |
| nev      | 0            | Integer   | events per query [^3]          |
| rad      | []           | Float64   | radius search                  |
| reg      | []           | Float64   | geographic search region       |
| src [^4] | "IRIS"       | String    | data source; `?seis_www` lists |
| to       | 30           | Int64     | timeout (s) for web requests   |
| v        | 0            | Integer   | verbosity                      |

[^1]: `Array{T, 1}` for `evw`, `mag`, `rad`, `reg`; `T` for others
[^2]: search range is always `ot-|evw[1]| ≤ t ≤ ot+|evw[2]|`
[^3]: if `nev=0`, all matches are returned.
[^4]: In an event query, keyword `src` can be a comma-delineated list, like `"IRIS, INGV, NCEDC"`.

See also: `SeisIO.KW`, `?seis_www`
"""
function FDSNevq(ot::String;
  evw::Array{Float64,1} = [600.0, 600.0],
  mag::Array{Float64,1} = [6.0, 9.9],
  nev::Integer = 0,
  rad::Array{Float64,1} = KW.rad,
  reg::Array{Float64,1} = KW.reg,
  src::String = KW.src,
  to::Int = KW.to,
  v::Integer = KW.v)

  if isempty(reg) && !isempty(rad)
    if length(rad) == 4
      append!(rad, [-30.0, 700.0])
    end
    search_coords = string( "&latitude=", rad[1], "&longitude=", rad[2],
                            "&minradius=", rad[3], "&maxradius=", rad[4],
                            "&mindepth=", rad[5], "&maxdepth=", rad[6] )
  else
    if isempty(reg)
      reg = Float64[-90.0, 90.0, -180.0, 180.0, -30.0, 700.0]
    elseif length(reg) == 4
      append!(reg, [-30.0, 700.0])
    end
    search_coords = string("&minlat=", reg[1], "&maxlat=", reg[2],
                           "&minlon=", reg[3], "&maxlon=", reg[4],
                           "&mindepth=", reg[5], "&maxdepth=", reg[6])
  end

  # Determine time window
  ot2::Float64 = try
    d2u(DateTime(ot))
  catch
    if length(ot) <= 14
      ot0 = string(ot[1:4],"-",ot[5:6],"-",ot[7:8],"T",ot[9:10],":",ot[11:12])
      if length(ot) > 12
        ot1 = string(ot0, ":", ot[13:14])
      else
        ot1 = string(ot0, ":00")
      end
    end
    d2u(DateTime(ot1))
  end
  d0 = string(u2d(ot2 - abs(evw[1])))
  d1 = string(u2d(ot2 + abs(evw[2])))
  oti = round(Int64, ot2*sμ)

  # multi-server query (most FDSN servers do NOT have an event service)
  sources = String.(strip.(split(lowercase(src) == "all" ? "EMSC, INGV, IRIS, LMU, NCEDC, NIEP, ORFEUS, SCEDC, USGS, USP" : src, ",")))
  sources = [strip(i) for i in sources]
  EvCat = Array{SeisHdr,1}(undef, 0)
  EvSrc = Array{SeisSrc,1}(undef, 0)
  origin_times = Array{Int64,1}(undef, 0)
  for k in sources
    v > 1 && println(stdout, "Querying ", k)
    url = string(fdsn_uhead(String(k)), "event/1/query?",
                                        "starttime=", d0, "&endtime=", d1,
                                        search_coords,
                                        "&minmag=", mag[1], "&maxmag=", mag[2],
                                        "&format=xml")
    v > 0 && println(stdout, "URL = ", url)
    req_info_str = "\nFDSN event query:"

    (R, parsable) = get_http_req(url, req_info_str, to)
    if parsable
      str_req = String(R)
      v > 1 && println(stdout, "REQUEST BODY:\n", str_req)
      xdoc = parse_string(str_req)
      event_xml!(EvCat, EvSrc, xdoc)
      v > 1 && println(stdout, "CATALOG:\n", EvCat)
    end
  end

  if nev > 0
    # Sort based on earliest origin time
    # sort!(EvCat, by = H -> abs(round(Int64, d2u(getfield(H, :ot))*sμ) - oti))
    for H in EvCat
      # push!(origin_times, round(Int64, d2u(getfield(H, :ot))*sμ))
      push!(origin_times, abs(round(Int64, d2u(getfield(H, :ot))*sμ) - oti))
    end
    # k = sortperm(abs.(origin_times.-oti))
    k = sortperm(origin_times)

    n0 = min(length(EvCat), nev)
    n0 < nev && @warn(string("Catalog only contains ", n0, " events (original request was ", nev,")"))
    return EvCat[k[1:n0]], EvSrc[k[1:n0]]
    # return EvCat[1:n0], EvSrc[1:n0]
  else
    return EvCat, EvSrc
  end
end

"""
    FDSNevt(ot::String, chans::String)

Get header and trace data for the event closest to origin time `ot` on channels
`chans`. Returns a SeisEvent structure.

### Keywords
| KW       | Default      | T [^1]    | Meaning                        |
|----------|:-----------  |:----------|:-------------------------------|
| evw      | [600., 600.] | Float64   | search window in seconds [^2]  |
| fmt      | "miniseed"   | String    | request data format            |
| len      | 120.0        | Float64   | desired trace length [s]       |
| mag      | [6.0, 9.9]   | Float64   | search magitude range          |
| model    | "iasp91"     | String    | velocity model for phases      |
| nd       | 1            | Real      | number of days per subrequest  |
| opts     | ""           | String    | user-specified options[^3]     |
| pha      | "P"          | String    | phases to get  [^4]            |
| rad      | []           | Float64   | radius search                  |
| reg      | []           | Float64   | geographic search region       |
| src      | "IRIS"       | String    | data source; `?seis_www` lists |
| to       | 30           | Int64     | timeout (s) for web requests   |
| v        | 0            | Integer   | verbosity                      |
| w        | false        | Bool      | write requests to disk?        |

[^1]: KW is `Array{T, 1}` for `evw`, `mag`, `rad`, `reg`, type `T` for others
[^2]: Search range is always `ot-|evw[1]| ≤ t ≤ ot+|evw[2]|`
[^3]: Format like an http request string, e.g. "szsrecs=true&repo=realtime" for FDSN. String shouldn't begin with an ampersand.
[^4]: Comma-separated String, like `"P, pP"`; use `"ttall"` for all phases

See also: `distaz!`, `FDSNevq`, `FDSNsta`
"""
function FDSNevt(ot::String, chans::ChanOpts;
  evw::Array{Float64,1} = [600.0, 600.0],
  fmt::String           = KW.fmt,
  len::Real             = 120.0,
  mag::Array{Float64,1} = [6.0, 9.9],
  model::String         = "iasp91",
  nd::Real              = KW.nd,
  opts::String          = KW.opts,
  pha::String           = "P",
  rad::Array{Float64,1} = KW.rad,
  reg::Array{Float64,1} = KW.reg,
  src::String           = KW.src,
  to::Int64             = KW.to,
  v::Integer            = KW.v,
  w::Bool               = KW.w
  )

  C = fdsn_chp(chans, v)

  # Create header
  v > 0 && println(stdout, now(), ": event query begins.")
  (H,R) = FDSNevq(ot, nev=1,
                  rad=rad,
                  reg=reg,
                  mag=mag,
                  src=src,
                  to=to,
                  evw=evw,
                  v=v
                  )
  H = H[1]
  R = R[1]

  # Create channel data
  v > 0 && println(stdout, now(), ": data query begins.")
  s = H.ot                                      # Start time for FDSNsta is event origin time
  t = u2d(d2u(s) + 60*len)                      # End time is len minutes later
  (d0, d1) = parsetimewin(s,t)
  S = SeisData()
  # FDSNget!(S, C, fmt=fmt, nd=nd, rad=rad, reg=reg, s=d0, si=true, src=src, t=d1, to=to, v=v, w=w)
  FDSNget!(S, C, d0, d1, false, fmt, false, nd, "", rad, reg, true, src, to, v, w, "FDSNsta.xml", false)

  v > 1 && println(stdout, now(), ": channels initialized.")
  v > 2 && println(stdout, S)

  # Initialize SeisEvent structure
  Ev = SeisEvent(hdr = H, source = R, data = S)
  v > 1 && println(stdout, now(), ": SeisEvent created.")
  v > 1 && println(stdout, S)

  # Update Ev with distance, azimuth
  distaz!(Ev)
  v > 1 && println(stdout, now(), ": Δ,Θ updated.")

  # Add phase arrivals to :data
  v > 0 && println(stdout, now(), ": phase query begins.")
  get_pha!(Ev, pha=pha, model=model, to=to, v=v)
  return Ev
end


text:
export GA

function selection(fits::Array{Float64}, n::Int64=3)
    # return the index of the winner of a n-way tournament
    fshuffle = randperm(length(fits))[1:n]
    winner = indmax(fits[fshuffle])
    fshuffle[winner]
end

function GA(nin::Int64, nout::Int64, fitness::Function;
            ctype::DataType=CGPChromo, seed::Int64=0, expert::Any=nothing,
            id::String="")
    population = Array{ctype}(Config.ga_population)
    fits = -Inf*ones(Float64, Config.ga_population)
    population = [ctype(nin, nout) for i in 1:Config.ga_population]
    best = population[1]
    max_fit = -Inf
    eval_count = 0

    nelite = Int64(round(Config.ga_population * Config.ga_elitism_rate))
    ncross = Int64(round(Config.ga_population * Config.ga_crossover_rate))
    nmutate = Int64(round(Config.ga_population * Config.ga_mutation_rate))
    if (nelite + ncross + nmutate) > Config.ga_population
        rates = [Config.ga_elitism_rate, Config.ga_crossover_rate,
                 Config.ga_mutation_rate]
        norm_rates = rates ./ sum(rates)
        nelite = Int64(floor(Config.ga_population * norm_rates[1]))
        ncross = Int64(floor(Config.ga_population * norm_rates[2]))
        nmutate = Config.ga_population - (nelite + ncross)
        ncopy = 0
    else
        ncopy = Config.ga_population - (nelite + ncross + nmutate)
    end

    @assert (nelite + ncross + nmutate + ncopy) == Config.ga_population

    while eval_count < Config.total_evals
        # evaluation
        log_gen = false
        for p in eachindex(population)
            if fits[p] == -Inf
                fit = fitness(population[p])
                fits[p] = fit
                eval_count += 1
                if fit > max_fit
                    max_fit = fit
                    best = clone(population[p])
                    log_gen = true
                end
                if eval_count == Config.total_evals
                    log_gen = true
                    break
                end
            end
        end

        eval(Config.log_function)(id, seed, eval_count, max_fit, best, GA,
                                  ctype, log_gen)

        if eval_count == Config.total_evals
            break
        end

        # create new population
        new_pop = Array{ctype}(Config.ga_population)
        new_fits = -Inf*ones(Float64, Config.ga_population)
        popi = 1

        # elites
        sinds = collect(eachindex(fits))
        sort!(sinds, by=x->fits[x], rev=true)
        for i in 1:nelite
            new_pop[popi] = clone(population[sinds[i]])
            new_fits[popi] = fits[sinds[i]]
            popi += 1
        end

        # crossover
        for i in 1:ncross
            p1 = population[selection(fits)]
            p2 = population[selection(fits)]
            new_pop[popi] = crossover(p1, p2)
            popi += 1
        end

        # mutation
        for i in 1:nmutate
            parent = population[selection(fits)]
            new_pop[popi] = mutate(parent)
            popi += 1
        end

        # copy
        for i in 1:ncopy
            new_pop[popi] = clone(population[selection(fits)])
            new_fits[popi] = fits[popi]
            popi += 1
        end

        population = new_pop
        fits = new_fits

        # size limit
        for i in eachindex(population)
            if length(population[i].nodes) > Config.node_size_cap
                population[i] = ctype(nin, nout)
                fits[i] = -Inf
            end
        end
    end

    max_fit, best.genes
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 036555ee-7758-11eb-11e0-194ebe3d6f15
html"<button onclick='present()'>present</button>"

# ╔═╡ 1fc8fccc-7758-11eb-3d10-b1c72df65ebd
md"""
# Programming scientifically in Julia

- Julia is easy to read and write (but so are Python or Matlab)
- Julia is performant and scalable* (but so are C or Fortran)
- Julia is interactive (but so are Python or Mathematica)

> Researchers often find themselves coding algorithms in one programming language, only to have to rewrite them in a faster one.

**Julia** solves the two-language problem.
"""

# ╔═╡ cf2303d2-7b45-11eb-3301-89fab00d5c9d
md"""
## Power of (a good) Language

Abelman and Sussman, Structure and Interpretation of Computer Programs
> Programs must be written for people to read, and only incidentally for machines to execute.

Bertrand Russell
> Language serves not only to express thoughts, but to make possible thoughts which could not exist without it.

Ludwig Wittgenstein
> The limits of my language mean the limits of my world.
"""

# ╔═╡ 1340a278-7b47-11eb-2958-cff4ccdd630a
md"""
## So why exactly is Julia great?

#### Julia is easy to learn

This easiness arrises from a clear and concise syntax (similar but quite more refined than Python), which results in small but clear code.
"""

# ╔═╡ be149b64-7b47-11eb-14cc-2ba364843de5
md"""
## So why exactly is Julia great?

#### Julia is (easily) fast
Tricky definition:
- Most languages are "equivalent" on a theoretical level but in practise greatly differ
- One can write terribly slow programs in languages that are supposed to be fast

The norm is that something in pure Julia will run a lot faster than in most dynamic languages.

**Native Julia can beat highly optimised libraries** (some of the fastest BLAS-like operations are in `LoopVectorization.jl`). No magic: Julia shares the same LLVM backend with e.g. C++ so the compiled performances can be very similar.

Don't focus on the fact that Julia is faster or not but **that it's simply much simpler to write high-performant and clear code in Julia than in other languages**.
"""

# ╔═╡ 528ae93c-7b49-11eb-12f2-63a759823973
md"""
## So why exactly is Julia great?

#### Julia is (easily) fast

Calculate the mean of an array ignoring NaNs [in Julia](https://github.com/mlubin/NaNMath.jl/blob/master/src/NaNMath.jl#L216)

```
function mean_count(x::AbstractArray{T}) where T<:AbstractFloat
    z = zero(eltype(x))
    sum = z
    count = 0
    @simd for i in x
        count += ifelse(isnan(i), 0, 1)
        sum += ifelse(isnan(i), z, i)
    end
    result = sum / count
    return (result, count)
end
```
"""

# ╔═╡ 8b1823f0-7b49-11eb-227f-ab9723138570
md"""
## So why exactly is Julia great?

#### Julia is (easily) fast

Calculate the mean of an array ignoring NaNs [in Python (Numpy)]
(https://github.com/pydata/bottleneck/blob/b62427287f61ad9c12b8101b5c94a6b7b90c9a14/bottleneck/src/reduce_template.c#L213-L314)

```
/* dtype = [['float64'], ['float32']] */
REDUCE_ALL(nanmean, DTYPE0) {
    Py_ssize_t count = 0;
    npy_DTYPE0 ai, asum = 0;
    INIT_ALL
    BN_BEGIN_ALLOW_THREADS
    WHILE {
        FOR {
            ai = AI(DTYPE0);
            if (!bn_isnan(ai)) {
                asum += ai;
                count += 1;
            }
        }
        NEXT
    }
    BN_END_ALLOW_THREADS
    if (count > 0) {
        return PyFloat_FromDouble(asum / count);
    } else {
        return PyFloat_FromDouble(BN_NAN);
    }
}

REDUCE_ONE(nanmean, DTYPE0) {
    Py_ssize_t count;
    npy_DTYPE0 ai, asum;
    INIT_ONE(DTYPE0, DTYPE0)
    BN_BEGIN_ALLOW_THREADS
    if (LENGTH == 0) {
        FILL_Y(BN_NAN)
    } else {
        WHILE {
            count = 0;
            asum = 0;
            FOR {
                ai = AI(DTYPE0);
                if (!bn_isnan(ai)) {
                    asum += ai;
                    count += 1;
                }
            }
            if (count > 0) {
                asum /= count;
            } else {
                asum = BN_NAN;
            }
            YPP = asum;
            NEXT
        }
    }
    BN_END_ALLOW_THREADS
    return y;
}
/* dtype end */

/* dtype = [['int64', 'float64'], ['int32', 'float64']] */
REDUCE_ALL(nanmean, DTYPE0) {
    Py_ssize_t total_length = 0;
    npy_DTYPE1 asum = 0;
    INIT_ALL
    BN_BEGIN_ALLOW_THREADS
    WHILE {
        FOR asum += AI(DTYPE0);
        total_length += LENGTH;
        NEXT
    }
    BN_END_ALLOW_THREADS
    if (total_length > 0) {
        return PyFloat_FromDouble(asum / total_length);
    } else {
        return PyFloat_FromDouble(BN_NAN);
    }
}

REDUCE_ONE(nanmean, DTYPE0) {
    npy_DTYPE1 asum;
    INIT_ONE(DTYPE1, DTYPE1)
    BN_BEGIN_ALLOW_THREADS
    if (LENGTH == 0) {
        FILL_Y(BN_NAN)
    } else {
        WHILE {
            asum = 0;
            FOR asum += AI(DTYPE0);
            if (LENGTH > 0) {
                asum /= LENGTH;
            } else {
                asum = BN_NAN;
            }
            YPP = asum;
            NEXT
        }
    }
    BN_END_ALLOW_THREADS
    return y;
}
/* dtype end */
```
"""

# ╔═╡ f9c5c958-7b47-11eb-05c5-59c9455b6ed1
md"""
## So why exactly is Julia great?

#### Julia is (uniquely) both strongly-typed and dynamic

- _static, compiled, user types_ (in C, Fortran, etc)
- _dynamic, interpreted, standard types_ (in Python, Mathematica, etc)
- **dynamic, compiled, user types** (Julia)

In dynamic languages the types of variables don't have to be known an runtime (type annotation is optional).

In compiled languages the functions are compiled (properly optimised) and don't suffer from interpretation overhead at runtime.

- Julia's type system has limits, no arrow types (f: Rⁿ → Rᵐ)
- The types has less dynamism than in Python (cannot add fields to a type)
"""
# . In Julia the number and type of fields are static such that their byte size and memory setups can be determined. So a concretely typed Vector{Float64} can be very fast (while in Python a List is essentially a Vector{Any} and one needs to resort to Numpy for speed).

# ╔═╡ 21ec2982-7b49-11eb-2fc3-a799a8b09712
md"""
## So why exactly is Julia great?

#### Julia has multiple dispatch

> People come to Julia for speed but stay for the multiple dispatch

Difference from OOP (single-dispatch)
```
class Car
    def __init__(self, color, brand):
		# instance variable unique to each instance
		self.color = color
		self.brand = brand

	def information(self):
		return "I'm a {} {}".format(self.color, self.brand)

car = Car("yellow", "mercedes")
car.information()
```

Philosophy: Methods shouldn’t belong to a specific data type:
> It makes little sense to artificially deem the operations to "belong" to one argument more than any of the others.

Since there are no classes in Julia, the structures only contain data but no methods
```julia
struct Car
	color
	brand
end

information(c::Car) = "I'm a $(c.color) $(c.brand)"

car = Car("yellow", "mercedes")
information(car)
```

##### Levels of dispatch
- none: `f(x1,x2,x3,...)` constant expressive power (e.g., Python functions)
- single: `(x1::T1).f(x2,x3,...)` linear expressive power (e.g., Python class methods)
- multiple `f(x1::T1,x2::T2,x3,...)` exponential expressive power (Julia's core  paradigm)

#### Sharing types
Much simpler to define **new** operations on existing types

- Extending existing operations
```julia
information(n::Number) = "I'm a number"
```
Give several functions the same name, because they perform conceptually similar tasks, but operate on different types.

- Defining new operations
```julia
is_beautiful(c::Car) = (c.color == "rose gold")
```

How to do it in OOP? With inheritance one needs to use another name (`NewCar`) or edit the original class or simply drop (single) dispatch.
"""

# ╔═╡ 2cf41e6e-7b49-11eb-02bf-23174a7e8ce3
md"""
## So why exactly is Julia great?

#### Julia is garbage-collected
No need to worry about memory management.

#### Julia supports concurrent, parallel and distributed computing
Natively and with simple and clear syntax.

#### Julia interfaces particularly well with other languages
Can directly call Python or C.
"""

# ╔═╡ Cell order:
# ╟─036555ee-7758-11eb-11e0-194ebe3d6f15
# ╟─1fc8fccc-7758-11eb-3d10-b1c72df65ebd
# ╟─cf2303d2-7b45-11eb-3301-89fab00d5c9d
# ╟─1340a278-7b47-11eb-2958-cff4ccdd630a
# ╟─be149b64-7b47-11eb-14cc-2ba364843de5
# ╟─528ae93c-7b49-11eb-12f2-63a759823973
# ╟─8b1823f0-7b49-11eb-227f-ab9723138570
# ╟─f9c5c958-7b47-11eb-05c5-59c9455b6ed1
# ╟─21ec2982-7b49-11eb-2fc3-a799a8b09712
# ╟─2cf41e6e-7b49-11eb-02bf-23174a7e8ce3

max_stars_repo_name: fmeirinhos/JuliaSeminar

max_issues_repo_name: Fromeworld/JuliaSeminar

max_forks_repo_name: Fromeworld/JuliaSeminar


text:
mutable struct Callback <: AbstractCallback
  every::Int64
  func::Function
  function Callback(every::Int64,func::Function)
    new(every,func)
  end
end
Callback(func::Function) = Callback(1,func)

function callback(cb::AbstractCallback,system::AbstractSystem,model::AbstractPotential,args...)
  cb.func(system,model,args)
end


text:
using ArgParse
using Formatting
using JSON,TOML
using DataStructures

include("runCmd.jl")

s = ArgParseSettings()
@add_arg_table s begin

    "--casename"
        help = "Prefix of the file."
        arg_type = String
        required = true

    "--input-dir"
        help = "The folder where case folders are contained"
        arg_type = String
        required = true

    "--output-dir"
        help = "Output folder"
        arg_type = String
        required = true

    "--year-rng"
        help = "Range of years being diagnosed"
        arg_type = Int64
        nargs = 2
        required = true

    "--domain-file"
        help = "Domain file of atmosphere"
        arg_type = String
        required = true

    "--overwrite"
        action = :store_true

end

parsed = DataStructures.OrderedDict(parse_args(ARGS, s))

JSON.print(parsed, 4)

for y = parsed["year-rng"][1]:parsed["year-rng"][2], m=1:12

    date_str = format("{:04d}-{:02d}", y, m)

    old_file = joinpath(parsed["input-dir"], "$(parsed["casename"]).cam.h0.$(date_str).nc")

    new_file1 = joinpath(parsed["output-dir"], "$(parsed["casename"]).cam_extra1.$(date_str).nc")
    if !isfile(new_file1) || parsed["overwrite"]
        println("Generating file: $(new_file1)")
        pleaseRun(`ncap2 -h -O -v -s 'PREC_TOTAL=PRECC+PRECL;' $old_file $new_file1`)
    end

    new_file2 = joinpath(parsed["output-dir"], "$(parsed["casename"]).cam_extra2_zmean.$(date_str).nc")
    if !isfile(new_file2) || parsed["overwrite"]
        println("Generating file: $(new_file2)")
        pleaseRun(`ncwa -O -a lon -v T,U,V,Z3,ilev $old_file $new_file2`)
    end

end

println("Diagnose streamfunction")
pleaseRun(`julia $(@__DIR__)/lib/streamfunction.jl
    --input-data-file-prefix $(parsed["output-dir"])/$(parsed["casename"]).cam_extra2_zmean.
    --output-data-file-prefix $(parsed["output-dir"])/$(parsed["casename"]).cam_extra3_streamfunction.
    --domain-file $(parsed["domain-file"])
    --beg-year $(parsed["year-rng"][1])
    --end-year $(parsed["year-rng"][2])
    --V-varname V
`)

#=
println("Diagnose ENSO index")
pleaseRun(`julia $(@__DIR__)/lib/ENSO.jl
    --input-data-file-prefix $(parsed["output-dir"])/$(parsed["casename"]).cam.h0.
    --output-data-file-prefix $(parsed["output-dir"])/$(parsed["casename"]).cam_extra4_ENSO.
    --domain-file $(parsed["domain-file"])
    --beg-year $(parsed["year-rng"][1])
    --end-year $(parsed["year-rng"][2])
    --SST-varname SST
`)

=#

max_stars_repo_name: meteorologytoday/CESM-diagnostic

max_issues_repo_name: meteorologytoday/CESM-diagnostic

max_forks_repo_name: meteorologytoday/CESM-diagnostic


text:
using Test
using Documenter
using Interp1d
doctest(Interp1d)


text:
# This file is a part of DataDeps.jl. License is MIT.

"""
    progress_update_period()

Returns the period between updated being logged on the progress.
This is used by the default `fetch_method` and is generally a good idea
to use it in any custom fetch method, if possible
"""
function progress_update_period()
    envvar = get(ENV, "DATADEPS_PROGRESS_UPDATE_PERIOD") do
        if haskey(ENV, "DATADEPS_ALWAYS_ACCEPT")
            # Running in a script, probably want minimal updates
            "Inf"
        else # default
            "5" # seconds
        end
    end
    parse(Float32, envvar)
end

"""
    fetch_default(remote_path, local_path)

The default fetch method.
It tries to be a little bit smart to work with things other than just HTTP.
See also [`fetch_base`](@ref) and [`fetch_http`](@ref).
"""
function fetch_default(remotepath, localdir)
    if remotepath isa AbstractString && occursin(r"^https?://", remotepath)
        # It is HTTP, use good HTTP method, that gets filename by HTTP rules
        return fetch_http(remotepath, localdir)
    else
        # More generic fallback, hopefully `Base.basename`
        return fetch_base(remotepath, localdir)
    end
end

"""
  fetch_base(remote_path, local_dir)

Download from `remote_path` to `local_dir`, via `Base` mechanisms.
The download is performed using `Base.download`
and `Base.basename(remote_path)` is used to determine the filename.
This is very limitted in the case of HTTP as the filename is not always encoded in the URL.
But it does work for simple paths like `"http://myserver/files/data.csv"`.
In general for those cases prefer `http_download`.

The more important feature is that this works for anything that has overloaded
`Base.basename` and `Base.download`, e.g. [`AWSS3.S3Path`](https://github.com/JuliaCloud/AWSS3.jl).
While this doesn't work for all transport mechanisms (so some datadeps will still a custom `fetch_method`),
it works for many.
"""
function fetch_base(remote_path, local_dir)
    localpath = joinpath(local_dir, basename(remote_path))
    return Base.download(remote_path, localpath)
    return string(localpath)
end

"""
    fetch_http(remotepath, localdir; update_period=5)

Pass in a HTTP[/S] URL  and a directory to save it to,
and it downloads that file, returing the local path.
This is using the HTTP protocol's method of defining filenames in headers,
if that information is present.
Returns the localpath that it was downloaded to.

`update_period` controls how often to print the download progress to the log.
It is expressed in seconds. It is printed at `@info` level in the log.
By default it is once per second, though this depends on configuration

"""
function fetch_http(remotepath, localdir; update_period=progress_update_period())
    @assert(isdir(localdir))
    return HTTP.download(remotepath, localdir; update_period=update_period)
end


text:
# Code to match https://www.youtube.com/watch?v=omz_NdFgWyU&list=PLQVvvaa0QuDcjD5BAw2DxE6OF2tius3V3&index=6

# Usings.
using Random
Random.seed!(0)

##############################
# Stuff from previous videos #
##############################

# Define the spiral_data function to get some data.
function spiral_data(points::Int, classes::Int)

    X = zeros(Float64, points * classes, 2)
    y = zeros(Int64, points * classes)

    for class_number = 1:classes

        ix = points*(class_number-1)+1:points*class_number
        r = range(0, 1, length=points)
        t = range((class_number-1)*4, class_number*4, length=points) .+ randn(points)*0.2

        @. X[ix, :] = [r*sin(t*2.5) r*cos(t*2.5)]
        @. y[ix] = class_number
    end

    return X, y
end

# Struct for a layer.
mutable struct Layer_Dense
    weights::Array{Float64, 2}
    biases::Vector{Float64}
end

# Make a constructor for this struct.
function Layer_Dense(n_inputs::Int, n_neurons::Int)
    weights = 0.10 * randn(n_inputs, n_neurons)
    biases = zeros(n_neurons)

    return Layer_Dense(weights, biases)
end

# And define the forward function for this struct.
function forward(layer::Layer_Dense, inputs)
    return inputs * layer.weights .+ layer.biases'
end

# The rectified linear unit.
function activation_relu(inputs)
    # Using the dot to take the element-wise maximum.
    return max.(0., inputs)
end

#############
# New stuff #
#############

# The softmax activation function.
function activation_softmax(x)

    # Protect against overflow.
    x .-= maximum(x, dims=2)

    exp_values = exp.(x)
    return exp_values ./ sum(exp_values, dims=2)
end

###############
# The example #
###############
X, y = spiral_data(100, 3)

dense1 = Layer_Dense(2, 3)
dense2 = Layer_Dense(3, 3)

output1 = forward(dense1, X)
output1_a = activation_relu(output1)

output2 = forward(dense2, output1_a)
output2_a = activation_softmax(output2)

println(output2_a[1:5, :])


text:
using Dates
using Test

using MarketData

using TimeSeries

@testset "modify" begin

@testset "update method works" begin
    new_cls  = update(cl, today(), 111.11)
    new_clv  = update(cl, today(), [111.11])
    new_ohlc = update(ohlc, today(), [111.11 222.22 333.33 444.44])
    empty1   = TimeArray(Vector{Date}(), Array{Int}(undef, 0,1))
    empty2   = TimeArray(Vector{Date}(), Array{Int}(undef, 0,2))

    @testset "empty time arrays can be constructed" begin
        @test length(empty1) == 0
        @test length(empty2) == 0
    end

    @testset "update an empty time array fails" begin
        @test_throws ArgumentError update(empty1, Date(2000,1,1), 10)
        @test_throws ArgumentError update(empty2, Date(2000,1,1), [10 11])
    end

    @testset "update a single column time array with single value vector" begin
        @test last(values(new_clv)) == 111.11
        @test_throws DimensionMismatch update(cl, today(), [111.11, 222.22])
    end

    @testset "update a multi column time array" begin
        @test values(tail(new_ohlc, 1)) == [111.11 222.22 333.33 444.44]
        @test_throws MethodError update(ohlc, today(), [111.11, 222.22, 333.33])
    end

    @testset "cannot update more than one observation at a time" begin
        @test_throws(
            MethodError,
            update(cl, [Date(2002,1,1), Date(2002,1,2)], [111.11, 222,22]))
    end

    @testset "cannot update oldest observations" begin
        @test_throws ArgumentError update(cl, Date(1999,1,1), [111.11])
        @test_throws ArgumentError update(cl, Date(1999,1,1), 111.11)
    end

    @testset "cannot update in-between observations" begin
        @test_throws ArgumentError update(cl, Date(2000,1,8), [111.11])
        @test_throws ArgumentError update(cl, Date(2000,1,8), 111.11)
    end
end

@testset "rename method works" begin
    re_ohlc = rename(ohlc, [:a, :b, :c, :d])
    re_cl   = rename(cl, [:vector])
    re_cls  = rename(cl, :symbol)

    @testset "change colnames with multi-member vector" begin
        @test colnames(re_ohlc) == [:a, :b, :c, :d]
        @test_throws DimensionMismatch rename(ohlc, [:a])
    end

    @testset "change colnames with single-member vector" begin
        @test colnames(re_cl) == [:vector]
        @test_throws DimensionMismatch rename(cl, [:a, :b])
    end

    @testset "change colnames with symbol" begin
        @test colnames(re_cls) == [:symbol]
        @test_throws MethodError rename(cl, :symbol_a, :symbol_b)
    end
end

end  # @testset "modify"


text:
using Distributions
using TensorFlow

# Generate some synthetic data
x = randn(100, 50)
w = randn(50, 10)
y_prob = exp(x * w)
y_prob ./= sum(y_prob, 2)

function draw(probs)
    y = zeros(size(probs))
    for i in 1:size(probs, 1)
        idx = rand(Categorical(probs[i, :]))
        y[i, idx] = 1
    end
    return y
end

y = draw(y_prob)

# Build the model
sess = Session(Graph())
X = placeholder(Float64)
Y_obs = placeholder(Float64)

variable_scope("logistic_model", initializer = Normal(0, 0.001)) do
    global W = get_variable("weights", [50, 10], Float64)
    global B = get_variable("bias", [10], Float64)
end

Y=nn.softmax(X * W + B)
Loss = -reduce_sum(log(Y) .* Y_obs)
optimizer = train.AdamOptimizer()
minimize_op = train.minimize(optimizer, Loss)
saver = train.Saver()
# Run training
run(sess, global_variables_initializer())
checkpoint_path = mktempdir()
info("Checkpoint files saved in $checkpoint_path")
for epoch in 1:10
    cur_loss, _ = run(sess, [Loss, minimize_op], Dict(X=>x, Y_obs=>y))
    println(@sprintf("Current loss is %.2f.", cur_loss))
    train.save(saver, sess, joinpath(checkpoint_path, "logistic"), global_step=epoch)
end


text:
include("../src/hypergraph-affinity-functions.jl")
using StatsBase
name = "primary"

dataset = "contact-$name-school-classes-gender"
r = 25
H, labels= read_hypergraph_data(dataset,r)

m = [2 3 6 7 9 11 14 15 17]
w = [1 4 5 8 10 12 13 16 18]

keep = findall(x->x<19,labels)
H = H[:,keep]
labels = labels[keep]
n = size(H,2)
classes = zeros(Int64,n)
women = findall(x->in(x,w),labels)
classes[women] .= 1

## Bootstrap
alpha, B1, B2, H1, H2, TR1, TR2, N = Hypergraph_to_Scores(H,classes,5)
for k = 2:4
B = 100
T1 = TR1[k,1:k]  # True affinity scores
T0 = TR2[k,1:k]

TA1 = H1[k,1:k]  # True affinity scores
TA0 = H2[k,1:k]

# Bootstraps
R0, R1, MR1, MR0, SR1, SR0, A0, A1, MA1, MA0, SA1, SA0 = Bootstrap_Affinities(H,k,B,classes)

E1 = vec(maximum(R1,dims = 1))
E1m = vec(minimum(R1,dims = 1))
E0 = vec(maximum(R0,dims = 1))
E0m = vec(minimum(R0,dims = 1))

# Plot Ratio Scores
lw = 2.5
ms = 8
gfs = 15
tfs = 15
titlesize = 14
s1 = 350
s2 = 300

sc = 1
color1 = :green
color2 = :blue
title = ""
leg = false
xlab = "Affinity type t"
ylab = "Affinity / Baseline "
p = plot()
p = Plots.plot([1; k],[1; 1],fontfamily = "helvetica",linewidth = 2,yscale = :identity,legend = false,grid = false, size = (s1,s2),color = :gray)
title!(p,title,titlefontsize = titlesize)
xaxis!(p,xlab, xticks = 1:k,tickfontsize=tfs,guidefontsize = gfs)
yaxis!(p,ylab,tickfontsize=tfs,guidefontsize = gfs)

# [M0-2*SE0 M0+2*SE]
# [E0m E0]
# [E1m E1]
p = plot!(1:k,[T0 T0], fillrange=[MR0-sc*SR0 MR0+sc*SR0], fillalpha=0.3, color = color2, linewidth = lw, markersize = ms,
    markershape = :circle, markerstrokecolor = color2, label = class2)
p = plot!(1:k,[T1 T1], fillrange=[MR1-sc*SR1 MR1+sc*SR1],fillalpha=0.3,color = color1, linewidth = lw, markersize = ms,
    markershape = :circle,markerstrokecolor = color1, label = class1)
savefig("Figures/Ratios_$(k)_$(name)_Bootstrap.pdf")

# Plot just affinity scores
title = ""
leg = false
xlab = "Affinity type t"
ylab = "Affinity"

p = Plots.plot(fontfamily = "helvetica",linewidth = 2,yscale = :identity,legend = false,grid = false, size = (s1,s2),color = :gray)
title!(p,title,titlefontsize = titlesize)
xaxis!(p,xlab, xticks = 1:k,tickfontsize=tfs,guidefontsize = gfs)
yaxis!(p,ylab,tickfontsize=tfs,guidefontsize = gfs)

# [M0-2*SE0 M0+2*SE]
# [E0m E0]
# [E1m E1]
p = plot!(1:k,[TA0 TA0], fillrange=[MA0-sc*SA0 MA0+sc*SA0], fillalpha=0.3, color = color2, linewidth = lw, markersize = ms,
    markershape = :circle, markerstrokecolor = color2, label = class2)
p = plot!(1:k,[TA1 TA1], fillrange=[MA1-sc*SA1 MA1+sc*SA1],fillalpha=0.3,color = color1, linewidth = lw, markersize = ms,
    markershape = :circle,markerstrokecolor = color1, label = class1)
    savefig("Figures/Affinities_$(k)_$(name)_Bootstrap.pdf")
end

max_stars_repo_name: nveldt/HypergraphHomophily

max_issues_repo_name: nveldt/HypergraphHomophily

max_forks_repo_name: nveldt/HypergraphHomophily


text:
"""
    PadFirstAxis(d, N1)

Zero-padding first axis of N-dimensional array `d`.

# Arguments
* `d::Array{Real, N}`: N-dimensional data.
* `N1::Int`: total number of samples for the first axis of the zero-padded data.

# Examples
```julia
julia> w = Ricker(); wpad = PadFirstAxis(w, 128); plot(w); plot(wpad+1.0)

julia> d, ext = SeisHypEvents(); dpad = PadFirstAxis(d, 512);
SeisPlot(d); SeisPlot(dpad)
```
"""

function PadFirstAxis{T<:Real, N}(a::Array{T,N}, N1::Int)
    dims = size(a)
    n1 = dims[1]
    N1 > n1 || error("N1 must be greater than size(a,1)")
    b = zeros(T,N1,dims[2:end]...)
    b[1:dims[1],:] = a
    return b
end

max_stars_repo_name: fercarozzi/myseismicjulia

max_issues_repo_name: fercarozzi/myseismicjulia

max_forks_repo_name: fercarozzi/myseismicjulia


text:
__precompile__(true)
module FixedSizeArrays

using Compat

if isdefined(Base, :normalize)
    import Base: normalize
else
    import Compat: normalize
    export normalize
end

importall Base
import Base.LinAlg.chol!

if VERSION >= v"0.5.0-dev+4677"
    import Base.LinAlg._chol!
end

# for 0.5 and 0.4 compat, use our own functor type
abstract Functor{N}

if VERSION < v"0.5.0-dev+698"
    macro pure(ex)
        esc(ex)
    end
else
    import Base: @pure
end

include("core.jl")
include("functors.jl")
include("constructors.jl")

# put them here due to #JuliaLang/julia#12814
# needs to be before indexing and ops, but after constructors
immutable Mat{Row, Column, T} <: FixedMatrix{Row, Column, T}
    values::NTuple{Column, NTuple{Row, T}}
end

# most common FSA types
immutable Vec{N, T} <: FixedVector{N, T}
    values::NTuple{N, T}
end
immutable Point{N, T} <: FixedVector{N, T}
    values::NTuple{N, T}
end

include("mapreduce.jl")
include("destructure.jl")
include("indexing.jl")
include("ops.jl")
include("expm.jl")
include("array_of_fixedsize.jl")
include("conversion.jl")

function show{R,C,T}(io::IO, m::Mat{R,C,T})
	println(io, typeof(m), "(")
	for i=1:R
		println(io, "    ", join(row(m, i), " "))
	end
	println(io, ")")
end

show(io::IO, v::FixedVector{0}) = print(io, typeof(v).name.name, "()")
function show{N,T}(io::IO, v::FixedVector{N,T})
    print(io, typeof(v).name.name, "(", v[1])
    for i = 2:N
        print(io, ",", v[i])
    end
    print(io, ")")
end

export FixedArray
export FixedVector
export FixedMatrix
export MutableFixedArray
export MutableFixedVector
export MutableFixedMatrix
export Mat, Vec, Point
export @fsa
export similar_type
export construct_similar

export unit

export row
export column
export MatMulFunctor
export setindex
export eltype_or, size_or, ndims_or
export @fslice
export destructure
export push, pop, shift, unshift, deleteat, insert

end


text:
using RDatasets, DataArrays, Gadfly

set_default_plot_size(10inch, 3inch)

plot(dataset("car", "Vocab"),
     xgroup="Year", x="Vocabulary",
     Geom.subplot_grid(Geom.histogram))


text:
### Function making an empty cell
function Empty_Cell!(v::Vehicle)
    v.speed = Int8(-1)
    v.tipo = Int8(-1)
    v.change = Int8(0)
    v.num = 0
    v.len = 0
end

## Function making an obstacle cell
function Obstacle!(v::Vehicle)
    v.speed = Int8(-1)
    v.tipo = Int8(-2)
    v.change = Int8(0)
    v.num = Int8(0)
    v.len = 0
end

## Function changing the values of a vehicle
function Change_Vehicle!(v, speed, pos, tipo, change, num, length)
    v.speed = Int8(speed)
    v.position = Int64(pos)
    v.tipo = Int8(tipo)
    v.change = Int8(change)
    v.num = num
    v.len = length
end

## Function inserting a vehicle with probability p of being a truck in the first six cells of the highway
function Insert_Vehicle!(highway, num::Int64, p::Float64, vmax::Array{Int8, 1} = Int8[3, 5], left_boundary::Int8 = Int8(5), Len::Array{Int8, 1} = Int8[2, 1])

    lim = left_boundary + Int8(1)

    while highway[lim].tipo == -1 && highway[lim+Int8(1)].tipo != 1 && lim < (left_boundary+vmax[2])
        lim += Int8(1)
    end

    v_new = ( rand() <= p ? vmax[1] : vmax[2] )
    length_new = (v_new == vmax[1] ? Len[1] : Len[2])
    tipo = (v_new == vmax[1] ? 1 : 2)

    pos_new = min( lim - v_new, v_new)

    if highway[pos_new + 1].tipo != 1
      if tipo == 2
        Change_Vehicle!(highway[pos_new], v_new, pos_new, tipo, 1, num, Int64(length_new))
      end

      if tipo == 1
        if highway[pos_new-1] == -1
          Change_Vehicle!(highway[pos_new], v_new, pos_new, tipo, 1, num, Int64(length_new))
        else # If there is no place, find one
          pos_new = Int8(6)
          found = false
          while pos_new > 1 && !found
            if highway[pos_new].tipo == -1 && highway[pos_new-1].tipo == -1 && highway[pos_new+1].tipo != 1
              found = true
            else
              pos_new -= Int8(1)
            end
          end
          if found
            Change_Vehicle!(highway[pos_new], v_new, pos_new, tipo, 1, num, Int64(length_new))
          end
          if !found && pos_new == Int8(1)
            if highway[pos_new].tipo == -1 && highway[pos_new+1].tipo != 1
              Change_Vehicle!(highway[pos_new], v_new, pos_new, tipo, 1, num, Int64(length_new))
            end
          end
        end
      end
    end
    lim = v_new = pos_new = tipo = left_boundary = vmax = length_new = 0
end

## Function giving the position of the vehicle ahead when used in an empty cell
## When used at a non empty cell, it returns the position of the cell
## The function Vehicle_Ahead counts obstacles
function Vehicle_Ahead( position, Highway, N)
    i = position
    while i < N && Highway[i].tipo == -1
        i += 1
    end
    return i
end

## Function giving the position of the vehicle behind, even when we call it at a non empty cell
## If there's no vehicle behing, it returns the last cell
## It does not count obstacles
function Vehicle_Behind( position, Highway, N, left_boundary::Int8 = Int8(5) )

    i = position - Highway[position].len
    while i > left_boundary && (Highway[i].tipo == -1 || Highway[i].tipo == -2)
        i -= 1
    end

    if i > left_boundary
        return i
    else
        return N
    end
end

## Function giving the position of the vehicle closest to the end
function Pos_left(highway)

    gap = highway[end].position

    while gap >= 1 && (highway[gap].tipo == -1 || highway[gap].tipo == -2)
        gap -= 1
    end
    gap
end

## Same as Pos_left but without counting the first 9 cells
function Pos_left9(highway)

    gap = highway[end].position

    while gap >= 9 && (highway[gap].tipo == -1 || highway[gap].tipo == -2)
        gap -= 1
    end
    gap
end

## Function giving the position of the vehicle closest to the begining of the highway
function Pos_right(highway, tipo::Int8 = Int8(0))
    x = highway[1].position
    i = 1

    if tipo == 0
      while x < highway[end].position && (highway[i].tipo == -1 ||highway[i].tipo == -2)
          x += 1
          i += 1
      end
    else
      while x < highway[end].position && highway[i].tipo != tipo
          x += 1
          i += 1
      end
    end
    x
    #if highway[x].tipo == 1
    #  return x + highway[x].len - 1
    #elseif highway[x].tipo == 2
    #  return x + highway[x].len - 1
    #end
end

## Function giving the last empty cell in a highway
function Pos_next(highway)
    x = highway[end].position
    i = length(highway)
    while x > highway[1].position && (highway[i].tipo == 1 || highway[i].tipo == 2)
        x -= 1
        i -= 1
    end
    x
end


text:
@testset "research_log.jl" begin

    data_surv = CSO.dataset("cso_raw")
    data_trudeau = CSO.dataset("Montreal-Trudeau_15min_precipitations")
    data = innerjoin(data_surv, data_trudeau, on=:Date)
    dropmissing!(data, :Duration)
    addColOverflow!(data, "Duration", "Overflow", 5)
    overflowOnly!(data, "Overflow")

    topredict = "Duration"

    train, test = train_test_year(data, 2020)

    result = CR_gammaReg_log(topredict, [:d15min, :d45min, :d1h, :d6h], train, test)
    varnames = result[1, :Variables]
    RMSD = result[1, :RMSD]
    variables = stepRMSD_gammaReg_log(topredict, [:d45min, :d15min, :d1h, :d6h], [], train, test, false)

    #The RMSD is well derived
    @test RMSD == gammaReg_log("Duration", varnames, train, test)

    #The df is well sorted
    @test RMSD == minimum(result.RMSD)

    #The higher RMSD of the CR_gammaReg_log corresponds to the higher found stepwise
    @test RMSD == gammaReg_log("Duration", variables, train, test)

end

max_stars_repo_name: AlexLsn/CSO

max_issues_repo_name: AlexLsn/CSO

max_forks_repo_name: AlexLsn/CSO


text:
using Turing
using Distributions
using Base.Test

num_samples = Int(1e6)

resSystematic = Turing.resampleSystematic( [0.3, 0.4, 0.3], num_samples )
resStratified = Turing.resampleStratified( [0.3, 0.4, 0.3], num_samples )
resMultinomial= Turing.resampleMultinomial( [0.3, 0.4, 0.3], num_samples )
resResidual   = Turing.resampleResidual( [0.3, 0.4, 0.3], num_samples )
Turing.resample( [0.3, 0.4, 0.3])
resSystematic2=Turing.resample( [0.3, 0.4, 0.3], num_samples )

@test_approx_eq_eps sum(resSystematic .== 2) (num_samples * 0.4) 1e-3*num_samples
@test_approx_eq_eps sum(resSystematic2 .== 2) (num_samples * 0.4) 1e-3*num_samples
@test_approx_eq_eps sum(resStratified .== 2) (num_samples * 0.4) 1e-3*num_samples
@test_approx_eq_eps sum(resMultinomial .== 2) (num_samples * 0.4) 1e-2*num_samples
@test_approx_eq_eps sum(resResidual .== 2) (num_samples * 0.4) 1e-2*num_samples


text:
#Initiailize ULTRAFAST according to the parameters set in main.jl
#Note that the use of const global variable is more efficient than global variables

@everywhere using ParallelDataTransfer; @everywhere using DelimitedFiles; @everywhere using IterativeSolvers
@everywhere using Statistics; @everywhere using LinearAlgebra
#Define as const global variables the parameters set in main.jl
@everywhere const alpha = Int64(@getfrom 1 Alpha)
@everywhere const L_x = Int64(@getfrom 1 Lx)
@everywhere const L_y = Int64(@getfrom 1 Ly)
@everywhere const nspins = L_x*L_y            #Number of spins in the lattice

@everywhere const J_x = Float64(@getfrom 1 Jx)
@everywhere const J_y = Float64(@getfrom 1 Jy)
@everywhere const dim = alpha*nspins + alpha  #Number of independent variational parameters

#Set magnetization sector for Monte Carlo sampling. (default mag0 = true, only sample zero magnetization states)
@everywhere const mag0 = true

@everywhere include("../neuralNetwork/init_hyperparameters.jl")
@everywhere include("../neuralNetwork/nn.jl")
@everywhere include("../neuralNetwork/WtoPar_ti.jl")
@everywhere include("../observables/S_matrix.jl")
@everywhere include("../observables/GradE.jl")
@everywhere include("../parallelization/parallelization.jl")
@everywhere include("../observables/energy.jl")
@everywhere include("../observables/wv_derivative.jl")
@everywhere include("../sampling/sampling.jl")
@everywhere include("../lattice/Heisenberg_lattice.jl")
@everywhere include("../groundState/ground_state.jl")
@everywhere include("../dynamics/dynamics.jl")
@everywhere include("../observables/spincorr.jl")
@everywhere cd("..\\ULTRAFAST\\src\\output")
#Define struct storing hyper-parameters set in main.jl
GS_HP = GS_hp(nSampleGs,nIter,stepSize)
DYN_HP = DYN_hp(nSampleDyn,totTime,stepSizeHeun,minres_)

#Store nearest neighbouring bonds
@everywhere const bonds_x, bonds_y = Lattice.InitLattice(L_x,L_y,nspins)

#Store filter matrix for efficient evaluation of wavefunction derivatives
@everywhere const Filter_ind = WvDerivativeTi.transl_inv_tool(nspins,alpha,L_x,L_y)

#Define number of samples handled by each worker
@everywhere const Nsample = partition_sample(@getfrom 1 GS_HP.samples)
@everywhere const NsampleDYN = partition_sample(@getfrom 1 DYN_HP.samples)

############################### Start simulations ##############################

#Start ground state optimization
if gs == true
    const Energy_ = Float64[]; const Variance_ = Float64[]
    @time GS_optimization(GS_HP,dim)
end

#Start dynamics
if dy == true
    runDynamics(tamedHeun,DYN_HP,WRBMInit[:])
end

max_stars_repo_name: GiamFbn/MyULTRAFAST

max_issues_repo_name: GiamFbn/MyULTRAFAST

max_forks_repo_name: GiamFbn/MyULTRAFAST


text:
#   This file is part of Reduce.jl. It is licensed under the MIT license
#   Copyright (C) 2017 Michael Reed

using Documenter, SyntaxTree, ForceImport, Reduce

makedocs(
    # options
    modules = [Reduce,SyntaxTree,ForceImport],
    doctest = false,
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    sitename = "Reduce.jl",
    authors = "Michael Reed",
    pages = Any[
        "Home" => "index.md",
        "Library" => "library.md",
        "User's Manual" => Any[
            "man/acknowledgement.md",
            "man/01-introduction.md",
            "man/02-structure.md",
            "man/03-expressions.md",
            "man/04-lists.md",
            "man/05-statements.md",
            "man/06-commands-declarations.md",
            "man/07-prefix-ops.md",
            "man/08-display.md",
            "man/09-polynomials.md",
            "man/10-properties.md",
            "man/11-substitution.md",
            "man/12-file-io.md",
            "man/13-interactive.md",
            "man/14-matrix.md",
            "man/15-procedures.md",
            "man/16-packages.md",
            "man/17-symbolic.md",
            "man/18-physics.md",
            "man/19-rlisp.md",
            "man/20-maintaining.md",
            "man/A-reserved.md",
            "man/B-bibliography.md",
            "man/C-changelog.md",
            "man/index.md"
            ]
        ]
)

deploydocs(
    repo   = "github.com/chakravala/Reduce.jl.git",
    target = "build",
    deps = nothing,
    make = nothing
)