text:
# Date: July 1, 2018
# Author: Bernhard König

############################################################
# Title: Boosting Model for the French MTPL Data
############################################################

# We suggest that you first go through the frenchMTPLdata_single_tree.jl tutorial before doing this tutorial

##############################
# Load packages
##############################
t0 = time_ns()
cd(joinpath(GLOBAL_julia_code_folder, "DecisionTrees.jl"))
# @info("You may want to run 'pkg> instantiate' when you first run this. Use ] to enter the package mode.")

using Distributed
# Distributed.@everywhere using Revise
Distributed.@everywhere import CSV
Distributed.@everywhere import DataFrames
Distributed.@everywhere import DataFrames: DataFrame

Distributed.@everywhere using DecisionTrees

tela = (time_ns() - t0) / 1e9
@show tela # precompilation may take some time

##############################
# Read the data
##############################
datafile = joinpath("data", "freMTPL2", "freMTPL2.csv")
@assert isfile(datafile);
@time fullData = CSV.read(datafile, DataFrame, strict=true);

# the data is described in the R package CASdatasets
# also the Noll, Salzmann, Wüthrich Paper has some descriptive graphs of the data.

# add AreaInteger as new variable, i.e. A:F -> 1:6 (as suggested in the paper of Wüthrich, Noll, Salzmann)
areasSorted = sort(unique(fullData[:Area]))
AreaInteger = map(x->findall((in)([x]), areasSorted)[1], fullData[:Area])
fullData[:AreaInteger] = AreaInteger

# correct for unreasonable observations
for i = 1:size(fullData, 1)
    fullData[:ClaimNb][i] = min(4, fullData[:ClaimNb][i])
    fullData[:Exposure][i] = min(1.0, fullData[:Exposure][i])
end

# set independent variables
selected_explanatory_vars = ["Area","AreaInteger","VehPower","VehAge","DrivAge","BonusMalus","VehBrand","VehGas","Density","Region"]
# note: we refrain from using the exposure as explanatory variable, as this is not meaningful in practical applications
# where one aims to predict the claim frequency of a given policy (without knowing how long that policy will remain within a certain portfolio)

# check type of each column
for x in selected_explanatory_vars
    println(x)
    println(eltype(fullData[Symbol(x)]))
    println("first ten values")
    print(fullData[1:10,Symbol(x)])
    println("")
end

# consider unique elements for each variable
for x in 1:size(fullData, 2)
    println(x)
    @show size(unique(fullData[:,x]))
end

##############################
# Prepare the data
##############################

# dtmtable::DTMTable is a custom data format for DecisionTrees
# sett are the model settings
# dfprepped is an intermediary dataframe which is generally not needed
dtmtable, sett, dfprepped = prepare_dataframe_for_dtm!(fullData, keycol="IDpol", trnvalcol="trnTest", numcol="ClaimNb", denomcol="Exposure", weightcol="Exposure", independent_vars=selected_explanatory_vars);
# consider
fieldnames(typeof(dtmtable))
dtmtable.key # an identifier (String type), ideally it is a unique identifier for each row
dtmtable.numerator # a vector
dtmtable.denominator # a vector
dtmtable.weight # a vector
# each of the above have the same length

dtmtable.trnidx # an index vector for the training data. This is an integer vector
# therefore we generally have
length(dtmtable.trnidx) < length(dtmtable.numerator)
# dtmtable.validx is the corresponding validation data
# we trnidx and validx should generally be a partition of 1:n (where n==length(dtmtable.weights))

dtmtable.features # explanatory variables
# note that the data is generally stored in a compressed format
# for numerical variables, the algorithm considers sett.maxSplittingPoints  (default 250) splitting points which are uniformly spaced
# you can increase this number and re-prepare the data if you want
#=
    sett.maxSplittingPoints=100
    dtmtable=prep_data_from_df(df_prepped,sett,fn) =#
# we realize that it may be sutiable to define the splitting in a different manner (than uniformly spaced).
# this feature might be added at a later time. You can consider the function add_coded_numdata!(...) to see how splitting points are chosen

originalTrnValIndex = deepcopy(fullData[:trnTest])

# Redefine trn and val data sets
# if you prefer train on X% of the data, you can use this function to adjust the training set. By default it is sampled randomly
# resample_trnvalidx!(dtmtable,.85)

##############################
# Define model SETTINGS
##############################

# the minimum weight of each individual tree shall be 3% of the exposure per leaf
# You can set minWeight to any positive value, e.g. sett.minWeight=20000
# if minWeight is set to a value in [-1,0] its absolute value is interpreted as a percentage
# thus with sett.minWeight=-0.03 we define the min weight of each leaf as 3% of the training data

# model_type: currently only "build_tree" and "boosted_tree" are supported
# we might add the implementation of bagging at a later stage

# iterations is the number of trees

# learningRate is the moderation factor which is also know as shrinkage or learning rate

# subsampling_features_prop must be between 0 and 1 and defines the subsampling probability of the predictors for each tree
# i.e. with subsampling_features_prop=0.5 each tree will only use half of the predictors

# calculatePoissonError, is a boolean which we need to enable here in order for output to show the poisson error of the estimates
updateSettingsMod!(sett,minWeight=-0.03,model_type="boosted_tree",iterations=40,learningRate=0.05,subsampling_features_prop=.7,calculatePoissonError=true)

##############################
# Run single Boosting Model
##############################

# This model might take few minutes to run.
# resM is the resulting Model
resultingFiles, resM = dtm(dtmtable, sett)
# consider the resulting Excel file for a descritipon of the model
# it should be located here:
@show resultingFiles[2]

# if you prefer to analyze the data in Julia instead of excel, consider
sheetnames = map(x->x.name, resM.exceldata.sheets)
# the following is the ModelStatistics sheets as a Dataframe
resM.exceldata.sheets[3].data
# this is the lift in the training data for each iteration of the boosting model
resM.exceldata.sheets[3].data[:x6][1:sett.iterations + 1]

############################################################
# Investigate the predictions (fitted values)
############################################################

# DecisionTrees Boosting Models provide a number of different estimates
# 1. Fitted values based on the "raw relativities"
# this is simply the combination of all decision trees

# 2. Smoothed fitted values per Score.
# the user can specify the number of scores for the output (see sett.nScores which defaults to 1000)
# DecisionTrees runs a moving average over the observed data to smooth the noise

# 3. Unsmoothed fitted values per Score.
# As the score is merely an aggregation of the resulting 'relativities' (i.e. the spread of the fitted frequencies by the ensbemble)
# one can calculate the observed frequency for each score. Here, no smoothing is performed

# If you want, you can consider the function below to get the fitted value
# rawFittedValues,smoothedFittedValuesPerScore,unsmoothedFittedValuesPerScore=getFittedValues(reM)
# in the following lines we derive the fitted values from the resM struct

# We will show the different fitted values in the following

# 1. Fitted values based on the "raw relativities"
# for the fitted values, consider:
fitted = resM.meanobserved .* resM.rawrelativities
# resM.meanobserved is the mean observed value of the training data
# rawrelativities is the estimated relative deviation from the mean for each observation
# Note:as for the input data the fitted vector corresponds to all the data.
# dtmtable.validx determines which observations correspond to the validation data

# Let us calculate the Poisson Error for these estimates
errs = poissonError(dtmtable.numerator, fitted .* dtmtable.denominator);
trnAvgError = sum(errs[dtmtable.trnidx]) / length(dtmtable.trnidx) # should be around 0.3104379280001656
valAvgError = sum(errs[dtmtable.validx]) / length(dtmtable.validx) # should be around 0.32060219905526594

# 2. Smoothed fitted values per Score.
scores = resM.scores
fittedThroughScores = zeros(Float64, length(scores))
for i = 1:length(scores)
    @inbounds sc = scores[i]
    @inbounds fittedThroughScores[i] = resM.ScoreToSmoothedEstimate[sc]
end

fitted = fittedThroughScores
# notably for these fitted values we have less unique values
@show unique(fitted)
# this should correspond to sett.nScores
errs = poissonError(dtmtable.numerator, fitted .* dtmtable.denominator);
trnAvgError = sum(errs[dtmtable.trnidx]) / length(dtmtable.trnidx) # 0.30909682928383053
valAvgError = sum(errs[dtmtable.validx]) / length(dtmtable.validx) # 0.3192856166040531

# 3. (raw) Observed frequency per Score
# As the score is merely an aggregation of the resulting 'relativities' (i.e. the spread of the fitted frequencies by the ensbemble)
# one can calculate the observed frequency for each score. Here, no smoothing is performed
scores = resM.scores
rawFittedThroughScores = zeros(Float64, length(scores))
for i = 1:length(scores)
    @inbounds sc = scores[i]
    @inbounds rawFittedThroughScores[i] = resM.rawObservedRatioPerScore[sc]
end

# Again, let us consider the Poisson Error
fitted = rawFittedThroughScores
errs = poissonError(dtmtable.numerator, fitted .* dtmtable.denominator);
trnAvgError = sum(errs[dtmtable.trnidx]) / length(dtmtable.trnidx) # 0.30523053930420263
valAvgError = sum(errs[dtmtable.validx]) / length(dtmtable.validx) # 0.31953356919603665

# Unseen data
# Note: For any out of sample data (which MUST have exactly the same structure as dtmtable.features!)
# you can also use the predict function
estimatedFrequencieForFirst20Obs = predict(resM, dtmtable.features[1:20,:])
# getting new/unseen data into the right format is not trivial at the moment. It is easiset to use prepare_dataframe_for_dtm! to the
# whole data set and then discrard possible hold out observations  (or mark as validation through validx)

# Let us try a more granular model with more iterations and a higher learning rate
updateSettingsMod!(sett,minWeight=-0.01,iterations=100,learningRate=0.1,subsampling_features_prop=1.0,calculatePoissonError=true)
resultingFiles, resMGranular = dtm(dtmtable, sett)

# todo write down error metrics of this model

# Let us calculate its error
fitted = resMGranular.meanobserved .* resMGranular.rawrelativities
# Let us calculate the Poisson Error for these estimates
errs = poissonError(dtmtable.numerator, fitted .* dtmtable.denominator);
trnAvgError = sum(errs[dtmtable.trnidx]) / length(dtmtable.trnidx) # should be around 0.296411226074372
valAvgError = sum(errs[dtmtable.validx]) / length(dtmtable.validx) # should be around 0.3129144168276613

# The model is slightly better than the previous model.
# Let us consider a cross validation of this model.
cvsampler = CVOptions(-5, 0.0, true)
# this will run 5 models and thus might take some time
statsdf, settsdf, cvModels = dtm(dtmtable, sett, cvsampler)

# todo continue here
# Let us attach the scores to the data and write it to disk.
# We will then fit a GLM in R with and without the scores as predictor
fullDataWithScores = deepcopy(fullData);
fullDataWithScores[:scores] = resMGranular.scores;
writetable(joinpath("data", "freMTPL2", "freMTPL2_withJuliaBoostingScores1000.csv"),fullDataWithScores)

############################################################
# Additional options
############################################################

# You can modify the settings directly, but should generally not do so.
sett.iterations = 4
# It is preferred to use the update function to avoid creating an 'invalid' setting (such as a negative number of iterations)
# e.g. this will throw an error (whereas sett.iterations=-3 does not)
# updateSettingsMod!(sett,iterations=-3)

# you can modify the number of scores which are derived
updateSettingsMod!(sett,nScores=10000)
# let us set it back to 1000
updateSettingsMod!(sett,nScores=1000)
# The resulting Excel ouput shows a number of graphs per score band
# if you are interested in the tails of the data, you many want to pick more granular score bands
# the default is this
sett.scorebandsstartingpoints
# you can set them manually by defining the 'start' of each band
updateSettingsMod!(sett,scorebandsstartingpoints=[1,25,50,75,100,300,500,800,900,925,950,975])
sett.scorebandsstartingpoints
# at the moment you will need to re run the hwole model if you amend nScores or scorebands
updateSettingsMod!(sett,iterations=20)
resultingFiles, resM = dtm(dtmtable, sett)

# SOME COMMENT ON MANY UNUSED AND EXPERMIANTAL options
# if you want to use the resutling tree in a different programming language you may want to
# consider these options which give you the core structure of the model in other languages
# updateSettingsMod!(sett,writeCsharpCode="true",writeVbaCode="true",writeSasCode="true")

# mention DOT graph

############################################################
# Prepare a grid search over the parameter space
############################################################

# set some default settings
updateSettingsMod!(sett,
iterations=250,
model_type="boosted_tree",
nScores="1000",
writeStatistics="true",
writeSasCode="false",
writeIterationMatrix="false",
writeResult="false",
writeCsharpCode="false",
writeVbaCode="false",
saveJLDFile="false",
saveResultAsJLDFile="false",

showProgressBar_time="true",
prroduceEstAndLeafMatrices="false",
write_dot_graph="false",

calculatePoissonError="true",
performanceMeasure="Average Poisson Error Val"
)

# performanceMeasure will determine which metric is considered for the output summary
# performanceMeasure must be in DecisionTrees.global_statsperiter_header
# you may want to consider the Excel file of any boosting to get an idea of the available options for performanceMeasure
# consider the first row of the sheet ModelStatistics of the Excel file which lists all measures
# We note that some missing (or constant) measures are either not fully implemented yet, or were not calculated (e.g. because a boolean in the settings was set to false)
@show sett.performanceMeasure
@assert in(sett.performanceMeasure, DecisionTrees.global_statsperiter_header)

# create a vector of ModelSettings
# note, we can provide an arbitrary number of vectors over which we want to loop
# the function createGridSearchSettings will copy the 'input' settings X times into a vector of settings
# X is the size of the cartesian product of all parameter vectors over which we want to loop.

settV = createGridSearchSettings(sett,
    minWeight=[-0.03,-0.01,-0.005,-0.0025]
    ,learningRate=[0.1,0.05],
    subsampling_features_prop=[.5,.7,1.0],
    crit=["poissondeviance", "difference"])
    # ,   subsampling_prop=[1.0,.5,.7]);

# consider the length(settV) which is the number of models that will be run
length(settV)

# depending on the size of settV, you may want to restart Julia with additional workers for parallelization, i.e.
# shell> julia -p 8
# to start julia with 8 workers. Depending on your CPU you may want to choose a different number

Sys.CPU_CORES # might be give you an indication of the number of workers() you could use

############################################################
# Run grid search
############################################################

@warn "This may take quiet some time depending on length(settV)=$(length(settV))"
@info "Starting grid search..."

tt0 = time_ns()
statsdf, settsdf, allmodels = dtm(dtmtable, settV)
@show ela = (-tt0 + time_ns()) / 1e9
@info ".....done"

# Consider the resulting vector of models 'allmodels'
# We apply addtional functions on these results, e.g.
fittedValuesRawRelativities = map(x->x.meanobserved .* x.rawrelativities, allmodels);

# in the following we vill evaluate the Poisson Error for the LAST iteration of the model
# Note: this is not the 'BEST' iteration, but the last one.
valAvgErrsMultirun = zeros(length(fittedValuesRawRelativities))
for i = 1:length(fittedValuesRawRelativities)
    fitted = fittedValuesRawRelativities[i]
    errs = poissonError(dtmtable.numerator, fitted .* dtmtable.denominator);
    trnAvgError = sum(errs[dtmtable.trnidx]) / length(dtmtable.trnidx)
    valAvgError = sum(errs[dtmtable.validx]) / length(dtmtable.validx)
    valAvgErrsMultirun[i] = valAvgError
end

@show bestModel = argmin(valAvgErrsMultirun)
bestSetting = settsdf[bestModel,:]
# this data should be consistent with the data in the Excel output

    # crit::SplittingCriterion # fn version of 4

############################################################
# Partial Dependence Plots
############################################################
# this is currently in development and in need of review
# NOTE: the function is very slow. This will take a few minutes
# pds=partialDependence(resM,dtmtable,listOfFeatures=[:Density,:VehGas])
# pds["Density"]


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

## dummy stub for https://github.com/JuliaBinaryWrappers/MozillaCACerts_jll.jl

module MozillaCACerts_jll

function __init__()
	global cacert = normpath(Sys.BINDIR::String, Base.DATAROOTDIR, "julia", "cert.pem")
end

end # module

ext: jl

max_stars_repo_name: vanillajonathan/julia

max_issues_repo_name: vanillajonathan/julia

max_forks_repo_name: vanillajonathan/julia


text:
bson_type(::Nothing) = null
bson_type(::Bool) = boolean
bson_type(::Int32) = int32
bson_type(::Int64) = int64
bson_type(::Float64) = double
bson_type(::String) = string
bson_type(::Vector{UInt8}) = binary
bson_type(::BSONDict) = document
bson_type(::BSONArray) = array

bson_primitive(io::IO, ::Nothing) = return
bson_primitive(io::IO, x::Union{Bool,Int32,Int64,Float64}) = write(io, x)
bson_primitive(io::IO, x::Float64) = write(io, x)
bson_primitive(io::IO, x::Vector{UInt8}) = write(io, Int32(length(x)), 0x00, x)
bson_primitive(io::IO, x::String) = write(io, Int32(sizeof(x)+1), x, 0x00)

bson_key(io::IO, k) = write(io, Base.string(k), 0x00)

function bson_pair(io::IO, k, v)
  write(io, bson_type(v))
  bson_key(io, k)
  bson_primitive(io, v)
end

function bson_doc(io::IO, doc)
  buf = IOBuffer()
  for (k, v) in doc
    bson_pair(buf, k, v)
  end
  write(buf, eof)
  bytes = read(seek(buf,0))
  write(io, Int32(length(bytes)+4), bytes)
  return
end

bson_primitive(io::IO, doc::BSONDict) = bson_doc(io, doc)
bson_primitive(io::IO, x::BSONArray) =
  bson_doc(io, [Base.string(i-1) => v for (i, v) in enumerate(x)])

# Lowering

lower(x::Primitive) = x

import Base: RefValue

ismutable(T) = !isbitstype(T)
ismutable(::Type{String}) = false

typeof_(x) = typeof(x)
typeof_(T::DataType) = T

function _lower_recursive(x, cache, refs)
  _lower(x) = applychildren!(x -> _lower_recursive(x, cache, refs), lower(x)::Primitive)
  ismutable(typeof_(x)) || return RefValue{Any}(_lower(x))
  if haskey(cache, x)
    if !any(y -> x === y, refs)
      push!(refs, cache[x])
    end
    return cache[x]
  end
  cache[x] = RefValue{Any}(nothing)
  val = applychildren!(x -> _lower_recursive(x, cache, refs), lower(x)::Primitive)
  cache[x].x == nothing && (cache[x].x = val)
  return cache[x]
end

stripref(x::RefValue) = stripref(x.x)
stripref(x) = applychildren!(stripref, x)

function lower_recursive(y)
  cache = IdDict()
  backrefs = []
  x = _lower_recursive(y, cache, backrefs).x
  isempty(backrefs) || (x[:_backrefs] = Any[x.x for x in backrefs])
  foreach((ix) -> (ix[2].x = BSONDict(:tag=>"backref",:ref=>ix[1])), enumerate(backrefs))
  return stripref(x)
end

# Interface

bson(io::IO, doc::AbstractDict) = bson_primitive(io, lower_recursive(doc))

bson(path::String, doc::AbstractDict) = open(io -> bson(io, doc), path, "w")

bson(path::String; kws...) = bson(path, Dict(kws))


text:
using JuDGE, DelimitedFiles, DataFrames, Glob, CSV

treeFileName = "7node-tree2.csv"
# treeFileName = "341node-tree.csv"

# Load the n-node tree and ASCII-print the tree to the REPL
mytree, data = tree_from_file(joinpath(@__DIR__,treeFileName))

# Define a function that gets the parent for a given node
parentfn = JuDGE.parent_builder(mytree)

# Get all nodes, and the corresponding parent of each node
rawNodes = JuDGE.collect(mytree)
rawParents = parentfn.(rawNodes)

# Get all nodenames
node = []
for n in rawNodes; push!(node, n.name) end

# Tree hierarchy
hierarchy = vcat(hcat("", reshape(node, (1,length(node)))), hcat(node, zeros(Int64, length(node),length(node))))

for i in 1:length(rawNodes)
    # Get the name, and the name of the parent of the 'i'th rawNode, accounting
    # for the root
    nodename = rawNodes[i].name
    if rawParents[i] != nothing
        parname = rawParents[i].name
        idx = findall(x->x==rawParents[i], rawNodes)[1]
        hierarchy[idx + 1, i + 1] = 1
    end
end

# Write to CSV. This can be imported to igraph in R as an adjacency matrix
CSV.write(joinpath(@__DIR__, "node_adj_matrix.csv"), DataFrame(hierarchy), header=false)

max_stars_repo_name: calroughan/EPSILON

max_issues_repo_name: calroughan/EPSILON

max_forks_repo_name: calroughan/EPSILON


text:
module ITensorVisualizationCore

using Compat

export @visualize,
  @visualize!,
  @visualize_noeval,
  @visualize_noeval!,
  @visualize_sequence,
  @visualize_sequence_noeval

# Visualizing ITensor networks
include("visualize_macro.jl")

end


text:
## Julia program for Truncated normal distribution
## author: weiya <szcfweiya@gmail.com>
## date: 2018-08-22

# Truncated normal distribution
function rtrunormal(T, mu, sigma, mu_down)
    x = ones(T)
    z = ones(T+1)
    # set initial value of z
    z[1] = rand()
    if mu < mu_down
        z[1] = z[1] * exp(-0.5 * (mu - mu_down)^2 / sigma^2)
    end
    for t = 1:T
        x[t] = rand() * (mu - mu_down + sqrt(-2*sigma^2*log(z[t]))) + mu_down
        z[t+1] = rand() * exp(-(x[t] - mu)^2/(2*sigma^2))
    end
    return(x)
end

## example
rtrunormal(1000, 1.0, 1.0, 1.2)

max_stars_repo_name: szcf-weiya/MCMC

max_issues_repo_name: szcf-weiya/MCMC

max_forks_repo_name: szcf-weiya/MCMC


text:
# This file is a part of GERDAMetadata.jl, licensed under the MIT License (MIT).

export FileKey

struct FileKey
    setup::Symbol
    run::Int
    time::Int64
    category::Symbol
end

const filekey_expr = r"^([A-Za-z_]+)-run([0-9]{4})-([0-9]{8}T[0-9]{6}Z)-([A-Za-z_]+)$"
const filekey_wildcard_expr = r"^(%|[A-Za-z_]+)-(%|run[0-9]{4})-(%|[0-9]{8}T[0-9]{6}Z)-(%|[A-Za-z_]+)$"

const timestamp_format = dateformat"yyyymmddTHHMMSSZ"

timestamp2unix(t::String) = Int64(Dates.datetime2unix(DateTime(t, timestamp_format)))

unix2timestamp(t::Integer) = Dates.format(Dates.unix2datetime(t), timestamp_format)

import Base.==
==(a::FileKey, b::FileKey) = a.setup == b.setup && a.run == b.run && a.time == b.time && a.category == b.category

filekey_run_str(key::FileKey) = "run$(lpad(string(key.run), 4, string(0)))"

function FileKey(s::AbstractString)
    m = match(filekey_expr, s)
    if (m == nothing)
        throw(ArgumentError("String \"$s\" does not represent a valid file key"))
    else
        x = (m::RegexMatch).captures
        FileKey(
            Symbol(x[1]),
            parse(Int, x[2]),
            Int64(Dates.datetime2unix(DateTime(x[3], timestamp_format))),
            Symbol(x[4])
        )
    end
end

Base.convert(::Type{FileKey}, s::AbstractString) = FileKey(s)

function Base.print(io::IO, key::FileKey)
    print(io, key.setup)
    print(io, "-", filekey_run_str(key))
    print(io, "-", unix2timestamp(key.time))
    print(io, "-", key.category)
end

Base.show(io::IO, key::FileKey) = print(io, "FileKey(\"$(string(key))\")")

function Base.occursin(key::FileKey, pattern::AbstractString)
    mk = match(filekey_wildcard_expr, string(key))
    mp = match(filekey_wildcard_expr, pattern)

    if (mp == nothing)
        throw(ArgumentError("Not a valid file key pattern: \"$pattern\""))
    else
        all((a == b || b == "%" for (a,b) in zip(mk.captures, mp.captures)))
    end
end


text:
# Chemfiles.jl, a modern library for chemistry file reading and writing
# Copyright (C) Guillaume Fraux and contributors -- BSD license

export remove!, bonds_count, angles_count, dihedrals_count, impropers_count,
bonds, angles, dihedrals, impropers, add_bond!, remove_bond!, add_residue!,
residues, are_linked, count_residues

"""
Create an empty ``Topology``.
"""
function Topology()
    return Topology(lib.chfl_topology())
end

"""
Get a copy of the ``Topology`` of the given ``frame``.
"""
function Topology(frame::Frame)
    return Topology(lib.chfl_topology_from_frame(frame.handle))
end

"""
Free the allocated memory for the ``Topology`` object.
"""
function _free(topology::Topology)
    lib.chfl_topology_free(topology.handle)
end

"""
Get the ``Topology`` size, i.e. the current number of atoms.
"""
function Base.size(topology::Topology)
    n = Ref{UInt64}(0)
    _check(
        lib.chfl_topology_atoms_count(topology.handle, n)
    )
    return n[]
end

"""
Add an ``atom`` at the end of a ``topology``.
"""
function add_atom!(topology::Topology, atom::Atom)
    _check(
        lib.chfl_topology_add_atom(topology.handle, atom.handle)
    )
    return nothing
end

"""
Remove the atom at the given ``index`` from a ``topology``.
"""
function remove!(topology::Topology, index::Integer)
    _check(
        lib.chfl_topology_remove(topology.handle, UInt64(index))
    )
    return nothing
end

"""
Get the number of bonds in the ``topology``.
"""
function bonds_count(topology::Topology)
    n = Ref{UInt64}(0)
    _check(
        lib.chfl_topology_bonds_count(topology.handle, n)
    )
    return n[]
end

"""
Get the number of angles in the ``topology``.
"""
function angles_count(topology::Topology)
    n = Ref{UInt64}(0)
    _check(
        lib.chfl_topology_angles_count(topology.handle, n)
    )
    return n[]
end

"""
Get the number of dihedral angles in the ``topology``.
"""
function dihedrals_count(topology::Topology)
    n = Ref{UInt64}(0)
    _check(
        lib.chfl_topology_dihedrals_count(topology.handle, n)
    )
    return n[]
end

"""
Get the number of improper angles in the ``topology``.
"""
function impropers_count(topology::Topology)
    n = Ref{UInt64}(0)
    _check(
        lib.chfl_topology_impropers_count(topology.handle, n)
    )
    return n[]
end

"""
Get the bonds in the ``topology``, in a ``2 x bonds_count(topology)`` array.
"""
function bonds(topology::Topology)
    count = bonds_count(topology)
    result = Array{UInt64}(2, count)
    _check(
        lib.chfl_topology_bonds(topology.handle, pointer(result), count)
    )
    return result
end

"""
Get the angles in the ``topology``, in a ``3 x angles_count(topology)`` array.
"""
function angles(topology::Topology)
    count = angles_count(topology)
    result = Array{UInt64}(3, count)
    _check(
        lib.chfl_topology_angles(topology.handle, pointer(result), count)
    )
    return result
end

"""
Get the dihedral angles in the ``topology``, in a ``4 x dihedrals_count(topology)``
array.
"""
function dihedrals(topology::Topology)
    count = dihedrals_count(topology)
    result = Array{UInt64}(4, count)
    _check(
        lib.chfl_topology_dihedrals(topology.handle, pointer(result), count)
    )
    return result
end

"""
Get the improper angles in the ``topology``, in a ``4 x impropers_count(topology)``
array.
"""
function impropers(topology::Topology)
    count = impropers_count(topology)
    result = Array{UInt64}(4, count)
    _check(
        lib.chfl_topology_impropers(topology.handle, pointer(result), count)
    )
    return result
end

"""
Add a bond between the atoms ``i`` and ``j`` in the ``topology``.
"""
function add_bond!(topology::Topology, i::Integer, j::Integer)
    _check(
        lib.chfl_topology_add_bond(topology.handle, UInt64(i), UInt64(j))
    )
    return nothing
end

"""
Remove any existing bond between the atoms ``i`` and ``j`` in the ``topology``.
"""
function remove_bond!(topology::Topology, i::Integer, j::Integer)
    _check(
        lib.chfl_topology_remove_bond(topology.handle, UInt64(i), UInt64(j))
    )
    return nothing
end

"""
Add a copy of ``residue`` to this ``topology``.

The residue id must not already be in the topology, and the residue must
contain only atoms that are not already in another residue.
"""
function add_residue!(topology::Topology, residue::Residue)
    _check(
        lib.chfl_topology_add_residue(topology.handle, residue.handle)
    )
    return nothing
end

"""
Get the number of residues in the ``topology``.
"""
function count_residues(topology::Topology)
    nresidues = Ref{UInt64}(0)
    _check(
        lib.chfl_topology_residues_count(topology.handle, nresidues)
    )
    return nresidues[]
end

"""
Check if the two residues ``first`` and ``second`` from the ``topology`` are
linked together. *i.e.* if there is a bond between one atom in the first
residue and one atom in the second one.
"""
function are_linked(topology::Topology, first::Residue, second::Residue)
    result = Ref{UInt8}(0)
    _check(
        lib.chfl_topology_residues_linked(topology.handle, first.handle, second.handle, result)
    )
    return convert(Bool, result[])
end

"""
Resize the ``topology`` to hold ``natoms`` atoms. If the new number of atoms is
bigger than the current number, new atoms will be created with an empty name
and type.

If it is lower than the current number of atoms, the last atoms will be removed,
together with the associated bonds, angles and dihedrals.
"""
function Base.resize!(topology::Topology, size::Integer)
    _check(
        lib.chfl_topology_resize(topology.handle, UInt64(size))
    )
end

"""
Make a deep copy of a ``topology``.
"""
function Base.deepcopy(topology::Topology)
    handle = lib.chfl_topology_copy(topology.handle)
    return Topology(handle)
end


text:
module GalleryTests

using Test
using CellwiseValues
using TensorValues
using StaticArrays

v = [ VectorValue(i,i) for i in 1:10 ]
cv = CellValueFromArray(v)
test_index_cell_number( cv, v )

v = [ SVector(i,i) for i in 1:10 ]
cv = CellValueFromArray(v)
test_index_cell_number( cv, v )

v = CartesianIndices((3,10))
cv = CellValueFromArray(v)
test_index_cell_value( cv, v )

v = [ [i,i,i] for i in 1:10 ]
cv = CellValueFromArray(v)
test_index_cell_array( cv, v )

v = [ [i,i,i] for i in 1:10 ]
cv = CellValueFromArray(v)
test_index_cell_array( cv, v )

end # module


text:
using Terming

function main()
    # set term size and clear
    Terming.displaysize(20, 75); Terming.clear()

    # enable raw mode
    Terming.raw!(true)
    event = nothing
    while event != Terming.KeyPressedEvent(Terming.ESC)
        # read in_stream
        sequence = Terming.read_stream()
        # parse in_stream sequence to event
        event = Terming.parse_sequence(sequence)
        @show event
    end
    # disable raw mode
    Terming.raw!(false)

    return
end

main()


text:
using Test
using RiemannComplexNumbers

a = 2+3IM
b = 4-5IM

@test a+b == 6-2im
@test a/0 == ComplexInf
@test isinf(a/0)
@test Complex(a)==a
@test hash(a) == hash(Complex(a))
@test isnan((a-a)/(b-b))
@test Inf + IM == Inf + 2*im
@test ComplexInf * ComplexInf == ComplexInf
@test isnan(ComplexInf + ComplexInf)
@test isnan(0*ComplexInf)

@test 1/a == inv(a)
@test a' == conj(a)

@test one(RC) == 1
@test zero(RC) == 0

z = RC(2+3im)
@test z*z == -5 + 12im
@test z/z == 1
@test z^3 == z*z*z
@test z^(-1) == 1/z
@test sqrt(RC(-1)) == IM
@test exp(2+3im) == exp(2+3IM)
@test angle(IM) == pi/2
z = 3 - 4IM
@test abs(z*z) == 25
@test abs(z*z) == z'*z

@test abs(exp(log(-1+0IM)) + 1) < 1e-10
@test abs(log(exp(-1+2IM)) - (-1+2IM)) < 1e-10

a = RC(3.0 + 0.0im)
b = a'
@test a==b
@test !isequal(a,b)

@test ComplexNaN != NaN
@test isequal(ComplexNaN,NaN)


text:
function to_pixel2(point, dims)
    lon, lat = point[1], point[2]
    x = ((dims[2]/180.0) * (90 + (lon / 10^7)))
    y = ((dims[1]/360.0) * (180 - (lat / 10^7)))
    Point2f0(y, x)
end

using MakiE, Images
scene = Scene()
scat = scatter(rand(Point2f0, 10^6) .* 1000f0, markersize = 0.2, color = (:blue, 0.4))
gpu_buff = scene[:screen].renderlist[1][1][:position]
scene[:screen].clear = false
scat[:markersize] = 0.1
w = Point2f0(widths(screen))
scat[:positions] = rand(Point2f0, 10^6)
scat[:color] = (:white, 0.1)
center!(scene)
using ModernGL

function render_no_clear(window)
    !isopen(window) && return
    fb = GLWindow.framebuffer(window)
    wh = widths(window)
    resize!(fb, wh)
    w, h = wh
    #prepare for geometry in need of anti aliasing
    glBindFramebuffer(GL_FRAMEBUFFER, fb.id[1]) # color framebuffer
    glDrawBuffers(2, [GL_COLOR_ATTACHMENT0, GL_COLOR_ATTACHMENT1])
    # setup stencil and backgrounds
    glEnable(GL_STENCIL_TEST)
    glStencilOp(GL_KEEP, GL_KEEP, GL_REPLACE)
    glStencilMask(0xff)
    glClearStencil(0)
    glClearColor(0,0,0,0)
    glClear(GL_DEPTH_BUFFER_BIT | GL_STENCIL_BUFFER_BIT)
    glEnable(GL_SCISSOR_TEST)
    GLWindow.setup_window(window, false)
    glDisable(GL_SCISSOR_TEST)
    glColorMask(GL_TRUE, GL_TRUE, GL_TRUE, GL_TRUE)
    # deactivate stencil write
    glEnable(GL_STENCIL_TEST)
    glStencilOp(GL_KEEP, GL_KEEP, GL_REPLACE)
    glStencilMask(0x00)
    GLAbstraction.render(window, true)
    glDisable(GL_STENCIL_TEST)

    # transfer color to luma buffer and apply fxaa
    glBindFramebuffer(GL_FRAMEBUFFER, fb.id[2]) # luma framebuffer
    glDrawBuffer(GL_COLOR_ATTACHMENT0)
    glViewport(0, 0, w, h)
    GLAbstraction.render(fb.postprocess[1]) # add luma and preprocess
    glBindFramebuffer(GL_FRAMEBUFFER, fb.id[1]) # transfer to non fxaa framebuffer
    glDrawBuffer(GL_COLOR_ATTACHMENT0)
    GLAbstraction.render(fb.postprocess[2]) # copy with fxaa postprocess

    #prepare for non anti aliased pass
    glDrawBuffers(2, [GL_COLOR_ATTACHMENT0, GL_COLOR_ATTACHMENT1])

    glEnable(GL_STENCIL_TEST)
    glStencilOp(GL_KEEP, GL_KEEP, GL_REPLACE)
    glStencilMask(0x00)
    GLAbstraction.render(window, false)
    glDisable(GL_STENCIL_TEST)
    # draw strokes
    glEnable(GL_SCISSOR_TEST)
    GLWindow.setup_window(window, true)
    glDisable(GL_SCISSOR_TEST)
    glViewport(0,0, wh...)
    #Read all the selection queries
    GLWindow.push_selectionqueries!(window)
    glBindFramebuffer(GL_FRAMEBUFFER, 0) # transfer back to window
    glClearColor(0,0,0,0)
    glClear(GL_COLOR_BUFFER_BIT)
    GLAbstraction.render(fb.postprocess[3]) # copy postprocess
    return
end

function render_points(points, scene, gpu_buff, n = typemax(Int))
    screen = scene[:screen]
    dims = (widths(screen)...,)
    N = min(n, length(points))
    NBuff = length(gpu_buff)
    slices = floor(Int, N / NBuff)
    for i = 1:NBuff:N
        GLWindow.poll_glfw()
        ptr, _ = start(gpu_buff)
        for j = 1:NBuff
            @inbounds p = to_pixel2(points[i + (j - 1)], dims)
            unsafe_store!(ptr, p, j)
        end
        done(gpu_buff, (ptr, length(gpu_buff) + 1)) # unmap buffer
        render_no_clear(screen)
        GLWindow.swapbuffers(screen)
        glFinish()
    end
end
path = render_points(points, scene, gpu_buff)
screen = scene[:screen]
dims = reverse(widths(screen))


text:
# ============================================================================ #
# CORE TYPES
# ============================================================================ #
"""
Basic dict-like thing that holds multiple values for each key in a set.
"""
struct DictSet{K, V}
    data::Dict{K, Set{V}}
end

DictSet{K, V}() where {K, V} = DictSet{K, V}(Dict{K, V}())

function Base.getindex(ds::DictSet{K, V}, k::K) where {K, V}
    if haskey(ds.data, k)
        ds.data[k]
    else
        Set{V}()
    end
end

function Base.insert!(ds::DictSet{K, V}, k::K, v::V) where {K, V}
    if haskey(ds.data, k)
        push!(ds.data[k], v)
    else
        ds.data[k] = Set{V}(v)
    end
    v
end

function Base.delete!(ds::DictSet{K, V}, k::K) where {K, V}
    delete!(ds.data, k)
    ds
end

function Base.delete!(ds::DictSet{K, V}, k::K, v::V) where {K, V}
    if haskey(ds.data, k)
        s = ds.data[k]
        if length(s) > 1
            delete!(s, v)
        else
            delete!(ds.data, k)
        end
    end
    ds
end

"""

"""
struct SALNode{L, T}
    next_lbl::L
    mxsx_lbl::L

    data::Nullable{T}
end

function Base.isless(n1::SALNode, n2::SALNode)
    for f in [:next_lbl, :mxsx_lbl]
        getfield(n1, f) < getfield(n2, f) && return true
        getfield(n1, f) > getfield(n2, f) && return false
    end
    n1.data < n2.data
end

SALNode{L, T}(next_lbl::L, mxsx_lbl::L) where {L, T} =
    SALNode(next_lbl, mxsx_lbl, Nullable{T}())

Base.show(io::IO, nd::SALNode) =
    print(io, "(" * (iszero(nd.next_lbl) ? "φ" : "$(nd.next_lbl)") *
          ", $(nd.mxsx_lbl)" * (isnull(nd.data) ? "" : ", $(nd.data)") * ")")

"""

"""
mutable struct SAL{L, T}
    nodes::Vector{SortedSet{SALNode{L, T}}}

    last_label::L
end

function SAL{L, T}() where {L, T}
    SAL(SortedSet{SALNode{L, T}}[], 0)
end
SAL(::Type{L}, ::Type{T}) where {L, T} = SAL{L, T}()
SAL() = SAL(Int, Float64)

Base.show(io::IO, sal::SAL) =
    print("SAL with $(n_vertices(sal)) vertices.")

# Pretty print, showing structure.
function Base.print(io::IO, sal::SAL)
    println("SAL:")
    for i in 1:n_vertices(sal)
        print("$(lpad(i, 5)): [")
        nds = collect(sal.nodes[i])
        for el in nds[1:end-1]
            print("$el, ")
        end
        println("$(nds[end])]")
    end
end

# ============================================================================ #
# GETTERS, SETTERS, UTILS
# ============================================================================ #
get_label!(sal::SAL) = sal.last_label += 1

n_vertices(sal::SAL) = length(sal.nodes)

"""
Create a new node and insert it into the SAL.
"""
function insert_node!(sal::SAL{L, T}, ndlbl, nxlbl, sxlbl) where {L, T}
    while n_vertices(sal) < ndlbl
        push!(sal.nodes, SortedSet{SALNode{L, T}}())
    end

    push!(sal.nodes[ndlbl], SALNode{L, T}(nxlbl, sxlbl))
end

# ============================================================================ #
# OPERATIONS
# ============================================================================ #
function insert_simplex!(sal::SAL, simplex)
    # TODO: rename to Base.push?
    # TODO: remove faces of simplex!
    # Find maximal simplices that are contained in σ and remove them

    check_simplex(simplex)
    n = length(simplex)

    sxlbl = get_label!(sal)

    for (i, li) in enumerate(simplex), (j, lj) in enumerate(simplex[i+1:end])
        insert_node!(sal, li, lj, sxlbl)
    end
    insert_node!(sal, simplex[end], 0, sxlbl)
end

function Base.in(simplex, sal::SAL)
    simplex[end] > n_vertices(sal) && return false
    check_simplex(simplex)

    # TODO: filter is inefficient!
    findmatches(node, next) =
        IntSet(map(x -> x.mxsx_lbl,
                   Iterators.filter(x -> x.next_lbl == next, sal.nodes[node])))

    s = findmatches(simplex[1], simplex[2])
    for l in zip(simplex[2:end-1], simplex[3:end])
        intersect!(s, findmatches(l...))
        isempty(s) && return false
    end
    true
end

function remove_simplex!(sal::SAL, simplex)
    # find maximal simplices which contain σ
    # remove them
    # insert facets that do not contain σ
end


text:
import Distributed: Future, @spawnat, RemoteChannel, RemoteException

# Monoid instances
# ================

# this is standard Applicative combine implementation
# there is no other sensible definition for combine and hence if this does fail, it fails correctly
TypeClasses.combine(a::Future, b::Future) = mapn(combine, a, b)

"""
    orelse(future1, future2, ...)
    future1 ⊘ future2

Runs both in parallel and collects which ever result is first. Then interrupts all other futures and returns the found result.
"""
TypeClasses.orelse(a::Future, b::Future, more::Future...) = @spawnat :any begin
  n = length(more) + 2
  c = Channel(n)  # all tasks should in principle be able to send there result without waiting
  futures = (a, b, more...)
  futures′ = map(futures) do future
    @spawnat :any begin
    # fetch returns an RemoteException in case something failed, it does not throw an error
      put!(c, fetch(future))
    end
  end

  exceptions = Vector{Exception}(undef, n)
  for i in 1:n
    result = take!(c)

    if isa(result, RemoteException)  # all exceptions are Thrown{TaskFailedException}
      exceptions[i] = if isa(result.captured.ex, MultipleExceptions)
        # extract MultipleExceptions for better readability
        result.captured.ex
      else
        result
      end
      # we found an error, hence need to wait for another result
      continue
    end

    # no error, i.e. we found the one result we wanted to find
    close(c)  # we only need one result
    for t in (futures′..., futures...)
      # all tasks can be interrupted now that we have the first result
      @spawnat :any Base.throwto(t, InterruptException())
    end
    return result  # break for loop
  end
  # only if all n tasks failed
  return throw(MultipleExceptions(exceptions...))
end

# FunctorApplicativeMonad
# =======================

function TypeClasses.foreach(f, x::Future)
  f(fetch(x))
  nothing
end

TypeClasses.map(f, x::Future) = @spawnat :any f(fetch(x))
TypeClasses.pure(::Type{<:Future}, x) = @spawnat :any x

# we use the default implementation of ap which follows from flatten
# TypeClasses.ap
TypeClasses.ap(f::Future, x::Future) = @spawnat :any fetch(f)(fetch(x))
# we don't use convert for typesafety, as fetch is more flexible and also enables typechecks
# e.g. this works seamlessly to combine a Task into a Future
TypeClasses.flatmap(f, x::Future) = @spawnat :any fetch(f(fetch(x)))

# FlipTypes
# =========

# does not make much sense as this would need to execute the Future, and map over its returned value,
# creating a bunch dummy Futures within.


text:
remdepot() = joinpath(homedir(), "remnote")
knowledgebase_dir(kb::String) = joinpath(remdepot(), kb)

assistant_dir() = joinpath(remdepot(), "RemAssistant")

local_backups_dir(kb::String) = joinpath(remdepot(), kb, "backups")

function _resolve_newest_local_bk!(kb::String)

    # find newest
    lbk_dir = local_backups_dir(kb)
    !isdir(lbk_dir) && error("local backup dir not found: '$(lbk_dir)'")

    newest = ""
    for file in readdir(lbk_dir; join = true)
        !endswith(file, ".zip") && continue
        if mtime(newest) < mtime(file)
            newest = file
        end
    end

    return newest
end


text:
println("正则表达式")

pattern = r"^\s*(?:#|$)"

matched = match(pattern, "# a comment")

if matched == nothing
    println("not a comment")
else
    println("blank or a comment")
end

max_stars_repo_name: JulialangOrgCN/howtojulia

max_issues_repo_name: JulialangOrgCN/howtojulia

max_forks_repo_name: JulialangOrgCN/howtojulia


text:
using PSCOPF

using Test
using DataStructures

@testset "test_check_ptdf" begin

    network = PSCOPF.Networks.Network()
    PSCOPF.add_new_buses!(network, ["bus_1", "bus_2", "bus_3"])
    PSCOPF.add_new_branch!(network, "branch_1_2", 1000.)
    PSCOPF.add_new_branch!(network, "branch_1_3", 1000.)

    @testset "ptdf_dict" begin
        @testset "working_ptdf" begin
            ptdf = PSCOPF.Networks.PTDFDict(
                "BASECASE" => PSCOPF.Networks.PTDFValues(
                                    "branch_1_2" => SortedDict("bus_1"=>0.,
                                                                "bus_2"=>0.,
                                                                "bus_3"=>0.),
                                    "branch_1_3" => SortedDict("bus_1"=>0.,
                                                                "bus_2"=>0.,
                                                                "bus_3"=>0.)
                ),
                "branch_1_2" => PSCOPF.Networks.PTDFValues(
                                    "branch_1_2" => SortedDict("bus_1"=>0.,
                                                                "bus_2"=>0.,
                                                                "bus_3"=>0.),
                                    "branch_1_3" => SortedDict("bus_1"=>0.,
                                                                "bus_2"=>0.,
                                                                "bus_3"=>0.)
                ),
                "branch_1_3" => PSCOPF.Networks.PTDFValues(
                                    "branch_1_2" => SortedDict("bus_1"=>0.,
                                                                "bus_2"=>0.,
                                                                "bus_3"=>0.),
                                    "branch_1_3" => SortedDict("bus_1"=>0.,
                                                                "bus_2"=>0.,
                                                                "bus_3"=>0.)
                ),
            )

            @test PSCOPF.check_ptdf_case_entries(ptdf, network)
            @test PSCOPF.check_ptdf(ptdf, network)
        end

        @testset "ptdf_with_missing_n-1" begin
            ptdf = PSCOPF.Networks.PTDFDict(
                "BASECASE" => PSCOPF.Networks.PTDFValues(
                                    "branch_1_2" => SortedDict("bus_1"=>0.,
                                                                "bus_2"=>0.,
                                                                "bus_3"=>0.),
                                    "branch_1_3" => SortedDict("bus_1"=>0.,
                                                                "bus_2"=>0.,
                                                                "bus_3"=>0.)
                ),
                "branch_1_2" => PSCOPF.Networks.PTDFValues(
                                    "branch_1_2" => SortedDict("bus_1"=>0.,
                                                                "bus_2"=>0.,
                                                                "bus_3"=>0.),
                                    "branch_1_3" => SortedDict("bus_1"=>0.,
                                                                "bus_2"=>0.,
                                                                "bus_3"=>0.)
                ),
            )

            @test_broken !PSCOPF.check_ptdf_case_entries(ptdf, network)
            @test_broken !PSCOPF.check_ptdf(ptdf, network)
        end

        @testset "ptdf_with_missing_branch_values" begin
            ptdf = PSCOPF.Networks.PTDFDict(
                "BASECASE" => PSCOPF.Networks.PTDFValues(
                                    "branch_1_2" => SortedDict("bus_1"=>0.,
                                                                "bus_2"=>0.,
                                                                "bus_3"=>0.),
                                    # "branch_1_3" => SortedDict("bus_1"=>0.,
                                    #                             "bus_2"=>0.,
                                    #                             "bus_3"=>0.)
                ),
                "branch_1_2" => PSCOPF.Networks.PTDFValues(
                                    "branch_1_2" => SortedDict("bus_1"=>0.,
                                                                "bus_2"=>0.,
                                                                "bus_3"=>0.),
                                    "branch_1_3" => SortedDict("bus_1"=>0.,
                                                                "bus_2"=>0.,
                                                                "bus_3"=>0.)
                ),
                "branch_1_3" => PSCOPF.Networks.PTDFValues(
                                    "branch_1_2" => SortedDict("bus_1"=>0.,
                                                                "bus_2"=>0.,
                                                                "bus_3"=>0.),
                                    "branch_1_3" => SortedDict("bus_1"=>0.,
                                                                "bus_2"=>0.,
                                                                "bus_3"=>0.)
                ),
            )

            @test PSCOPF.check_ptdf_case_entries(ptdf, network)
            @test !PSCOPF.check_ptdf(ptdf, network)
        end

        @testset "ptdf_with_missing_bus_value" begin
            ptdf = PSCOPF.Networks.PTDFDict(
                "BASECASE" => PSCOPF.Networks.PTDFValues(
                                    "branch_1_2" => SortedDict("bus_1"=>0.,
                                                                "bus_2"=>0.,
                                                                "bus_3"=>0.),
                                    "branch_1_3" => SortedDict("bus_1"=>0.,
                                                                "bus_2"=>0.,
                                                                "bus_3"=>0.)
                ),
                "branch_1_2" => PSCOPF.Networks.PTDFValues(
                                    "branch_1_2" => SortedDict("bus_1"=>0.,
                                                                "bus_2"=>0.,
                                                                "bus_3"=>0.),
                                    "branch_1_3" => SortedDict("bus_1"=>0.,
                                                                "bus_2"=>0.,
                                                                "bus_3"=>0.)
                ),
                "branch_1_3" => PSCOPF.Networks.PTDFValues(
                                    "branch_1_2" => SortedDict("bus_1"=>0.,
                                                                "bus_2"=>0.,
                                                                "bus_3"=>0.),
                                    "branch_1_3" => SortedDict("bus_1"=>0.,
                                                                # "bus_2"=>0.,
                                                                "bus_3"=>0.)
                ),
            )

            @test PSCOPF.check_ptdf_case_entries(ptdf, network)
            @test !PSCOPF.check_ptdf(ptdf, network)
        end

    end

    @testset "check_ptdf_values" begin
        @testset "definition_example" begin
            ptdf = PSCOPF.Networks.PTDFValues(
                            "branch_1_2" => SortedDict("bus_1"=>0.,
                                                        "bus_2"=>0.,
                                                        "bus_3"=>0.),
                            "branch_1_3" => SortedDict("bus_1"=>0.,
                                                        "bus_2"=>0.,
                                                        "bus_3"=>0.)
            )

            @test PSCOPF.check_ptdf_branch_entries(ptdf, network)
            @test PSCOPF.check_ptdf_bus_entries(ptdf, network)
            @test PSCOPF.check_ptdf(ptdf, network)
        end

        @testset "missing_branch_1_3_entry" begin
            ptdf = PSCOPF.Networks.PTDFValues(
                            "branch_1_2" => SortedDict("bus_1"=>0.,
                                                        "bus_2"=>0.,
                                                        "bus_3"=>0.)
            )

            @test !PSCOPF.check_ptdf_branch_entries(ptdf, network)
            @test PSCOPF.check_ptdf_bus_entries(ptdf, network)
            @test !PSCOPF.check_ptdf(ptdf, network)
        end

        @testset "missing_bus_2_in_branch_1_2" begin
            ptdf = PSCOPF.Networks.PTDFValues(
                            "branch_1_2" => SortedDict("bus_1"=>0.,
                                                        "bus_3"=>0.),
                            "branch_1_3" => SortedDict("bus_1"=>0.,
                                                        "bus_2"=>0.,
                                                        "bus_3"=>0.)
            )

            @test PSCOPF.check_ptdf_branch_entries(ptdf, network)
            @test !PSCOPF.check_ptdf_bus_entries(ptdf, network)
            @test !PSCOPF.check_ptdf(ptdf, network)
        end

        @testset "missing_value_for_branch_1_2_bus_3" begin
            ptdf = SortedDict(
                                "branch_1_2" => SortedDict("bus_1"=>0.,
                                                            "bus_2"=>0.,
                                                            "bus_3"=>missing),
                                "branch_1_3" => SortedDict("bus_1"=>0.,
                                                            "bus_2"=>0.,
                                                            "bus_3"=>0.)
                            )
            @test_throws MethodError PSCOPF.check_ptdf(ptdf, network)
        end
    end

end

max_stars_repo_name: klorel/pscopf

max_issues_repo_name: klorel/pscopf

max_forks_repo_name: klorel/pscopf


text:
@testset "cuda/msgpass" begin
    in_channel = 10
    out_channel = 5
    N = 6
    T = Float32
    adj = [0 1 0 0 0 0
        1 0 0 1 1 1
        0 0 0 0 0 1
        0 1 0 0 1 0
        0 1 0 1 0 1
        0 1 1 0 1 0]

    struct NewCudaLayer{G} <: GNNLayer
        weight
    end
    NewCudaLayer{GRAPH_T}(m, n) = NewCudaLayer{GRAPH_T}(randn(T, m, n))
    Flux.@functor NewCudaLayer{GRAPH_T}

    (l::NewCudaLayer{GRAPH_T})(g, X) = GraphNeuralNetworks.propagate(l, g, +, X)[1]
    GraphNeuralNetworks.compute_message(n::NewCudaLayer{GRAPH_T}, x_i, x_j, e_ij) = n.weight * x_j
    GraphNeuralNetworks.update_node(::NewCudaLayer{GRAPH_T}, m, x) = m

    X = rand(T, in_channel, N) |> gpu
    g = GNNGraph(adj, ndata=X, graph_type=GRAPH_T)
    l = NewCudaLayer{GRAPH_T}(out_channel, in_channel) |> gpu

    g_ = l(g)
    @test size(node_features(g_)) == (out_channel, N)
end


text: fact(n) = prod(one(n):n)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using Test, YaoSym

@testset "constructors" begin
    @test ket"011" isa SymReg
    @test (ket"011")' isa AdjointSymReg
end


text:
function Coords2StringList(R,Z;SigDigits = 4)
    L = length(R)
    Str = " "
    for i in 1:L
        Str = Str*"$(round(R[i];sigdigits=SigDigits)), $(round(Z[i];sigdigits=SigDigits)) "
    end
    return Str
end

function DCoreParams2SVG(ID,OD,FileName;NPts = 100)
    r, z =DCoreGeom(ID/2, OD/2;NPts = NPts,PlotOn=false,UpsamplePoints = 1e4)
    z = z .+ maximum(z)
    Str = Coords2StringList(r,z;SigDigits = 4)
    open(FileName,"w") do file
        write(file,
"""
<?xml version="1.0" encoding="UTF-8" standalone="no"?>
<svg
   xmlns:dc="http://purl.org/dc/elements/1.1/"
   xmlns:cc="http://creativecommons.org/ns#"
   xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#"
   xmlns:svg="http://www.w3.org/2000/svg"
   xmlns="http://www.w3.org/2000/svg"
   xmlns:sodipodi="http://sodipodi.sourceforge.net/DTD/sodipodi-0.dtd"
   xmlns:inkscape="http://www.inkscape.org/namespaces/inkscape"
   width="$(round(OD/2))mm"
   height="$(round(maximum(abs.(z))))mm"
   viewBox="0 0 $(round(OD/2)) $(round(maximum(abs.(z))))"
   version="1.1"
   id="svg1501"
   inkscape:version="1.0.2 (e86c870879, 2021-01-15, custom)"
   sodipodi:docname="drawing.svg">
  <defs
     id="defs1495" />
  <sodipodi:namedview
     id="base"
     pagecolor="#ffffff"
     bordercolor="#666666"
     borderopacity="1.0"
     inkscape:pageopacity="0.0"
     inkscape:pageshadow="2"
     inkscape:zoom="0.35"
     inkscape:cx="400"
     inkscape:cy="560"
     inkscape:document-units="mm"
     inkscape:current-layer="layer1"
     inkscape:document-rotation="0"
     showgrid="false"
     inkscape:window-width="1624"
     inkscape:window-height="1228"
     inkscape:window-x="190"
     inkscape:window-y="190"
     inkscape:window-maximized="0" />
  <metadata
     id="metadata1498">
    <rdf:RDF>
      <cc:Work
         rdf:about="">
        <dc:format>image/svg+xml</dc:format>
        <dc:type
           rdf:resource="http://purl.org/dc/dcmitype/StillImage" />
        <dc:title></dc:title>
      </cc:Work>
    </rdf:RDF>
  </metadata>
  <g
     inkscape:label="Layer 1"
     inkscape:groupmode="layer"
     id="layer1">
    <path
       style="fill:none;stroke:#000000;stroke-width:0.5px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1;image-rendering:auto"
       d="M  $(Str) z"
       id="path2064" />
  </g>
</svg> """)
    end
end


text:
@testset "Test beads segmentation" begin
    A = [
        1 0 0 1 0
        1 0 1 1 1
    ]

    tgt = [
        0 0 0 1 0
        0 0 1 1 1
    ]

    threshold = 3
    @test beads_segment(A, threshold) == tgt

end


text:
@testset "DrivenProblem - Absence of Dissipator" begin
	# ----------------------------------------
	# DrivenProblem - empty Dissipator
  testBasis = SpinBasis(1//2)
  σm, σp = sigmam(testBasis), sigmap(testBasis)

	WP, mode = WavePacket(SoftBoxMode(), Coherent(1.0)), SoftBoxMode()

  problemA = WaveguideProblem((σm+σp,     σm, spindown(testBasis)), WP, 1.)
  solA     = solve(problemA)[2]

  problemB = WaveguideProblem((σm+σp, [], σm, spindown(testBasis)), WP, 1.)
  solB     = solve(problemB)[2]

  @test sum(tracedistance.(solA, solB)) ≈ 0
end

@testset "ScatterProblem - Absence of Dissipator" begin
	# ----------------------------------------
	# ScatterProblem - empty Dissipator
  testBasis = SpinBasis(1//2)
  σm, σp = sigmam(testBasis), sigmap(testBasis)

	WP, mode = WavePacket(SoftBoxMode(), Coherent(1.0)), SoftBoxMode()

  problemA = WaveguideProblem((σm+σp,     σm, spindown(testBasis)), WP, mode, 1.)
  solA     = solve(problemA)[2]

  problemB = WaveguideProblem((σm+σp, [], σm, spindown(testBasis)), WP, mode, 1.)
  solB     = solve(problemB)[2]

  @test sum(tracedistance.(solA, solB)) ≈ 0
end

@testset "DrivenProblem - ContinuousWave" begin
	# ----------------------------------------
	# ScatterProblem - empty Dissipator
  testBasis = SpinBasis(1//2)
  σm, σp = sigmam(testBasis), sigmap(testBasis)

  problemA = WaveguideProblem((σm+σp, σm, spindown(testBasis)), ContinuousWave(     1.0), 1.)
  solA     = solve(problemA)[2]

  problemB = WaveguideProblem((σm+σp, σm, spindown(testBasis)), ContinuousWave(t -> 1.0), 1.)
  solB     = solve(problemB)[2]

  @test sum(tracedistance.(solA, solB)) ≈ 0
end

@testset "ScatterProblem - ContinuousWave" begin
	# ----------------------------------------
	# ScatterProblem - empty Dissipator
  testBasis = SpinBasis(1//2)
  σm, σp = sigmam(testBasis), sigmap(testBasis)

	mode = SoftBoxMode()

  problemA = WaveguideProblem((σm+σp, σm, spindown(testBasis)), ContinuousWave(     1.0), mode, 1.)
  solA     = solve(problemA)[2]

  problemB = WaveguideProblem((σm+σp, σm, spindown(testBasis)), ContinuousWave(t -> 1.0), mode, 1.)
  solB     = solve(problemB)[2]

  @test sum(tracedistance.(solA, solB)) ≈ 0
end

@testset "Chirality of ScatterProblem" begin
	# ----------------------------------------
	# DrivenProblem - empty Dissipator
  testBasis = SpinBasis(1//2)
  σm, σp = sigmam(testBasis), sigmap(testBasis)

	WP, mode = WavePacket(SoftBoxMode(), Coherent(1.0)), SoftBoxMode()

	problemA = WaveguideProblem((σm+σp, [], σm, spindown(testBasis)), WP, 1.)
  solA     = solve(problemA, reltol = 1e-10)[2]

  problemB = WaveguideProblem((σm+σp, [], σm, spindown(testBasis)), WP, mode, 1.)
  solB     = ptrace.(solve(problemB, reltol = 1e-10)[2], 2)

  @test sum(tracedistance.(solA, solB)) < 1e-5
end

@testset "DrivenProblem - Mollow transform" begin
	# ----------------------------------------
	# Compare the Mollow transform to a manual coherent state
  testBasis = SpinBasis(1//2)
  σm, σp = sigmam(testBasis), sigmap(testBasis)

  problemA = WaveguideProblem((σm+σp, [], σm, spindown(testBasis)),
    WavePacket(SoftBoxMode(τ = 1.0), ArbitraryState(coherentstate(FockBasis(12), 1.0).data)), 3.
  )
  solA     = ptrace.(solve(problemA, reltol = 1e-10)[2], 1)

  problemB = WaveguideProblem((σm+σp, [], σm, spindown(testBasis)),
		WavePacket(SoftBoxMode(τ = 1.0), Coherent(1.0)), 3.
	)
  solB     = solve(problemB, reltol = 1e-10)[2]

  @test sum(tracedistance.(solA, solB)) < 1e-5
end

@testset "ScatterProblem - Mollow transform" begin
	# ----------------------------------------
	# Compare the Mollow transform to a manual coherent state
  testBasis = SpinBasis(1//2)
  σm, σp = sigmam(testBasis), sigmap(testBasis)

  problemA = WaveguideProblem((σm+σp, [], σm, spindown(testBasis)),
    WavePacket(SoftBoxMode(τ = 1.0), ArbitraryState(coherentstate(FockBasis(12), 1.0).data)),
		SoftBoxMode(τ = 2.0), 3.
  )
  solA     = ptrace.(solve(problemA, Nouts = [6], reltol = 1e-10)[2], 1)

  problemB = WaveguideProblem((σm+σp, [], σm, spindown(testBasis)),
		WavePacket(SoftBoxMode(τ = 1.0), Coherent(1.0)),
		SoftBoxMode(τ = 2.0), 3.
	)
  solB     = solve(problemB, Nouts = [6], reltol = 1e-10)[2]

  @test sum(tracedistance.(solA, solB)) < 1e-4
end

@testset "ContinuesWave(::Function) - normalisable" begin
	# ----------------------------------------
	# Show that ContinuesWave(::Function) creates correct state
	testBasis = SpinBasis(1//2)
	σm, σp = sigmam(testBasis), sigmap(testBasis)

	ts = range(-5.0, 5.0, length = 101)
	problemA = WaveguideProblem((σm+σp, [], σm, spindown(testBasis)),
    WavePacket(GaussMode(), Coherent( (π)^(1/4) )), ts
  )
	solA = solve(problemA, Nouts = [6], reltol = 1e-10)[2]

	problemB = WaveguideProblem((σm+σp, [], σm, spindown(testBasis)),
    ContinuousWave(t -> exp(-t^2 / 2)), ts
  )
	solB = solve(problemB, Nouts = [6], reltol = 1e-10)[2]

	@test sum(tracedistance.(solA, solB)) < 1e-12
end

using GSL, SpecialFunctions
@testset "ContinuesWave(::Function) - displace output" begin
	# ----------------------------------------
	# Show that displacing the output generates correct result
	function displaceProper(basis, α)
	  operator = DenseOperator(basis)

	  αc, α2, eα2 = -conj(α), abs2(α), exp(-0.5abs2(α))

	  for m ∈ 1:basis.N
	    operator.data[m, m] = eα2 * sf_laguerre_n(m - 1, 0, α2)
	    for n ∈ m:basis.N+1
	      c = √(gamma(m)/gamma(n)) * eα2 * sf_laguerre_n(m - 1, n-m, α2)
	      operator.data[m, n] = c * αc^(n-m)
	      operator.data[n, m] = c * α^(n-m)
	    end
	  end
	  return operator
	end

	testBasis = SpinBasis(1//2)
	σm, σp = sigmam(testBasis), sigmap(testBasis)
	α = 1.0 + rand()

	problemA = WaveguideProblem((σm+σp, [], σm, spindown(testBasis)),
    ContinuousWave(α), HardBoxMode(t₀ = 1.0), 2.0
  )
	solA = ptrace(solve(problemA, Nouts = [ceil(Int, α^2 + 10α)], reltol = 1e-10)[2][end], 1)

	problemB = WaveguideProblem((σm+σp, [], σm, spindown(testBasis)),
    ContinuousWave(α), HardBoxMode(t₀ = 1.0), 2.0, true
  )
	solB = ptrace(solve(problemB, Nouts = [ceil(Int, α^2 + 10α)], reltol = 1e-10)[2][end], 1)

	D = displaceProper(basis(solB), α)
	solB = D * solB * D'

	@test tracedistance(solA, solB) |> real < 1e-5
end

@testset "Controlling the output Fock space" begin
	# ----------------------------------------
	# Setting the cutoff vs. setting `Nouts`
  testBasis = SpinBasis(1//2)
  σm, σp = sigmam(testBasis), sigmap(testBasis)

  problemA = WaveguideProblem((σm+σp, [], σm, spindown(testBasis)),
    WavePacket(SoftBoxMode(τ = 1.0), Coherent(1.0, 10)),
		SoftBoxMode(τ = 2.0), 3.
  )
  solA     = solve(problemA, reltol = 1e-10)[2]

  problemB = WaveguideProblem((σm+σp, [], σm, spindown(testBasis)),
		WavePacket(SoftBoxMode(τ = 1.0), Coherent(1.0)),
		SoftBoxMode(τ = 2.0), 3.
	)
  solB     = solve(problemB, Nouts = [10], reltol = 1e-10)[2]

  @test sum(tracedistance.(solA, solB)) < 1e-10
end

@testset "DrivenProblem - DisplacedFock" begin
	# ----------------------------------------
	# Compare the `DisplacedFock` implementation to the manually created version
  testBasis = SpinBasis(1//2)
  σm, σp = sigmam(testBasis), sigmap(testBasis)

	aᵀ, ψα = create(FockBasis(13)), coherentstate(FockBasis(13), 1.0)
	data = (aᵀ * ψα - ψα).data
  problemA = WaveguideProblem((σm+σp, [], σm, spindown(testBasis)),
    WavePacket(SoftBoxMode(τ = 1.0), ArbitraryState(data)), 3.
  )
  solA     = ptrace.(solve(problemA, reltol = 1e-10)[2], 1)

  problemB = WaveguideProblem((σm+σp, [], σm, spindown(testBasis)),
		WavePacket(SoftBoxMode(τ = 1.0), DisplacedFock(1.0, 1)), 3.
	)
  solB     = ptrace.(solve(problemB, reltol = 1e-10)[2], 1)

  @test sum(tracedistance.(solA, solB)) < 1e-5
end

@testset "DrivenProblem - DisplacedArbitraryState" begin
	# ----------------------------------------
	# Compare the `DisplacedArbitraryState` implementation to the manually created version
  testBasis = SpinBasis(1//2)
  σm, σp = sigmam(testBasis), sigmap(testBasis)

	aᵀ, ψα = create(FockBasis(13)), coherentstate(FockBasis(13), 1.0)
	data = (aᵀ * ψα - ψα).data
  problemA = WaveguideProblem((σm+σp, [], σm, spindown(testBasis)),
    WavePacket(SoftBoxMode(τ = 1.0), ArbitraryState(data)), 3.
  )
  solA     = ptrace.(solve(problemA, reltol = 1e-10)[2], 1)

  problemB = WaveguideProblem((σm+σp, [], σm, spindown(testBasis)),
		WavePacket(SoftBoxMode(τ = 1.0), DisplacedArbitraryState(1.0, [0.0im, 1.0])), 3.
	)
  solB     = ptrace.(solve(problemB, reltol = 1e-10)[2], 1)

  @test sum(tracedistance.(solA, solB)) < 1e-5
end

@testset "ScatterProblem - DisplacedFock" begin
	# ----------------------------------------
	# Compare the `DisplacedFock` implementation to the manually created version
  testBasis = SpinBasis(1//2)
  σm, σp = sigmam(testBasis), sigmap(testBasis)

	aᵀ, ψα = create(FockBasis(13)), coherentstate(FockBasis(13), 1.0)
	data = (aᵀ * ψα - ψα).data
  problemA = WaveguideProblem((σm+σp, [], σm, spindown(testBasis)),
    WavePacket(SoftBoxMode(τ = 1.0), ArbitraryState(data)), SoftBoxMode(τ = 2.0), 3.
  )
  solA     = ptrace.(solve(problemA, reltol = 1e-10)[2], 1)

  problemB = WaveguideProblem((σm+σp, [], σm, spindown(testBasis)),
		WavePacket(SoftBoxMode(τ = 1.0), DisplacedFock(1.0, 1)), SoftBoxMode(τ = 2.0), 3.
	)
  solB     = ptrace.(solve(problemB, Nouts = [13], reltol = 1e-10)[2], 1)

  @test sum(tracedistance.(solA, solB)) < 1e-4
end

@testset "DrivenProblem - DisplacedArbitraryState" begin
	# ----------------------------------------
	# Compare the `DisplacedArbitraryState` implementation to the manually created version
  testBasis = SpinBasis(1//2)
  σm, σp = sigmam(testBasis), sigmap(testBasis)

	aᵀ, ψα = create(FockBasis(13)), coherentstate(FockBasis(13), 1.0)
	data = (aᵀ * ψα - ψα).data
  problemA = WaveguideProblem((σm+σp, [], σm, spindown(testBasis)),
    WavePacket(SoftBoxMode(τ = 1.0), ArbitraryState(data)), SoftBoxMode(τ = 2.0), 3.
  )
  solA     = ptrace.(solve(problemA, reltol = 1e-10)[2], 1)

  problemB = WaveguideProblem((σm+σp, [], σm, spindown(testBasis)),
		WavePacket(SoftBoxMode(τ = 1.0), DisplacedArbitraryState(1.0, [0.0im, 1.0])),
		SoftBoxMode(τ = 2.0), 3.
	)
  solB     = ptrace.(solve(problemB, Nouts = [13], reltol = 1e-10)[2], 1)

  @test sum(tracedistance.(solA, solB)) < 1e-4
end


text:
module Operations

using HyperRectangles

export update!, points

function update!{T, N}(b::HyperRectangle{T, N}, v)
    for i = 1:N
        b.max[i] = max(b.max[i], v[i])
        b.min[i] = min(b.min[i], v[i])
    end
end

function points{T}(rect::HyperRectangle{T, 3})
    Vector{T}[rect.min,
        [rect.min[1],rect.min[2],rect.max[3]],
        [rect.min[1],rect.max[2],rect.min[3]],
        [rect.min[1],rect.max[2],rect.max[3]],
        [rect.max[1],rect.min[2],rect.min[3]],
        [rect.max[1],rect.min[2],rect.max[3]],
        [rect.max[1],rect.max[2],rect.min[3]],
        rect.max]
end

function points{T}(rect::HyperRectangle{T, 2})
    Vector{T}[rect.min,
        [rect.min[1],rect.max[2]],
        [rect.max[1],rect.min[2]],
        rect.max]
end

end # module


text:
using Random
using LightGraphs

mutable struct Agent
    id::Int32
    values::AbstractArray
    vote::Int8
    energy::Float32
end

function initNetwork(N)
    #Initialize agents
    global nodes = []
    for i in 1:N
        push!(nodes,Agent(i,rand(0:10,(5)),rand(1:10),0))
    end

    #=
    #Generate graph with N
    A = bitrand((N,N)) #Adjacency matrix #TODO:Make sparse
    for i in 1:N
        A[i,i] = 0 #Make sure nodes are not be connected to themselves
    end
    global Network = DiGraph(A) #Generate graph
    A = nothing #Clear A

    #We are now implementing a Barabasi-Albert graph, instead of a random one
    =#

    #Generate Barabasi graph with N nodes, 3 conntections each, 10 initial nodes
    global Network = barabasi_albert(N, 10, 3, seed=1, is_directed=true)
    #global Network = erdos_renyi(N, 4, is_directed=true)

    #Initialize energy
    for i in 1:N
        dE(i)
    end
    computeEnergy()

    #Count initial preferences
    trackPreference()
end

max_stars_repo_path: Scripts/Prototype 2/init.jl

max_stars_repo_name: luca-mm/Bachelor-Thesis

max_issues_repo_path: Scripts/Prototype 2/init.jl

max_issues_repo_name: luca-mm/Bachelor-Thesis

max_forks_repo_path: Scripts/Prototype 2/init.jl

max_forks_repo_name: luca-mm/Bachelor-Thesis


text:
# Residual block from Putzky and Welling (2019): https://arxiv.org/abs/1911.10914
# Author: Philipp Witte, pwitte3@gatech.edu
# Date: January 2020

export ResidualBlock

"""
    RB = ResidualBlock(nx, ny, n_in, n_hidden, batchsize; k1=3, k2=3, p1=1, p2=1, s1=1, s2=1, fan=false) (2D)

    RB = ResidualBlock(nx, ny, nz, n_in, n_hidden, batchsize; k1=3, k2=3, p1=1, p2=1, s1=1, s2=1, fan=false) (3D)

or

    RB = ResidualBlock(nx, ny, n_in, n_hidden, batchsize; k1=3, k2=3, p1=1, p2=1, s1=1, s2=1, fan=false) (2D)

    RB = ResidualBlock(nx, ny, nz, n_in, n_hidden, batchsize; k1=3, k2=3, p1=1, p2=1, s1=1, s2=1, fan=false) (3D)

 Create a (non-invertible) residual block, consisting of three convolutional layers and activation functions.
 The first convolution is a downsampling operation with a stride equal to the kernel dimension. The last
 convolution is the corresponding transpose operation and upsamples the data to either its original dimensions
 or to twice the number of input channels (for `fan=true`). The first and second layer contain a bias term.

 *Input*:

 - `nx`, `ny`, `nz`: spatial dimensions of input

 - `n_in`, `n_hidden`: number of input and hidden channels

 - `k1`, `k2`: kernel size of convolutions in residual block. `k1` is the kernel of the first and third
    operator, `k2` is the kernel size of the second operator.

 - `p1`, `p2`: padding for the first and third convolution (`p1`) and the second convolution (`p2`)

 - `s1`, `s2`: stride for the first and third convolution (`s1`) and the second convolution (`s2`)

 - `fan`: bool to indicate whether the ouput has twice the number of input channels. For `fan=false`, the last
    activation function is a gated linear unit (thereby bringing the output back to the original dimensions).
    For `fan=true`, the last activation is a ReLU, in which case the output has twice the number of channels
    as the input.

or

 - `W1`, `W2`, `W3`: 4D tensors of convolutional weights

 - `b1`, `b2`: bias terms

 - `nx`, `ny`: spatial dimensions of input image

 *Output*:

 - `RB`: residual block layer

 *Usage:*

 - Forward mode: `Y = RB.forward(X)`

 - Backward mode: `ΔX = RB.backward(ΔY, X)`

 *Trainable parameters:*

 - Convolutional kernel weights `RB.W1`, `RB.W2` and `RB.W3`

 - Bias terms `RB.b1` and `RB.b2`

 See also: [`get_params`](@ref), [`clear_grad!`](@ref)
"""
struct ResidualBlock <: NeuralNetLayer
    W1::Parameter
    W2::Parameter
    W3::Parameter
    b1::Parameter
    b2::Parameter
    fan::Bool
    cdims1::DenseConvDims
    cdims2::DenseConvDims
    cdims3::DenseConvDims
end

@Flux.functor ResidualBlock

#######################################################################################################################
# 2D Constructors

# Constructor 2D
function ResidualBlock(nx, ny, n_in, n_hidden, batchsize; k1=3, k2=3, p1=1, p2=1, s1=1, s2=1, fan=false)

    # Initialize weights
    W1 = Parameter(glorot_uniform(k1, k1, n_in, n_hidden))
    W2 = Parameter(glorot_uniform(k2, k2, n_hidden, n_hidden))
    W3 = Parameter(1f-2*glorot_uniform(k1, k1, 2*n_in, n_hidden))
    b1 = Parameter(zeros(Float32, n_hidden))
    b2 = Parameter(zeros(Float32, n_hidden))

    # Dimensions for convolutions
    cdims1 = DenseConvDims((nx, ny, n_in, batchsize), (k1, k1, n_in, n_hidden);
        stride=(s1, s1), padding=(p1, p1))
    cdims2 = DenseConvDims((Int(nx/s1), Int(ny/s1), n_hidden, batchsize),
        (k2, k2, n_hidden, n_hidden); stride=(s2, s2), padding=(p2, p2))
    cdims3 = DenseConvDims((nx, ny, 2*n_in, batchsize), (k1, k1, 2*n_in, n_hidden);
        stride=(s1, s1), padding=(p1 ,p1))

    return ResidualBlock(W1, W2, W3, b1, b2, fan, cdims1, cdims2, cdims3)
end

# Constructor for given weights 2D
function ResidualBlock(W1, W2, W3, b1, b2, nx, ny, batchsize; p1=1, p2=1, s1=1, s2=1, fan=false)

    # Make weights parameters
    W1 = Parameter(W1)
    W2 = Parameter(W2)
    W3 = Parameter(W3)
    b1 = Parameter(b1)
    b2 = Parameter(b2)

    # Dimensions for convolutions
    k1, n_in, n_hidden = size(W1)[2:4]
    k2 = size(W2)[1]
    cdims1 = DenseConvDims((nx, ny, n_in, batchsize), (k1, k1, n_in, n_hidden);
        stride=(s1, s1), padding=(p1, p1))
    cdims2 = DenseConvDims((Int(nx/s1), Int(ny/s1), n_hidden, batchsize),
        (k2, k2, n_hidden, n_hidden); stride=(s2, s2), padding=(p2, p2))
    cdims3 = DenseConvDims((nx, ny, 2*n_in, batchsize), (k1, k1, 2*n_in, n_hidden);
        stride=(s1, s1), padding=(p1, p1))

    return ResidualBlock(W1, W2, W3, b1, b2, fan, cdims1, cdims2, cdims3)
end

#######################################################################################################################
# 3D Constructors

# Constructor 3D
function ResidualBlock(nx, ny, nz, n_in, n_hidden, batchsize; k1=3, k2=3, p1=1, p2=1, s1=1, s2=1, fan=false)

    # Initialize weights
    W1 = Parameter(glorot_uniform(k1, k1, k1, n_in, n_hidden))
    W2 = Parameter(glorot_uniform(k2, k2, k2, n_hidden, n_hidden))
    W3 = Parameter(1f-2*glorot_uniform(k1, k1, k1, 2*n_in, n_hidden))
    b1 = Parameter(zeros(Float32, n_hidden))
    b2 = Parameter(zeros(Float32, n_hidden))

    # Dimensions for convolutions
    cdims1 = DenseConvDims((nx, ny, nz, n_in, batchsize), (k1, k1, k1, n_in, n_hidden);
        stride=(s1, s1, s1), padding=(p1, p1, p1))
    cdims2 = DenseConvDims((Int(nx/s1), Int(ny/s1), Int(nz/s1), n_hidden, batchsize),
        (k2, k2, k2, n_hidden, n_hidden); stride=(s2, s2, s2), padding=(p2, p2, p2))
    cdims3 = DenseConvDims((nx, ny, nz, 2*n_in, batchsize), (k1, k1, k1, 2*n_in, n_hidden);
        stride=(s1, s1, s1), padding=(p1, p1, p1))

    return ResidualBlock(W1, W2, W3, b1, b2, fan, cdims1, cdims2, cdims3)
end

# Constructor for given weights 3D
function ResidualBlock(W1, W2, W3, b1, b2, nx::Int64, ny::Int64, nz::Int64, batchsize::Int64;
                      p1=1, p2=1, s1=1, s2=1, fan=false)

    # Make weights parameters
    W1 = Parameter(W1)
    W2 = Parameter(W2)
    W3 = Parameter(W3)
    b1 = Parameter(b1)
    b2 = Parameter(b2)

    # Dimensions for convolutions
    k1, n_in, n_hidden = size(W1)[3:5]
    k2 = size(W2)[1]
    cdims1 = DenseConvDims((nx, ny, nz, n_in, batchsize), (k1, k1, k1, n_in, n_hidden);
        stride=(s1, s1, s1), padding=(p1, p1, p1))
    cdims2 = DenseConvDims((Int(nx/s1), Int(ny/s1), Int(nz/s1), n_hidden, batchsize),
        (k2, k2, k2, n_hidden, n_hidden); stride=(s2, s2, s2), padding=(p2, p2, p2))
    cdims3 = DenseConvDims((nx, ny, nz, 2*n_in, batchsize), (k1, k1, k1, 2*n_in, n_hidden);
        stride=(s1, s1, s1), padding=(p1, p1, p1))

    return ResidualBlock(W1, W2, W3, b1, b2, fan, cdims1, cdims2, cdims3)
end

#######################################################################################################################
# Functions

# Forward 2D
function forward(X1::AbstractArray{Float32, 4}, RB::ResidualBlock; save=false)

    Y1 = conv(X1, RB.W1.data, RB.cdims1) .+ reshape(RB.b1.data, 1, 1, :, 1)
    X2 = ReLU(Y1)

    Y2 = X2 + conv(X2, RB.W2.data, RB.cdims2) .+ reshape(RB.b2.data, 1, 1, :, 1)
    X3 = ReLU(Y2)

    Y3 = ∇conv_data(X3, RB.W3.data, RB.cdims3)
    RB.fan == true ? (X4 = ReLU(Y3)) : (X4 = GaLU(Y3))

    if save == false
        return X4
    else
        return Y1, Y2, Y3, X2, X3
    end
end

# Forward 3D
function forward(X1::AbstractArray{Float32, 5}, RB::ResidualBlock; save=false)

    Y1 = conv(X1, RB.W1.data, RB.cdims1) .+ reshape(RB.b1.data, 1, 1, 1, :, 1)
    X2 = ReLU(Y1)

    Y2 = X2 + conv(X2, RB.W2.data, RB.cdims2) .+ reshape(RB.b2.data, 1, 1, 1, :, 1)
    X3 = ReLU(Y2)

    Y3 = ∇conv_data(X3, RB.W3.data, RB.cdims3)
    RB.fan == true ? (X4 = ReLU(Y3)) : (X4 = GaLU(Y3))

    if save == false
        return X4
    else
        return Y1, Y2, Y3, X2, X3
    end
end

# Backward 2D
function backward(ΔX4::AbstractArray{Float32, 4}, X1::AbstractArray{Float32, 4}, RB::ResidualBlock)

    # Recompute forward states from input X
    Y1, Y2, Y3, X2, X3 = forward(X1, RB; save=true)

    # Backpropagate residual ΔX4 and compute gradients
    RB.fan == true ? (ΔY3 = ReLUgrad(ΔX4, Y3)) : (ΔY3 = GaLUgrad(ΔX4, Y3))
    ΔX3 = conv(ΔY3, RB.W3.data, RB.cdims3)
    ΔW3 = ∇conv_filter(ΔY3, X3, RB.cdims3)

    ΔY2 = ReLUgrad(ΔX3, Y2)
    ΔX2 = ∇conv_data(ΔY2, RB.W2.data, RB.cdims2) + ΔY2
    ΔW2 = ∇conv_filter(X2, ΔY2, RB.cdims2)
    Δb2 = sum(ΔY2, dims=[1,2,4])[1,1,:,1]

    ΔY1 = ReLUgrad(ΔX2, Y1)
    ΔX1 = ∇conv_data(ΔY1, RB.W1.data, RB.cdims1)
    ΔW1 = ∇conv_filter(X1, ΔY1, RB.cdims1)
    Δb1 = sum(ΔY1, dims=[1,2,4])[1,1,:,1]

    # Set gradients
    RB.W1.grad = ΔW1
    RB.W2.grad = ΔW2
    RB.W3.grad = ΔW3
    RB.b1.grad = Δb1
    RB.b2.grad = Δb2

    return ΔX1
end

# Backward 3D
function backward(ΔX4::AbstractArray{Float32, 5}, X1::AbstractArray{Float32, 5}, RB::ResidualBlock)

    # Recompute forward states from input X
    Y1, Y2, Y3, X2, X3 = forward(X1, RB; save=true)

    # Backpropagate residual ΔX4 and compute gradients
    RB.fan == true ? (ΔY3 = ReLUgrad(ΔX4, Y3)) : (ΔY3 = GaLUgrad(ΔX4, Y3))
    ΔX3 = conv(ΔY3, RB.W3.data, RB.cdims3)
    ΔW3 = ∇conv_filter(ΔY3, X3, RB.cdims3)

    ΔY2 = ReLUgrad(ΔX3, Y2)
    ΔX2 = ∇conv_data(ΔY2, RB.W2.data, RB.cdims2) + ΔY2
    ΔW2 = ∇conv_filter(X2, ΔY2, RB.cdims2)
    Δb2 = sum(ΔY2, dims=[1,2,3,5])[1,1,1,:,1]

    ΔY1 = ReLUgrad(ΔX2, Y1)
    ΔX1 = ∇conv_data(ΔY1, RB.W1.data, RB.cdims1)
    ΔW1 = ∇conv_filter(X1, ΔY1, RB.cdims1)
    Δb1 = sum(ΔY1, dims=[1,2,3,5])[1,1,1,:,1]

    # Set gradients
    RB.W1.grad = ΔW1
    RB.W2.grad = ΔW2
    RB.W3.grad = ΔW3
    RB.b1.grad = Δb1
    RB.b2.grad = Δb2

    return ΔX1
end

# Clear gradients
function clear_grad!(RB::ResidualBlock)
    RB.W1.grad = nothing
    RB.W2.grad = nothing
    RB.W3.grad = nothing
    RB.b1.grad = nothing
    RB.b2.grad = nothing
end

"""
    P = get_params(NL::NeuralNetLayer)

 Returns a cell array of all parameters in the network layer. Each cell
 entry contains a reference to the original parameter; i.e. modifying
 the paramters in `P`, modifies the parameters in `NL`.
"""
get_params(RB::ResidualBlock) = [RB.W1, RB.W2, RB.W3, RB.b1, RB.b2]

"""
    P = put_params!(NL::NeuralNetLayer, Params::Array{Parameter,1})

 Inputs a initialized ResidualBlock and desired parameters values to load.
"""
function put_params!(RB::ResidualBlock, Params::Array{Any,1})
    RB.W1.data = Params[1].data
    RB.W2.data = Params[2].data
    RB.W3.data = Params[3].data
    RB.b1.data = Params[4].data
    RB.b2.data = Params[5].data
end


text:
using PICAnalysisTools
using Test

@testset "PICAnalysisTools.jl" begin
    # Write your tests here.
end


text:
function NewFindMarker(r)
    counter = 1
    while !ismarker(r)
        for side in 2:5
            movementsbyfixmovesandcheckmarker(r, HorizonSide(side % 4), counter)
            if ismarker(r)
                break
            end
            if side == 3 || side == 5
                counter += 1
            end
        end
    end
end

function movementsbyfixmovesandcheckmarker(r, side, moves)
    while moves != 0
        trymove!(r, side)
        if ismarker(r)
            break
        end
        moves -= 1
    end
end

function trymove!(r::Robot, side::HorizonSide)
    if isborder(r,side)
        moves = 1
        end_side = HorizonSide((Int(side) + 1) % 4)
        while isborder(r, side)
            movementsbyfixmoves(r,HorizonSide((Int(side) + 1) % 4), moves)
            if !isborder(r, side)
                end_side = HorizonSide((Int(side) + 3) % 4)
                break
            end
            moves += 1
            movementsbyfixmoves(r,HorizonSide((Int(side) + 3) % 4), moves)
            moves += 1
        end
        move!(r,side)
        movementsbyfixmoves(r, end_side, Int(floor((moves+1)/2)))
    else
        move!(r,side)
    end
end

function movementsbyfixmoves(r, side, moves)
    while moves != 0
        move!(r, side)
        moves -= 1
    end
end


text:
# Autogenerated wrapper script for SDL2_net_jll for x86_64-unknown-freebsd
export libSDL2_net

using SDL2_jll
JLLWrappers.@generate_wrapper_header("SDL2_net")
JLLWrappers.@declare_library_product(libSDL2_net, "libSDL2_net-2.0.so.0")
function __init__()
    JLLWrappers.@generate_init_header(SDL2_jll)
    JLLWrappers.@init_library_product(
        libSDL2_net,
        "lib/libSDL2_net.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
function get_time_series(
    psi_container::PSIContainer,
    component::PSY.Component,
    forecast_label::String,
)
    initial_time = model_initial_time(psi_container)
    @debug initial_time
    use_forecast_data = model_uses_forecasts(psi_container)
    time_steps = model_time_steps(psi_container)
    if use_forecast_data
        forecast = PSY.get_time_series(
            PSY.Deterministic,
            component,
            forecast_label;
            start_time = initial_time,
            count = 1,
        )
        ts_vector = IS.get_time_series_values(
            component,
            forecast,
            initial_time;
            len = length(time_steps),
            ignore_scaling_factors = true,
        )
        return ts_vector
    else
        return ones(time_steps[end])
    end
end


text:
"""
    GEEBRA_results(results::Union{NLsolve.SolverResults, Optim.MultivariateOptimizationResults, Optim.UnivariateOptimizationResults}, theta::Vector, data::Any,  data::Any, template::Union{objective_function_template, estimating_function_template}, br::Bool, has_objective::Bool, br_method::String)

Composite type for the output of [`fit`](@ref) for an [`objective_function_template`](@ref) or an [`estimating_function_template`](@ref).
"""
struct GEEBRA_results
    results::Union{NLsolve.SolverResults, Optim.MultivariateOptimizationResults, Optim.UnivariateOptimizationResults}
    theta::Vector
    data::Any
    template::Union{objective_function_template, estimating_function_template}
    br::Bool
    has_objective::Bool
    br_method::String
end

"""
    vcov(results::GEEBRA_results)

Compute an esitmate of the variance-covariance matrix of the `M`-estimator or its reduced-bias version from the output of [`fit`](@ref) for an [`objective_function_template`](@ref) or an [`estimating_function_template`](@ref).
"""
function vcov(results::GEEBRA_results)
    if (results.has_objective)
        obj_quantities(results.theta, results.data, results.template, false)[1]
    else
        ef_quantities(results.theta, results.data, results.template, false)[1]
    end
end

"""
    tic(results::GEEBRA_results)

Compute the Takeuchi Information Criterion at the `M`-estimator or its reduced-bias version from the output of [`fit`](@ref) for an [`objective_function_template`](@ref). `nothing` is returned if `results` is the output of [`fit`](@ref) for an [`estimating_function_template`](@ref).
"""
function tic(results::GEEBRA_results)
    if (results.has_objective)
        obj = objective_function(results.theta, results.data, results.template, false)
        quants = obj_quantities(results.theta, results.data, results.template, true)
        -2 * (obj + 2 * quants[1])
    end
end

"""
    aic(results::GEEBRA_results)

Compute the Akaike Information Criterion at the `M`-estimator or its reduced-bias version from the output of [`fit`](@ref) for an [`objective_function_template`](@ref). `nothing` is returned if `results` is the output of [`fit`](@ref) for an [`estimating_function_template`](@ref).
"""
function aic(results::GEEBRA_results)
    if (results.has_objective)
        obj = objective_function(results.theta, results.data, results.template, false)
        p = length(results.theta)
        -2 * (obj - p)
    end
end

"""
    coef(results::GEEBRA_results)

Extract the `M`-estimates or their reduced-bias versions from the output of [`fit`](@ref) for an [`objective_function_template`](@ref) or an [`estimating_function_template`](@ref).
"""
function coef(results::GEEBRA_results)
    results.theta
end

"""
    show(io::IO, results::GEEBRA_results; digits::Real = 4)

`show` method for `GEEBRA_results` objects. If `GEEBRA_results.has_object == true`, then the result of `aic(results)` and `tic(results)` are also printed.
"""
function Base.show(io::IO, results::GEEBRA_results;
                   digits::Real = 4)
    theta = results.theta
    p = length(theta)
    v = vcov(results)
    if results.has_objective
        println(io,
                (results.br ? "RBM" : "M") * "-estimation with objective contributions ",
                results.template.obj_contribution)
    else
        println(io,
                (results.br ? "RBM" : "M") * "-estimation with estimating function contributions ",
                results.template.ef_contribution)
    end
    if (results.br)
        println(io, "Bias reduction method: ", results.br_method)
    end
    println(io)
    show(io, coeftable(results))
    if results.has_objective
        objfun = objective_function(results.theta, results.data, results.template, results.br)
        if results.br
            print(io, "\nMaximum penalized objetive:\t", round(objfun, digits = digits))
        else
            print(io, "\nMaximum objetive:\t\t", round(objfun, digits = digits))
        end
        print(io, "\nTakeuchi information criterion:\t", round(tic(results), digits = digits))
        print(io, "\nAkaike information criterion:\t", round(aic(results), digits = digits))
    else
        estfun = estimating_function(results.theta, results.data, results.template, results.br)
        if results.br
            print(io, "\nAdjusted estimating functions:\t", estfun)
        else
            print(io, "\nEstimating functions:\t", estfun)
        end
    end
end

"""
    stderror(results::GEEBRA_results)

Compute esitmated standard errors for the `M`-estimator or its reduced-bias version from the output of [`fit`](@ref) for an [`objective_function_template`](@ref) or an [`estimating_function_template`](@ref).
"""
function stderror(results::GEEBRA_results)
    sqrt.(diag(vcov(results)))
end

"""
    coeftable(results::GEEBRA_results; level::Real=0.95)

Return a `StatsBase.CoefTable` for the `M`-estimator or its reduced-bias version from the output of [`fit`](@ref) for an [`objective_function_template`](@ref) or an [`estimating_function_template`](@ref). `level` can be used to set the level of the reported Wald-type confidence intervals (using quantiles of the standard normal distribution).
"""
function coeftable(results::GEEBRA_results; level::Real=0.95)
    cc = coef(results)
    se = stderror(results)
    zz = cc ./ se
    p = 2 * ccdf.(Ref(Normal()), abs.(zz))
    ci = se * quantile(Normal(), (1-level)/2)
    levstr = isinteger(level*100) ? string(Integer(level*100)) : string(level*100)
    CoefTable(hcat(cc, se, zz, p, cc + ci, cc - ci),
              ["Estimate","Std. Error","z value","Pr(>|z|)","Lower $levstr%","Upper $levstr%"],
              ["theta[$i]" for i = 1:length(cc)], 4)
end


text:
module Telescope

using Telescope_jll
export Telescope_jll

using CEnum

struct TS_PositionInfo
    x::Cfloat
    y::Cfloat
    z::Cfloat
end

struct TS_VelocityInfo
    x::Cfloat
    y::Cfloat
    z::Cfloat
end

struct TS_CollisionEvent
    id1::Cint
    id2::Cint
    colliding::Bool
end

function TS_BtAddRigidBox(id, hx, hy, hz, m, px, py, pz, isKinematic)
    ccall((:TS_BtAddRigidBox, libtelescope), Cvoid, (Cint, Cfloat, Cfloat, Cfloat, Cfloat, Cfloat, Cfloat, Cfloat, Bool), id, hx, hy, hz, m, px, py, pz, isKinematic)
end

function TS_BtAddStaticBox(id, hx, hy, hz, px, py, pz)
    ccall((:TS_BtAddStaticBox, libtelescope), Cvoid, (Cint, Cfloat, Cfloat, Cfloat, Cfloat, Cfloat, Cfloat), id, hx, hy, hz, px, py, pz)
end

function TS_BtAddTriggerBox(id, hx, hy, hz, px, py, pz)
    ccall((:TS_BtAddTriggerBox, libtelescope), Cvoid, (Cint, Cfloat, Cfloat, Cfloat, Cfloat, Cfloat, Cfloat), id, hx, hy, hz, px, py, pz)
end

function TS_BtRemovePhysicsObject(id)
    ccall((:TS_BtRemovePhysicsObject, libtelescope), Cvoid, (Cint,), id)
end

function TS_BtSetLinearVelocity(id, vx, vy, vz)
    ccall((:TS_BtSetLinearVelocity, libtelescope), Cvoid, (Cint, Cfloat, Cfloat, Cfloat), id, vx, vy, vz)
end

function TS_BtGetLinearVelocity(id)
    ccall((:TS_BtGetLinearVelocity, libtelescope), TS_VelocityInfo, (Cint,), id)
end

function TS_BtSetGravity(gx, gy, gz)
    ccall((:TS_BtSetGravity, libtelescope), Cvoid, (Cfloat, Cfloat, Cfloat), gx, gy, gz)
end

function TS_BtSetCollisionMargin(id, margin)
    ccall((:TS_BtSetCollisionMargin, libtelescope), Cvoid, (Cint, Cfloat), id, margin)
end

# no prototype is found for this function at telescope.h:51:6, please use with caution
function TS_BtStepSimulation()
    ccall((:TS_BtStepSimulation, libtelescope), Cvoid, ())
end

# no prototype is found for this function at telescope.h:53:26, please use with caution
function TS_BtGetNextCollision()
    ccall((:TS_BtGetNextCollision, libtelescope), TS_CollisionEvent, ())
end

function TS_BtGetPosition(id)
    ccall((:TS_BtGetPosition, libtelescope), TS_PositionInfo, (Cint,), id)
end

# no prototype is found for this function at telescope.h:57:14, please use with caution
function TS_SDLGetError()
    ccall((:TS_SDLGetError, libtelescope), Ptr{Cchar}, ())
end

function TS_VkCmdDrawRect(r, g, b, a, x, y, w, h)
    ccall((:TS_VkCmdDrawRect, libtelescope), Cvoid, (Cfloat, Cfloat, Cfloat, Cfloat, Cfloat, Cfloat, Cfloat, Cfloat), r, g, b, a, x, y, w, h)
end

function TS_VkCmdDrawSprite(img, r, g, b, a, rx, ry, rw, rh, cw, ch, ci, cj, px, py, sx, sy)
    ccall((:TS_VkCmdDrawSprite, libtelescope), Cvoid, (Ptr{Cchar}, Cfloat, Cfloat, Cfloat, Cfloat, Cint, Cint, Cint, Cint, Cint, Cint, Cint, Cint, Cfloat, Cfloat, Cfloat, Cfloat), img, r, g, b, a, rx, ry, rw, rh, cw, ch, ci, cj, px, py, sx, sy)
end

function TS_VkCmdClearColorImage(r, g, b, a)
    ccall((:TS_VkCmdClearColorImage, libtelescope), Cvoid, (Cfloat, Cfloat, Cfloat, Cfloat), r, g, b, a)
end

# no prototype is found for this function at telescope.h:65:6, please use with caution
function TS_VkBeginDrawPass()
    ccall((:TS_VkBeginDrawPass, libtelescope), Cvoid, ())
end

function TS_VkEndDrawPass(r, g, b, a)
    ccall((:TS_VkEndDrawPass, libtelescope), Cvoid, (Cfloat, Cfloat, Cfloat, Cfloat), r, g, b, a)
end

function TS_Init(ttl, wdth, hght)
    ccall((:TS_Init, libtelescope), Cvoid, (Ptr{Cchar}, Cint, Cint), ttl, wdth, hght)
end

# no prototype is found for this function at telescope.h:71:6, please use with caution
function TS_Quit()
    ccall((:TS_Quit, libtelescope), Cvoid, ())
end

function TS_PlaySound(sound_file, loops, ticks)
    ccall((:TS_PlaySound, libtelescope), Cvoid, (Ptr{Cchar}, Cint, Cint), sound_file, loops, ticks)
end

# exports
const PREFIXES = ["TS_"]
for name in names(@__MODULE__; all=true), prefix in PREFIXES
    if startswith(string(name), prefix)
        @eval export $name
    end
end

end # module


text:
#!/usr/bin/env julia

Base.eval(:(have_color=true))

function die(msg)
    print_with_color(:red, STDERR, "ERROR: " * msg * "\n")
    exit(1)
end

if length(ARGS) != 1
    println("Usage: ", PROGRAM_FILE, "<folder>")
    println("    where <folder> is the output of a run in Janus.")
    println("    The program will output in stdout a csv with the")
    println("    temporal correlation of the configurations.")
    quit()
end

folder = ARGS[1]

if folder[end] == '/'
    thefolder = folder |> dirname |> basename
else
    thefolder = folder |> basename
end

if !isfile(folder * "/S00001/run.log")
    die("No run.log in folder.")
end

Jhash = @sprintf("%d",read(folder * "/S00001/conf_Js") |> hash)

cd(folder * "/S00001/confs")

confs = readdir(".")

isconf(f) = ismatch(r"conf_[0-9A-F]{10}", f)

if !all(map(isconf, confs))
    warn("Some files don't match conf_\\x{10}, filtering them.")
    filter!(isconf, confs)
end

function uncompact(load)
    N = length(load)
    data = BitArray(N*8)
    for ib in 1:N
        is = 8(ib-1)
        for j in 0:7
            data[1+is+j] = (load[ib]>>j) & 0x1
        end
    end
    return data
end

function conf2data(file)
    data = Dict()
    open(file,"r") do f
        data["itime"]           = read(f,Cint)
        data["Dt"]              = read(f,Clonglong)
        data["t"]               = read(f,Clonglong) # The MC step
        data["timewindow"]      = read(f,Clonglong)
        data["flagJ"]           = read(f,Cint)
        data["flagS"]           = read(f,Cint)
        data["flagh"]           = read(f,Cint)

        data["L"]               = read(f,Cint)
        data["β"]               = read(f,Cdouble)
        data["H"]               = read(f,Cdouble)
        data["tHon"]            = read(f,Culonglong)

        data["ir"]              = read(f,Cint)
        data["nfiles_dir"]      = read(f,Cint)

        data["energy_sigma"]    = read(f,Cdouble)
        data["energy_tau"]      = read(f,Cdouble)
        data["energymag_sigma"] = read(f,Cdouble)
        data["energymag_tau"]   = read(f,Cdouble)
        data["mag_sigma"]       = read(f,Cdouble)
        data["mag_tau"]         = read(f,Cdouble)

        num                     = data["L"]^3 ÷ 8
        data["spins_σ"]         = uncompact(read(f,Cchar,num))
        data["spins_τ"]         = uncompact(read(f,Cchar,num))
        @assert eof(f)
    end
    return data
end

function getC(bconfA, bconfB)
    spinsA_σ = bconfA["spins_σ"]
    spinsA_τ = bconfA["spins_τ"]
    spinsB_σ = bconfB["spins_σ"]
    spinsB_τ = bconfB["spins_τ"]
    N = length(spinsA_σ)
    differing_spins_σ = map(≠,spinsA_σ, spinsB_σ) |> sum
    differing_spins_τ = map(≠,spinsA_τ, spinsB_τ) |> sum
    # C = 1*N↑↑ + (-1)*N↑↓
    Cσ = (N-differing_spins_σ) - differing_spins_σ
    Cσ /= N
    Cτ = (N-differing_spins_τ) - differing_spins_τ
    Cτ /= N
    return (Cσ+Cτ)/2
end

function gettimes()

    function rectify_times!(array)
        L = length(array)
        for i in 2:L
            if array[i] == array[i-1]
                for j in i:L
                    array[j] += 1
                end
            end
        end
    end

    const tw1list = map(Int64,[1e4,1e6,1e8])
    tws = logspace(log10(1),log10(3e8+10),100)
    tws = map(x->floor(Int64,x), tws)
    tws = vcat(tws, tw1list) |> sort
    rectify_times!(tws)

    tws = unique(tws) |> sort
    ts = logspace(4,8,10)
    ts = map(x->floor(Int64,x), ts) |> sort

    return tws, ts
end

function measureandprint(confs, Jhash)
    confs = sort(confs)

    twlist, tlist = gettimes()
    ℓ   = length(confs)
    bconfs = map(conf2data, confs)

    L = bconfs[1]["L"]
    Ts = [1/c["β"] for c in bconfs]
    if Ts |> unique |> length != 1
        die("Multiple temperatures in quench!")
    end
    T = Ts[1]

    timeslist = [c["t"] for c in bconfs]
    realmag_σ_list = [c["mag_sigma"] for c in bconfs]
    realmag_τ_list = [c["mag_tau"] for c in bconfs]
    E_σ_list = [c["energy_sigma"] for c in bconfs]
    E_τ_list = [c["energy_tau"] for c in bconfs]

    E_list = 0.5*(E_σ_list + E_τ_list)
    M_list = 0.5*(realmag_σ_list + realmag_τ_list)

    # Consistency check
    measmag_σ_list = [(2sum(c["spins_σ"])-L^3)/L^3 for c in bconfs]
    measmag_τ_list = [(2sum(c["spins_τ"])-L^3)/L^3 for c in bconfs]
    if measmag_σ_list ≉ realmag_σ_list
        die("Readed magnetization is not ≈ calculated magnetization (σ spins)")
    end
    if measmag_τ_list ≉ realmag_τ_list
        die("Readed magnetization is not ≈ calculated magnetization (τ spins)")
    end

    for i in 1:ℓ
        tw = bconfs[i]["t"]
        tw ∉ twlist && continue
        for j in i+1:ℓ
            tw_plus_t = bconfs[j]["t"]
            t = tw_plus_t - tw
            t ∉ tlist && continue
            # Wohoo! measure!
            T = 1/bconfs[i]["β"]
            C = getC(bconfs[i],bconfs[j])
            println("$L,$tw,$T,\"$Jhash\",\"Energy\",0,$(E_list[i])")
            println("$L,$tw,$T,\"$Jhash\",\"Magnetization\",0,$(M_list[i])")
            println("$L,$tw,$T,\"$Jhash\",\"Correlation\",$t,$C")
        end
    end

end

println("\"L\",\"mc\",\"T\",\"Jhash\",\"observable\",\"parameter\",\"value\"")
measureandprint(confs, Jhash)

max_stars_repo_name: RedPointyJackson/tfg

max_issues_repo_name: RedPointyJackson/tfg

max_forks_repo_name: RedPointyJackson/tfg


text:
function Base.iterate(::Type{Mod{m}}) where {m}
    return Mod{m}(0), 0
end

function Base.iterate(::Type{Mod{m}}, s) where {m}
    if s == m - 1
        return nothing
    end
    s += 1
    return Mod{m}(s), s
end

Base.IteratorSize(::Type{Mod{m}}) where {m} = Base.HasLength()

Base.length(::Type{Mod{m}}) where {m} = m

function Base.iterate(::Type{GaussMod{m}}) where {m}
    return GaussMod{m}(0), 0
end

function Base.iterate(::Type{GaussMod{m}}, s) where {m}
    if s == m * m - 1
        return nothing
    end
    s += 1
    a, b = divrem(s, m)
    return GaussMod{m}(a + b * im), s
end

Base.IteratorSize(::Type{GaussMod{m}}) where {m} = Base.HasLength()

Base.length(::Type{GaussMod{m}}) where {m} = m * m


text:
# Taken from https://github.com/cmcaine/advent2020/blob/main/scripts/get-data.jl
aoc_cookie = try
    read(joinpath(@__DIR__, "cookie"), String)
catch
    println("No cookie provided for the AoC website. Save your cookie in a file called 'cookie' in the /src/EventUtils folder")
    nothing
end

function get_data(day)
    response = HTTP.get("https://adventofcode.com/2020/day/$day/input", cookies = Dict("session" => aoc_cookie))
    response.status != 200 && error(response)
    return response.body
end

function download_input(day)
    fname = joinpath("scripts", "Day$day", "input.txt")
    if !isfile(fname)
        data = get_data(day)
        write(fname, data)
    end
    return nothing
end


text:
#Load the required Packages
using JLD2, Random, LinearAlgebra, Statistics, CSV, DataFrames, FreqTables, Distributions
#Set the seed
Random.seed!(1234)
#(a)Create the following four matrices of random numbers
#(i) A 10 by 7 matrix drawn randomly from a uniform distribution U[-5,10]
A = rand(Uniform(-5,10),10,7)
#(ii) A 10 by 7 matrix drawn randomly from a normal distribution N[-2,15]
B = rand(Normal(-2,15),10,7)


text:
# #-----------------------------------------------------------------------------# XPixel
# """
# A pixel split into four values.  Adjacent pixels can be joined in multiple ways to create
# different shapes.

# ```
# (x, y + h)                       (x + w, y + h)
#             █████████████████
#             █ █           █ █
#             █   █   T   █   █
#             █     █   █     █
#             █  L    █    R  █
#             █     █   █     █
#             █   █   B   █   █
#             █ █           █ █
#             █████████████████
# (x, y)                           (x + w, y)
# ```
# """
# mutable struct XPixel{T, S}
#     x::S
#     y::S
#     h::S
#     w::S

#     top::T
#     right::T
#     bottom::T
#     left::T
# end
# Base.sum(x::XPixel) = x.top + x.right + x.bottom + x.left

# #-----------------------------------------------------------------------------# HexGrid
# struct HexGrid{EX, EY}
#     xedges::EX
#     yedges::EY
#     counts::Matrix{XPixel}
#     out::Int
# end