text:
@testset "io" begin
    fname = "testdata/n3-window-1m-LXe-20mm-1-5.h5"
    pdf   = read_abc(fname)

    @test length(pdf.volume_names ) ==    88
    @test length(pdf.process_names) ==     3
    @test nrow(pdf.sensor_xyz     ) == 47428
    @test ncol(pdf.sensor_xyz     ) ==     4
    @test nrow(pdf.primaries      ) ==     5
    @test ncol(pdf.primaries      ) ==     7
    @test nrow(pdf.vertices       ) ==     5
    @test ncol(pdf.vertices       ) ==    13
    @test nrow(pdf.total_charge   ) ==  1080
    @test ncol(pdf.total_charge   ) ==     3
    @test nrow(pdf.waveform       ) ==  8098
    @test ncol(pdf.waveform       ) ==     3

    evt_range = in(4995000:4995004)
    @test all(evt_range.(pdf.primaries.event_id))
    @test all(evt_range.(pdf.vertices.event_id))
    @test all(evt_range.(pdf.waveform.event_id))

    @test sum(pdf.total_charge.charge) == 8098

    outfile = "testdata/evtPar_example.h5"
    norm, rmin, rmax, evtdf = read_evtpar([outfile])

    @test norm == 5000
    @test isapprox(rmin, 351.5)
    @test isapprox(rmax, 372.9)
    @test nrow(evtdf) == 538
    exp_fields = fieldnames(ATools.EventParameters)
    @test ncol(evtdf) == length(exp_fields)
    @test all(in(exp_fields).(propertynames(evtdf)))
end

max_stars_repo_name: andLaing/ATools

max_issues_repo_name: andLaing/ATools

max_forks_repo_name: andLaing/ATools


text:
using CSV

@time a = CSV.read(
           "C:/data/Performance_All/Performance_2010Q3.txt",
               delim = '|',
           header = false,
           typemap = Dict(
               Int => Int32,
               Float64 => Float32,
               Union{Missing, Int} => Union{Missing, Int32},
               Union{Missing, Float64} => Union{Missing, Float32}
               ),
           copycols=true
       );

ext: jl


text:
module TestMOIWrapper

using Test
using YasolSolver
using JuMP
using MathOptInterface
const MOI = MathOptInterface

function runtests()
    for name in names(@__MODULE__; all = true)
        if !startswith("$(name)", "test_")
            continue
        end
        @testset "$(name)" begin
            getfield(@__MODULE__, name)()
        end
    end
end

function test_rawOptimizerAttributes()
    model = Model(() -> YasolSolver.Optimizer("C:/Yasol/Yasol_CLP"))

    set_optimizer_attribute(model, "time limit", 60)
    set_optimizer_attribute(model, "output info", 1)
    set_optimizer_attribute(model, "problem file name", "FullTest.qlp")

    @test get_optimizer_attribute(model, "time limit") == 60
    @test get_optimizer_attribute(model, "output info") == 1
    @test get_optimizer_attribute(model, "problem file name") == "FullTest.qlp"
end

function test_solverName()
    model = Model(() -> YasolSolver.Optimizer("C:/Yasol/Yasol_CLP"))
    @test MOI.get(model, MOI.SolverName()) == "YASOL"
end

function test_numberOfVariables()
    model = Model(() -> YasolSolver.Optimizer("C:/Yasol/Yasol_CLP"))
    @test MOI.get(model, MOI.NumberOfVariables()) == 0
    @variable(model, x1, integer=true, lower_bound=0, upper_bound=2)
    @test MOI.get(model, MOI.NumberOfVariables()) == 1
    @variable(model, x2, binary=true, lower_bound=0, upper_bound=1, YasolVariable, quantifier="all", block=2)
    @test MOI.get(model, MOI.NumberOfVariables()) == 2
end

function test_OptSense()
    model = Model(() -> YasolSolver.Optimizer("C:/Yasol/Yasol_CLP"))
    @variable(model, x1, integer=true, lower_bound=0, upper_bound=2)
    @variable(model, x2, binary=true, lower_bound=0, upper_bound=1, YasolVariable, quantifier="all", block=2)
    @constraint(model, con1, x1 + x2 <= 10)
    @objective(model, Max, x1+2x2)
    @test MOI.get(model, MOI.ObjectiveSense()) == MAX_SENSE
end

function test_OptFunctionType()
    model = Model(() -> YasolSolver.Optimizer("C:/Yasol/Yasol_CLP"))
    @variable(model, x1, integer=true, lower_bound=0, upper_bound=2)
    @variable(model, x2, binary=true, lower_bound=0, upper_bound=1, YasolVariable, quantifier="all", block=2)
    @constraint(model, con1, x1 + x2 <= 10)
    @objective(model, Max, x1+2x2)
    @test MOI.get(model, MOI.ObjectiveFunctionType()) == MathOptInterface.ScalarAffineFunction{Float64}
end

function test_variableAttributes()
    model = Model(() -> YasolSolver.Optimizer("C:/Yasol/Yasol_CLP"))
    @variable(model, x1, binary=true, lower_bound=0, upper_bound=1, YasolVariable, quantifier="all", block=2)
    @test MOI.get(model, YasolSolver.VariableAttribute("quantifier"), x1) == "all"
    @test MOI.get(model, YasolSolver.VariableAttribute("block"), x1) == 2
    @variable(model, x2, binary=true, lower_bound=0, upper_bound=1, YasolVariable, quantifier="exists", block=1)
    @test MOI.get(model, YasolSolver.VariableAttribute("quantifier"), x2) == "exists"
    @test MOI.get(model, YasolSolver.VariableAttribute("block"), x2) == 1
end

function test_constraintAttributes()
    model = Model(() -> YasolSolver.Optimizer("C:/Yasol/Yasol_CLP"))
    @variable(model, x1, binary=true, lower_bound=0, upper_bound=1, YasolVariable, quantifier="all", block=2)
    @variable(model, x2, binary=true, lower_bound=0, upper_bound=1, YasolVariable, quantifier="exists", block=1)
    @constraint(model, con1, -1*x1 +2*x2<=10, YasolConstraint, quantifier="all")
    @constraint(model, con2, -1*x1 +2*x2<=20, YasolConstraint, quantifier="exists")
    @test MOI.get(model, YasolSolver.ConstraintAttribute("quantifier"), con1) == "all"
    @test MOI.get(model, YasolSolver.ConstraintAttribute("quantifier"), con2) == "exists"
end

function test_status()
    model = Model(() -> YasolSolver.Optimizer("C:/Yasol/Yasol_CLP"))
    @test MOI.get(model, MOI.TerminationStatus()) == OPTIMIZE_NOT_CALLED
end

function test_status()
    model = Model(() -> YasolSolver.Optimizer("C:/Yasol/Yasol_CLP"))
    @test MOI.get(model, MOI.TerminationStatus()) == OPTIMIZE_NOT_CALLED
end

end

TestMOIWrapper.runtests()

"""
For local testing

# change yasol initial parameter
yasol.setInitialParameter("C:/Yasol", "writeOutputFile", 1)
# get initial parameter
@show yasol.getInitialParameter("C:/Yasol")

MAXIMIZE
1x1 +1x2 +1x3
SUBJECT TO
-2x2 -1x3 <= -2
-1x1 +2x2 +1x3 <= 2
2x1 + 4x2 <= 6
BOUNDS
0 <= x1 <= 2
0 <= x2 <= 1
0 <= x3 <= 2
GENERAL
x1
BINARY
x2
EXISTS
x1 x3
ALL
x2
ORDER
x1 x2 x3
END

# define model
cd("C:/Yasol")

model = Model(() -> yasol.Optimizer("C:/Yasol/Yasol_CLP"))

set_optimizer_attribute(model, "time limit", 60)
set_optimizer_attribute(model, "output info", 1)
set_optimizer_attribute(model, "problem file name", "Test08122021.qlp")

@variable(model, x1, integer=true, lower_bound=0, upper_bound=2)
MOI.set(model, yasol.VariableAttribute("quantifier"), x1, "exists")
MOI.set(model, yasol.VariableAttribute("block"), x1, 1)

@variable(model, x2, binary=true, lower_bound=0, upper_bound=1)
MOI.set(model, yasol.VariableAttribute("quantifier"), x2, "all")
MOI.set(model, yasol.VariableAttribute("block"), x2, 2)

@variable(model, x3, lower_bound=0, upper_bound=2)
MOI.set(model, yasol.VariableAttribute("quantifier"), x3, "exists")
MOI.set(model, yasol.VariableAttribute("block"), x3, 3)

@constraint(model, con1, -2*x2 -1x3 <= -2)

@constraint(model, con2, -1*x1 +2*x2 +1*x3 <= 2)

@constraint(model, con3, 2x1 + 4x2 <= 6)

@objective(model, Max, 1*x1 +1*x2 +1*x3)

optimize!(model)

# import solution
solution = yasol.importSolution("C:/Yasol/Test08122021.qlp.sol")
@show solution
yasol.printSolution(solution)

@show termination_status(model)

@show value(x1)

@show objective_value(model)

@show solve_time(model)

# access model parameters
print(model.moi_backend.optimizer.model)
print(model.moi_backend.optimizer.model.optimizer.is_objective_set)
print(model.moi_backend.model_cache.optattr.keys)
print(model.moi_backend.model)

MOI.get(model, yasol.VariableAttribute("quantifier"), x1)
MOI.get(model, yasol.VariableAttribute("block"), x1)
MOI.get(model, yasol.ConstraintAttribute("quantifier"), con3)

get_optimizer_attribute(model, "time limit")
get_optimizer_attribute(model, "output info")

# using constraints that have attributes

MINIMIZE
-x1 -2x2 +2x3 +x4
SUBJECT TO
E_Constraint1: x1 -2x2 +x3 -x4 <= 1
E_Constraint2: x1 +x2 +x3 -x4 <= 2
U_Constraint1: x1 +x2 +x3 <= 2
BOUNDS
0 <= x1 <= 1
0 <= x2 <= 1
0 <= x3 <= 1
0 <= x4 <= 1
BINARIES
x1 x2 x3 x4
EXISTS
x1 x2 x4
ALL
x3
ORDER
x1 x2 x3 x4
END

# define model
cd("C:/Yasol")

model = Model(() -> yasol.Optimizer("C:/Yasol/Yasol_CLP"))

set_optimizer_attribute(model, "time limit", 60)
set_optimizer_attribute(model, "output info", 1)
set_optimizer_attribute(model, "problem file name", "ConstraintExt.qlp")

@variable(model, x1, binary=true, lower_bound=0, upper_bound=1)
MOI.set(model, yasol.VariableAttribute("quantifier"), x1, "exists")
MOI.set(model, yasol.VariableAttribute("block"), x1, 1)

@variable(model, x2, binary=true, lower_bound=0, upper_bound=1)
MOI.set(model, yasol.VariableAttribute("quantifier"), x2, "exists")
MOI.set(model, yasol.VariableAttribute("block"), x2, 2)

@variable(model, x3, binary=true, lower_bound=0, upper_bound=1)
MOI.set(model, yasol.VariableAttribute("quantifier"), x3, "all")
MOI.set(model, yasol.VariableAttribute("block"), x3, 3)

@variable(model, x4, binary=true, lower_bound=0, upper_bound=1)
MOI.set(model, yasol.VariableAttribute("quantifier"), x4, "exists")
MOI.set(model, yasol.VariableAttribute("block"), x4, 4)

@constraint(model, con1, 1*x1 -2*x2 +1*x3 -1*x4 <= 1)
MOI.set(model, yasol.ConstraintAttribute("quantifier"), con1, "exists")

@constraint(model, con2, 1*x1 + 1*x2 +1*x3 -1*x4 <= 2)
MOI.set(model, yasol.ConstraintAttribute("quantifier"), con2, "exists")

@constraint(model, con3, 1*x1 + 1*x2 +1*x3 <= 2)
MOI.set(model, yasol.ConstraintAttribute("quantifier"), con3, "all")

@objective(model, Min, -1*x1 -2*x2 +2*x3 +1x4)

optimize!(model)

# import solution
solution = yasol.importSolution("C:/Yasol/ConstraintExt.qlp.sol")
@show solution
yasol.printSolution(solution)

#print(solver_name(model))

#@show model

#print(model)

#write_to_file(model, "model.mps")

# using JuMP extension
cd("C:/Yasol")

model = Model(() -> yasol.Optimizer("C:/Yasol/Yasol_CLP"))

set_optimizer_attribute(model, "time limit", 60)
set_optimizer_attribute(model, "output info", 1)
set_optimizer_attribute(model, "problem file name", "Test2-23_11.qlp")

@variable(model, x1, integer=true, lower_bound=0, upper_bound=2, YasolVariable, quantifier="exists", block=1)

@variable(model, x2, binary=true, lower_bound=0, upper_bound=1, YasolVariable, quantifier="all", block=2)

@variable(model, x3, lower_bound=0, upper_bound=2, YasolVariable, quantifier="exists", block=3)

@constraint(model, con1, -2*x2 -1x3 + 8 <= -2 + 2*x2)

@constraint(model, con2, -1*x1 +2*x2 +1*x3 >= -6 + 2*x2, YasolConstraint, quantifier="all")

@constraint(model, con3, 2x1 + 4x2 <= 6)

@objective(model, Max, 1*x1 +1*x2 +1*x3)

optimize!(model)

# import solution
solution = yasol.importSolution("C:/Yasol/Test25112021.qlp.sol")
@show solution
yasol.printSolution(solution)

@show termination_status(model)

@show value(x3)

@show objective_value(model)

@show solve_time(model)

# use full JuMP extension

MINIMIZE
-x1 -2x2 +2x3 +x4
SUBJECT TO
E_Constraint1: x1 -2x2 +x3 -x4 <= 1
E_Constraint2: x1 +x2 +x3 -x4 <= 2
U_Constraint1: x1 +x2 +x3 <= 2
BOUNDS
0 <= x1 <= 1
0 <= x2 <= 1
0 <= x3 <= 1
0 <= x4 <= 1
BINARIES
x1 x2 x3 x4
EXISTS
x1 x2 x4
ALL
x3
ORDER
x1 x2 x3 x4
END

# define model
cd("C:/Yasol")

model = Model(() -> yasol.Optimizer("C:/Yasol/Yasol_CLP"))

set_optimizer_attribute(model, "time limit", 60)
set_optimizer_attribute(model, "output info", 1)
set_optimizer_attribute(model, "problem file name", "ConstraintExt.qlp")

@variable(model, x1, binary=true, lower_bound=0, upper_bound=1, YasolVariable, quantifier="exists", block=1)

@variable(model, x2, binary=true, lower_bound=0, upper_bound=1, YasolVariable, quantifier="exists", block=2)

@variable(model, x3, binary=true, lower_bound=0, upper_bound=1, YasolVariable, quantifier="all", block=3)

@variable(model, x4, binary=true, lower_bound=0, upper_bound=1, YasolVariable, quantifier="exists", block=4)

@constraint(model, con1, 1*x1 -2*x2 +1*x3 -1*x4 +5 <= -1*x1 -3, YasolConstraint, quantifier="exists")

@constraint(model, con2, 1*x1 + 1*x2 +1*x3 -1*x4 -2 <= 2 - 2*x2, YasolConstraint, quantifier="exists")

@constraint(model, con3, 1*x1 + 1*x2 +1*x3 <= 1 + x3, YasolConstraint, quantifier="all")

@objective(model, Min, -1*x1 -2*x2 +2*x3 +1x4 - 5)

optimize!(model)

@show termination_status(model)

@show value(x1)

@show objective_value(model)

@show solve_time(model)
"""


text:
"""
The FADC supertype corresponds to the FADC component
    of the electronics chain in the DAQ setup.
"""
abstract type FADC end

"""
GenericFADC is a dummy FADC model that performs sampling
    of the input waveform based on the sampling interval <Δt>
"""
@with_kw struct GenericFADC <: FADC
    "Sampling interval of FADC, inverse of sampling rate"
    Δt::typeof(1.0*ns_unit) = uconvert(u"ns", inv(250u"MHz"))

    # to add: number of bits
end

"""
FADC with Flashcam algorithm (ToDo)
"""
struct Flashcam <: FADC end

"""
FADC with Struck algorithm (ToDo)
"""
struct Struck <: FADC end

"""
    GenericFADC(sim_conf)

PropDict -> GenericFADC

Construct a GenericFADC instance based on simulation
    configuration <sim_conf>
"""
function GenericFADC(sim_conf::PropDict)
    GenericFADC(
        Δt = haskey(sim_conf.setup.fadc, :sampling_rate) ? uconvert(u"ns", inv(T(sim_conf.setup.fadc.sampling_rate)u"MHz")) : T(sim_conf.setup.fadc.sampling_interval)u"ns",
    )
end

"""
    FADC(sim_conf)

PropDict -> <FADC>

Construct an FADC supertype instance based on simulation
    configuration <sim_conf>.
Type of <FADC> depends on the type specified in <sim_conf>
    (e.g. generic, Flashcam, Struck)
"""
function FADC(sim_conf::PropDict)
    if sim_conf.setup.fadc.type == "generic"
        GenericFADC(sim_conf)
    else
        @info "FADC type $(sim_conf.setup.fadc.type) not implemented!\n
        Available types: generic\n
        Planned in the future: Flashcam, Struck"
    end
end

# -------------------------------------------------------------------

"""
    simulate(wf, fadc)

RDWaveform, GenericFADC -> RDWaveform

Simulate effects of FADC module <fadc> on the waveform <wf>
    (sampling and digitization)
"""
function simulate(wf::RDWaveform, fadc::GenericFADC)
    # resample -> currently simply pick every Nth sample, more elaborate simulation later
    wf_sampled = wf.value[begin : Int(fadc.Δt/step(wf.time)) : end]
    t_sampled = range(T(0)u"ns", step = fadc.Δt, length = length(wf_sampled))
    wf_daq = RDWaveform(t_sampled, wf_sampled)

    # digitize
    wf_daq = RDWaveform(wf_daq.time, UInt16.(round.(wf_daq.value, digits = 0)))

    wf_daq
end


text:
using Plots
using LinearAlgebra

function plot_plane(center, spanning_vector_a, spanning_vector_b, color=:grey, ε = 1/2, frame=false, frame_color=:blue)
	x = [center[1] + ε * spanning_vector_a[1] + ε * spanning_vector_b[1],
	     center[1] - ε * spanning_vector_a[1] + ε * spanning_vector_b[1],
	     center[1] + ε * spanning_vector_a[1] - ε * spanning_vector_b[1],
	     center[1] - ε * spanning_vector_a[1] - ε * spanning_vector_b[1]]

	y = [center[2] + ε * spanning_vector_a[2] + ε * spanning_vector_b[2],
	     center[2] - ε * spanning_vector_a[2] + ε * spanning_vector_b[2],
	     center[2] + ε * spanning_vector_a[2] - ε * spanning_vector_b[2],
	     center[2] - ε * spanning_vector_a[2] - ε * spanning_vector_b[2]]

	z = [center[3] + ε * spanning_vector_a[3] + ε * spanning_vector_b[3],
	     center[3] - ε * spanning_vector_a[3] + ε * spanning_vector_b[3],
	     center[3] + ε * spanning_vector_a[3] - ε * spanning_vector_b[3],
	     center[3] - ε * spanning_vector_a[3] - ε * spanning_vector_b[3]]

	if frame
		path3d!(x,y,z;color=frame_color)
	end
	mesh3d!(x,y,z;color,connections=([0,3],[1,2],[3,0]))
end


text:
int main() {
    boolean b=false;
    if(b);
    else return 0;
}

max_stars_repo_name: carlostome/Javalette-Compiler

max_issues_repo_name: carlostome/Javalette-Compiler

max_forks_repo_name: carlostome/Javalette-Compiler


text:
module ChainRulesCore
using Base.Broadcast: broadcasted, Broadcasted, broadcastable, materialize, materialize!
using LinearAlgebra: LinearAlgebra
using MuladdMacro: @muladd

export on_new_rule, refresh_rules  # generation tools
export frule, rrule  # core function
export @non_differentiable, @scalar_rule, @thunk  # definition helper macros
export canonicalize, extern, unthunk  # differential operations
# differentials
export Composite, DoesNotExist, InplaceableThunk, One, Thunk, Zero, AbstractZero, AbstractThunk
export NO_FIELDS

include("compat.jl")
include("debug_mode.jl")

include("differentials/abstract_differential.jl")
include("differentials/abstract_zero.jl")
include("differentials/one.jl")
include("differentials/thunks.jl")
include("differentials/composite.jl")

include("differential_arithmetic.jl")

include("rules.jl")
include("rule_definition_tools.jl")
include("ruleset_loading.jl")

end # module


text:
# # `Palindromes<:CombinedParser`: a Tutorial for writing your combinable Parser
# Palindromes are an interesting example for parsing because
# intuitively programmers as well as laymen understand the problem:
#
# **The text is identical when read from left to right, as we are used to do,
# or when read from right to left in reverse,
# when we read only the letters and discard all non-word characters.**
#
# This example enables you to write and optimize your custom `CombinedParser` based off a minimal template.
using CombinedParsers
using CombinedParsers.Regexp
# ## 1. Regular Expression
# The PCRE test case contains nice examples of non-trivial palindromes.
#
## Defines parsers and output for pcre tests:
CombinedParsers.Regexp.@pcre_tests;

## The regular expression in the first quoted line is cryptic.
pt = pcre_test"""
/^\W*+(?:((.)\W*+(?1)\W*+\2|)|((.)\W*+(?3)\W*+\4|\W*+.\W*+))\W*+$/i
    1221
 0: 1221
 1: 1221
 2: 1
    Satan, oscillate my metallic sonatas!
 0: Satan, oscillate my metallic sonatas!
 1: <unset>
 2: <unset>
 3: Satan, oscillate my metallic sonatas
 4: S
    A man, a plan, a canal: Panama!
 0: A man, a plan, a canal: Panama!
 1: <unset>
 2: <unset>
 3: A man, a plan, a canal: Panama
 4: A
    Able was I ere I saw Elba.
 0: Able was I ere I saw Elba.
 1: <unset>
 2: <unset>
 3: Able was I ere I saw Elba
 4: A
\= Expect no match
    The quick brown fox
No match
"""

# As to why and how this PCRE pattern matches palindromes requires arcane reasoning even to the initiated:
re = Regex(pt.pattern...)
# I figure the expression is hard to construct and come up with.
# The easy part is that the pattern needs to ignore case `"i` and whitespace `\W`.
# The pattern makes intense use of backreferences and subroutines.
#
s=pt.test[3].sequence
match(re, s)
# The matched captures are purely technical (Pattern 4 is the first character).

# ### Tree display of regex
# I find it hard to understand the compact captures `(.)`, even in a nested tree display:
cp = Regcomb(pt.pattern...)
# Why no backreference `\1`, why no subroutine `(?2)`?
# Theoretical linguists, I wonder, is the minimum number of capture groups 4, for a regular expression matching palindromes?

# Writing a palindrome parser should be easier.
# And with julia compiler it should be faster.
#
# ### Regular Expression performance
# PCRE matching example 3 is fast
using BenchmarkTools
@benchmark match(re, s)
# In practice `CombinedParsers` [`Regcomb`](@ref) of the regular expression will detect palindromes too.
# Palindrome matching provides an interesting cross-parser performance benchmark.
@benchmark match(cp, s)
# `CombinedParsers.Regexp.Subroutine` matching is slow because the current implementation is using state-copies of captures.
# (TODO: Capture Subroutines could be implemented as a stack?).

# ## 2. A non-word skipping `Palindrome<:CombinedParser`
# This example of `Palindrome<:CombinedParser` is a much faster palindrome parser and more interesting and more easy to write.
# It mimics the human readable palindrome rule that is clear and quite easy to comprehend:
#
# the text is identical when read from left to right, as we are used to do,
# or when read from right to left in reverse,
# when we read only the letters and skip all non-word characters.
#
# This rule is efficient programming in natural language.
# After defining the parser, the third part of the example discusses the design of match iteration in `CombinedParsers`.
#
# ### Parsing strategy
# A custom parser needs a method to determine if there is a match and its extent at a position.
# How can this be implemented for a palindrome?
# There are two strategies:
# 1. inside-out: start at a position as `center`
#    - expand `left` and `right` from `center` until they are at word characters
#    - until word character left does not match word character at right.
#    - Succeed if a minimal length is met. Fail otherwise.
# 2. outside-in: start `left` and `right`,
#    - move positions towards `center` until they are at word characters and
#    - succeed if left and right positions meet at the center,
#    - compare these characters, and proceed to the next positions if the word characters match or fail if there is a mismatch.
#    (This might be [the-fastest-method-of-determining-if-a-string-is-a-palindrome](https://stackoverflow.com/questions/21403782/the-fastest-method-of-determining-if-a-string-is-a-palindrome).  But I figure finding all palindrome matches in a string is slow because you would be required to test for all possible substrings.)
# The inside out strategy seems easier and faster.
#
# #### Prerequisite: Skipping whitespace
# For the string `"two   words"`,  from the point of index 4 (`' '` after "from") the next word character after skipping whitespace left and right are indices of 3 (tw`o`) and 7 (`w`ords).
# In Julia syntax, this is expressed in terms of `direction` (functions `Base.prevind` and `Base.nextind` return next index to left or right), and `word_char::T`, what makes up a word character (provided method `CombinedParser.ismatch(char,parser::T)::Bool`.)

@inline function seek_word_char(direction,str,i,
                                till=lastindex(str),
                                word_char=UnicodeClass(:L))
    i=direction(str,i)
    while i>0 && i<=till && !CombinedParsers.ismatch((@inbounds str[i]),word_char)
        i=direction(str,i)
    end
    return i
end
( prev_index=seek_word_char(prevind, "two   words", 4),
  next_index=seek_word_char(nextind, "two   words", 4) )

# ### Subtyping `<: CombinedParser{STATE,RESULT}`.
# Subtyping requires you to define the type of the parsing state (for julia compiler optimizations)
# and the type of the parsing result.
STATE = NamedTuple{(:left,:center,:right),Tuple{Int,Int,Int}}
RESULT = SubString
struct Palindrome{P} <: CombinedParser{STATE,RESULT}
    word_char::P
end
Palindrome() = Palindrome(UnicodeClass(:L))

# ### Matching: `CombinedParsers._iterate`
# With the inside-out stratedy, the implementation greedily expands over non-word characters.
# Computing the first match at `posi`tion is done by this method dispatch
function CombinedParsers._iterate(x::Palindrome,
                                  str, till,
                                  posi, after,
                                  state::Nothing)
    right_ = left_ = left = right = posi
    while left>0 && right<=till &&
          lowercase(@inbounds str[left])==lowercase(@inbounds str[right])
        ## if we cannot expand, (left_,right_) succeeded
        right_ = right
        left_ = left
        left =  seek_word_char(
            prevind,str,
            left,till,x.word_char)
        right = seek_word_char(
            nextind,str,
            right,till,x.word_char)
    end
    left, left_, right_, right
    if left_ == right_
        nothing
    else
        tuple(nextind(str,right_),
              (left=left_, center=posi, right=right_))
    end
end

# Let's not iterate through next matches for now by dispatching `_iterate` for a match state to nothing.
CombinedParsers._iterate(x::Palindrome, str, till, posi, after, state) =
    nothing

# `_iterate` matches the right part of the palindrome if and only if `posi` at the center of a palindrome.
#
# The internal API calls (for the center index 18):
state = _iterate(Palindrome(),s,lastindex(s),18,18,nothing)

import CombinedParsers: _leftof, _rightof
# ### `Base.prevind` and `Base.nextind`
# `CombinedParsers` iterates through matches based on the parsing position and state.
CombinedParsers._rightof(str,i::Int,p::Palindrome,state) =
    nextind(str,state.right)
# Note that for the inside-out strategy the `Palindrome<:CombinedParser` matches from `center` and looks behind until `right`, possibly overlapping with the last match(es).
# The start index of a palindrome match is its center.
CombinedParsers._leftof(str,after::Int,p::Palindrome,state) =
    state.center

# ### `match` and `get`
# [`_iterate`](@ref) is called when the public API `match` or `parse` is used.
# Match searches for a center index and then matched the `state.center:state.right` part of the palindrome.
p = Palindrome()
m = match(p,s)

# The result of a parsing is the matching substring from `state.left:state.right`,
# implementing `Base.get` with the full argument range:
Base.get(x::Palindrome, str, till, after, posi, state) =
    SubString(str,state.left,state.right)

# The match result is matching the first palindrome, which is short and simple - but not yet what we want.
get(m)

# ### Iterating through matches
# The longest palindrome is matched too:
p = Palindrome()
[ get(m) for m in match_all(p,s) ]

# To skip trivial short palindromes we can use `Base.filter`
islong(sequence, till, posi, after, state) =
    state.right-state.left+1 > 5
long_palindrome = filter(islong,Palindrome())
get(match(long_palindrome,s))

# #### Iteration of smaller Sub-palindromes
# The set of all palindromes in a text includes the shorter palindromes contained in longer ones.
# Provide a method to iterate the previous state:
"Shrinking `state` match"
function CombinedParsers._iterate(x::Palindrome, str, till, posi, after, state)
    left_, posi_, right_ = state
    left =  seek_word_char(
        nextind,str,
        left_,till,x.word_char)
    right = seek_word_char(
        prevind,str,
        right_,till,x.word_char)
    if left >= right # left == posi
        nothing
    else
        tuple(nextind(str,right),
              (left=left, center=posi_, right=right))
    end
end

[ get(m) for m in match_all(p,s) ]

# Note that the greedy-only behaviour was atomic in terms of regular expression, which can be restored with [`Atomic`](@ref)
p = Atomic(Palindrome())
get.(match_all(p,s)) |> collect

# ### Performance Optimization
@benchmark match(long_palindrome,s)

# After writing a `CombinedParser`
## using Profile
# is recommended to investigate options to optimize.
## f(n) = for _ in 1:n; match(fast_palindrome,s); end
## @profile f(10000)
## Profile.clear()
## @profile f(1000000)
## ProfileView.view()
# ![](profile-palindrome-unicodeclass.png)\n\n
# `Combinedparsers` was optimized for minimal garbage collection (Red bars indicating garbage collection are in the `f` function).
#
# Here, replacing the `UnicodeClass` matcher (calling `categorycode`) with a fast tuple `==` check shaves off some extra time:
fast_palindrome = filter(islong,Palindrome(CharNotIn(tuple(" ,!:"...))))
@benchmark match(fast_palindrome,s)

# Note, that `CharNotIn(tuple(" ,!:"...))` is faster than the default `CharNotIn(" ,!:"...)` using `Base.Set` in case of very few `Char`s..
# Further possible optimization are
# - Caching prevind, nextind
# - Memoization
# ### Padding and combining
# Note that the PCRE pattern included outside non-words, specifically the tailing `!`.
re = Regex(pt.pattern...)
match(re,"  "*s)

# ``CombinedParsers` are designed with iteration in mind, and a small match set reduces computational time when iterating through all matches.
# `Palindrome` matches palindromes with word-char boundaries.
# The PCRE pattern includes non-words matches in the padding of palindromes, a superset of `Palindrome`.
# PCRE-equivalent matching can be achieved by combining the stricly matching `Palindrome` with parsers for the padding.
padding=Repeat(CharNotIn(fast_palindrome.parser.word_char))
match(fast_palindrome*padding*AtEnd(),s) |> get

# `Palindrome` matches from center to right, like a lookbehind parser.
# Padding to the left is not matched:
match(fast_palindrome*padding*AtEnd(),"  "*s) |> get

# Also, the left part can be parsed separately.
match(!Repeat(AnyChar())*fast_palindrome*padding*AtEnd(),"  "*s) |> get

# A prefix parser to the left requires a parser for the left-part coupled by filter:
palindrome = filter(
    Sequence(
        2,
        Lazy(Repeat(AnyChar())),
        Atomic(Palindrome()))) do sequence, till, posi, after, state
            ## posi is the start of the sequence
            posi==state[2].left
        end

# Now we can express the full pattern
p = padding * (palindrome) * padding * AtEnd()
match(p,"skipped: "*s)

# This extra complexity comes at a cost:
@benchmark match(p,$("skipped: "*s))

# ## Next...
# - optimize with memoization
# - match also palindromes with odd number of letters
# - elaborate on iteration documentation
# - comparative benchmarking, conditional on palindrome length


text:
#****************************************************************************
#                     Molecular Dynamics Potentials (MDP)
#                            CESMIX-MIT Project
#
#  Contributing authors: Ngoc-Cuong Nguyen (cuongng@mit.edu, exapde@gmail.com)
#****************************************************************************

function mdp(app)

# create exec folder if it does not exist
bindir = "exec";
cd(app.sourcepath);
if isdir(bindir) == 0
    mkdir(bindir);
end
cd(app.currentdir);

# preprocess input data
app, config = Preprocessing.preprocessing(app);

# generate code
Gencode.gencode(app); # Generate C++ code

# compile C++ source code
# cd(app.sourcepath * "C++/Main");
# Gencode.compile(app);
cd(app.sourcepath * bindir);
if isfile("CMakeCache.txt")
    rm("CMakeCache.txt", force=true);
end
if app.platform == "gpu"
    if isfile("libcpuCore.a") && isfile("libgpuCore.a") && isfile("gpuMDP")
        run(Gencode.string2cmd("cmake -D MDP_POTENTIALS=ON -D MDP_CUDA=ON ../Installation"));
    else
        run(Gencode.string2cmd("cmake -D MDP_POTENTIALS=ON -D MDP_CORES=ON -D MDP_EXECUTABLES=ON -D MDP_CUDA=ON ../Installation"));
    end
elseif app.platform == "cpu"
    if isfile("libcpuCore.a") && isfile("cpuMDP")
        run(Gencode.string2cmd("cmake -D MDP_POTENTIALS=ON ../Installation"));
    else
        run(Gencode.string2cmd("cmake -D MDP_POTENTIALS=ON -D MDP_CORES=ON -D MDP_EXECUTABLES=ON ../Installation"));
    end
end
run(Gencode.string2cmd("cmake --build ."));

# run code
#Gencode.runcode(app);
if app.platform == "cpu"
    run(Gencode.string2cmd("./cpuMDP " * app.appname * " out"));
else app.platform == "gpu"
    run(Gencode.string2cmd("./gpuMDP " * app.appname * " out"));
end
cd(app.currentdir);

if app.training > 0
    filename = app.sourcepath * "C++/Main/coefficients.bin";
    tmp = reinterpret(Float64,read(filename));
    app.coeff = reshape(tmp, 1, length(tmp));
end

return app, config

end

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
using Base.Test, Colors

@test isleaftype(eltype(colormap("Grays")))

@test_throws ArgumentError colormap("Grays", N=10)

col = distinguishable_colors(10)
@test isleaftype(eltype(col))
mindiff = Inf
for i = 1:10
    for j = i+1:10
        mindiff = min(mindiff, colordiff(col[i], col[j]))
    end
end
@test mindiff > 8


text:
# constants
const Ry2eV = 13.6058         # Ry to eV
const ANG2BOHR = 1.889725989  # angstrom to bohr


text:
using NDTools, Documenter

 # set seed fixed for documentation
DocMeta.setdocmeta!(NDTools, :DocTestSetup, :(using NDTools); recursive=true)
makedocs(modules = [NDTools],
         sitename = "NDTools.jl",
         pages = Any[
            "NDTools.jl" => "index.md",
            "Center Methods" => "center.md",
            "Select Regions" => "select_region.md",
            "Dimensionality Functions" => "dims.md",
            "Reverse views" => "reverse_view.md",
            "Base Extensions for working with ND Data" => "base_extensions.md",
         ]
        )

deploydocs(repo = "github.com/bionanoimaging/NDTools.jl.git", devbranch="main")


text:
module InterfaceTests

using Contour, Test

function setup()
    nx, ny = 10, 10

    xs = sort!(rand(nx))
    ys = sort!(rand(ny))
    zs = rand(nx, ny)

    xs, ys, zs
end

xs, ys, zs = setup()

cs = @inferred contours(xs, ys, zs)
for c in levels(cs)
    for l in lines(c)
        x, y = coordinates(l)
        @assert typeof(x) == typeof(y) == Vector{Float64}
    end
end

end # InterfaceTests


text:
using Test
using FileIO
using LibSndFile
using SampledSignals
using LinearAlgebra

extensions = (".wav",
              ".ogg",
              ".flac")

formats = [
           ("_wav", format"WAV"),
           ("_ogg", format"OGG"),
           ("_flac", format"FLAC"),
          ]

include("utils.jl")

@testset "LibSndFile Tests" begin
  @testset "Version" begin
    v = LibSndFile.version()
    @test v[1:10] == "libsndfile"
  end
  @testset "File reading" begin
    include("file_reading.jl")
  end
  @testset "Streaming and IO" begin
    include("stream.jl")
  end
  @testset "File writing" begin
    include("file_writing.jl")
  end
  @testset "FileIO Integration" begin
    include("fileio.jl")
  end
  @testset "Display" begin
    include("display.jl")
  end
  # TODO: check out what happens when samplerate, channels, etc. are wrong
  # when reading/writing
end


text:
#=
fbp-par.jl
Simple interfaces to parallel-beam FBP for user convenience
=#

using MIRT: sino_geom, image_geom # todo
#using Sinograms: fbp2, RealU

export fbp, fbp!

"""
    fbp(sino::AbstractMatrix ; kwargs...)

FBP reconstruction from a parallel-beam sinogram of size `[nr × nϕ]`.
Returns an image of size `[nx × ny]`.

# Options
* `nx` : default `nr`
* `ny` : default `nx`
* `kwargs` : passed to `fbp2`
"""
function fbp(
    sino::AbstractMatrix{T} ;
    nx::Int = size(sino, 1),
    ny::Int = nx,
    kwargs...
) where {T <: Number}
    image = zeros(promote_type(T, Float32), nx, ny)
    fbp!(image, sino; kwargs...)
end

"""
    fbp!(image, sino ; orbit::Real = 180, kwargs...)

FBP reconstruction from a parallel-beam sinogram of size `[nr × nϕ]`.
Writes result into `image` matrix.

# Input
* `sino::AbstractMatrix`

# Options for `sino_geom`
* `dr` : sinogram radial spacing; default 1
* `orbit` : angular range in degrees; default 180
* `orbit_start` : angular range in degrees; default 180

# Options for `image_geom`
* `dx`, `dy`, `offset_x`, `offset_y`

# Options
* `kwargs` : passed to `fbp2`

# Output
* `image::AbstractMatrix` is mutated
"""
function fbp!(
    image::AbstractMatrix{<:Number},
    sino::AbstractMatrix{<:Number} ;
    dr::RealU = 1,
    orbit::Real = 180,
    orbit_start::Real = 0,
    dx::RealU = dr,
    dy::RealU = dx,
    offset_x::Real = 0,
    offset_y::Real = 0,
    kwargs...
)
    nb, na = size(sino)
    nx, ny = size(image)
    sg = sino_geom(:par ; nb, na, d = dr, orbit, orbit_start)
    ig = image_geom(; nx, ny, dx, dy, offset_x, offset_y)
#   plan = FBPplan(sg, ig)
    plan = fbp2(sg, ig) # todo FBPplan ??
#   fbp2!(image, sino, plan) # todo
    tmp, _ = fbp2(plan, sino)
    copyto!(image, tmp)
    return image
end


text:
"""
    JDF.load(indir, verbose = true)

    JDF.load(indir, cols = Vector{Symbol}, verbose = true)

Load a `Tables.jl` table from JDF saved at `outdir`. On Julia > v1.3, a multithreaded
version is used.
"""
load(indir; cols = Symbol[], verbose = false) = begin
    # starting from DataFrames.jl 0.21 the colnames are strings
    cols = collect(string.(cols))
    metadatas = jdfmetadata(indir)

    # TODO simplify this this is duplicated in load_columns
    if length(cols) == 0
        cols =collect(string.(metadatas.names))
    else
        scmn = setdiff(cols,collect(string.(metadatas.names)))
        if length(scmn) > 0
            throw("columns $(reduce((x,y) -> string(x) * ", " * string(y), scmn)) are not available, please ensure you have spelt them correctly")
        end
    end

    cols_in_loaded_order, result_vectors = load_columns(indir; cols = cols, verbose = verbose)

    # reorders to specified order
    reorder_idx = indexin(cols_in_loaded_order, cols)

    Table(NamedTuple{Tuple(Symbol.(cols))}(@view result_vectors[reorder_idx]))
end

load(jdf::JDFFile; args...) = load(path(jdf); args...)

loadjdf(args...; kwargs...) = load(args...; kwargs...)

max_stars_repo_name: xiaodaigh/JDF

max_issues_repo_name: xiaodaigh/JDF

max_forks_repo_name: xiaodaigh/JDF


text:
using BangBang

struct UndefArray{N, T <: NTuple{N, Int}} <: AbstractArray{Union{}, N}
    size::T
end

Base.size(arr::UndefArray) = arr.size
Base.isassigned(::UndefArray, ::Integer...) = false

function BangBang.NoBang._setindex(x::UndefArray{N}, v, I...) where N
    ys = Array{typeof(v), N}(undef, x.size...)
    ys[I...] = v
    return ys
end

mapbb(f, xs::Array, inbounds::Val = Val(true)) =
    __map!!(f, UndefArray(size(xs)), xs, 1, inbounds)

@inline function __map!!(f, ys, xs, n, inbounds)
    @simd for i in n:length(ys)
        if inbounds === Val(true)
            @inbounds begin
                ys′ = setindex!!(ys, f(xs[i]), i)
            end
        else
            ys′ = setindex!!(ys, f(xs[i]), i)
        end
        if !(ys′ isa typeof(ys))
            return __map!!(f, ys′, xs, i + 1, inbounds)
        end
    end
    return ys
end


text:
function AngularTroposphericDelayModifier(arg0::DiscreteTroposphericModel)
    return AngularTroposphericDelayModifier((DiscreteTroposphericModel,), arg0)
end

function get_parameters_drivers(obj::AngularTroposphericDelayModifier)
    return jcall(obj, "getParametersDrivers", List, ())
end

function modify(obj::AngularTroposphericDelayModifier, arg0::EstimatedMeasurement)
    return jcall(obj, "modify", void, (EstimatedMeasurement,), arg0)
end


text:
"""
Module to produce incidence plots

* Developer: Hilding Elmqvist, Mogram AB
* First version: December 2020
* License: MIT (expat)

"""
module IncidencePlot

using GLMakie
using AbstractPlotting
using GeometryBasics
using Colors

# include("...ModiaBase/src/Tearing.jl")
# using .Tearing

export showIncidence

removeBlock(ex) = ex
function removeBlock(ex::Expr)
    if ex.head in [:block]
        ex.args[1]
    else
        Expr(ex.head, [removeBlock(arg) for arg in ex.args]...)
    end
end

AbstractPlotting.inline!(false)

"""
    showIncidence(equations, unknowns, G, vActive, substitutions, assigned, assignedVar, blocks; sortVariables=false)

Shows the incidence matrix using the plotting package Makie. If `assigned` is not [], the assigned variable is marked. `blocks` is a vector of vectors
enabling showing of BLT-blocks. After each incidence matrix is output, the user is promted for a RETURN, not to immediately overwriting the plot.

"""
function showIncidence(equations, unknowns, G, vActive, substitutions, assigned, assignedVar, blocks; sortVariables=false)
    scene = Scene(camera = campixel!, show_axis = false, resolution = (1000, 800))
    margin = 20
    size = 600

    nVar = length(unknowns)
    if sortVariables
        sortedVariables = fill(0, nVar)
        varCount = 0
        for b in blocks
            for e in b
                varCount += 1
                sortedVariables[varCount] = assignedVar[e]
            end
        end
        for i in 1:length(vActive)
            if ! vActive[i]
                varCount += 1
                sortedVariables[varCount] = i
            end
        end
    else
        sortedVariables = 1:nVar
    end

    n=length(equations)
    m=length(unknowns)

    width = min(div(size, n), 20)
    height = min(div(size, m), 20)

    function drawText(s, row, column; align=:left, vertical=false)
        rot = if ! vertical; 0.0 else pi/2 end
        text!(scene, s, position = Vec(margin + column*width + width/2, margin + (n+1-row)*height + height/2),
            color = :green, textsize = 15, align = (align, :center), rotation = rot, show_axis = false)
    end

    # -----------------------------------------------------------------------------

    # Draw matrix grid
    for i in 1:n+1
        poly!([Point2f0(margin + width, margin + i*height), Point2f0(margin + (m+1)*width, margin + i*height)], strokecolor = :black)
    end

    for j in 1:m+1
        poly!([Point2f0(margin + j*width, margin + height), Point2f0(margin + j*width, margin + (n+1)*height)], strokecolor = :black)
    end

    for j in 1:length(unknowns)
        # Draw variable number under matrix
        drawText(string(sortedVariables[j]), n+1, j, align=:right, vertical=true)

        # Draw variable names above matrix
        v = unknowns[sortedVariables[j]]
        var = string(v)
        if v in keys(substitutions)
            subs = replace(replace(replace(string(substitutions[v]), "(" => ""), ")" => ""), "--" => "")
            var *= " = " * subs
        end
        drawText(var, 0, j, vertical=true)
    end

    row = 0
    for b in blocks
        for e in b
            row += 1
            # Draw equation number to the left
            drawText(string(e), row, 0, align=:right)

            # Draw equation to the right
            equ = equations[e]
            Base.remove_linenums!(equ)
            equ = removeBlock(equ)
            equ = string(equ)
            equ = replace(equ, "\n" => ";")
            if equ != :(0 = 0)
                drawText(equ, row, m+1)
            end
        end
    end

    row = 0
    for b in blocks
        if length(b) > 1
#=
            println("\nTearing of BLT block")
            @show b
            es::Array{Int64,1} = b
            vs = assignedVar[b]
            td = EquationTearing(G, length(unknowns))
            (eSolved, vSolved, eResidue, vTear) = tearEquations!(td, (e,v) -> v in G[e], es, vs)
            @show eSolved vSolved eResidue vTear
            @show equations[eSolved] unknowns[vSolved] equations[eResidue] unknowns[vTear]
            b = vcat(eSolved, eResidue)
            @show b
            blocksVars = vcat(vSolved, vTear)
=#
            color = :lightgrey
            poly!(Rect(Vec(margin + (row+1)*width, margin + (n+1-row)*height-length(b)*height),
                Vec(length(b)*width, length(b)*height)), color = color, strokecolor = :black, transperancy=true)
        end
        for e in b
            row += 1
            incidence = fill(false, nVar)
            for v in G[e]
                if sortedVariables[v] > 0
                    incidence[sortedVariables[v]] = true
                end
            end

            col = 0
            for v in sortedVariables
                col += 1
                ins = if v == 0 || ! vActive[v]; 0 elseif sortedVariables[v] > 0 && ! incidence[sortedVariables[v]]; 0 elseif v > 0 && v<= length(assigned) && assigned[v] == e; 1 else 2 end
                if ins > 0
                    color = if ins == 1; :red else :blue end
                    poly!(Rect(Vec(margin + col*width, margin + (n+1-row)*height), Vec(width, height)), color = color, strokecolor = :black)
                end
            end
        end
    end

    display(scene)

    println("Continue?")
    readline()
end

function testShowIncidence()
    # RCRL circuit
    equations = Expr[:(ground.p.v = V.n.v), :(V.p.v = Ri.n.v), :(Ri.p.v = R.p.v), :(R.n.v = C.p.v), :(C.n.v = ground.p.v), :(ground.p.i + V.n.i + C.n.i = 0), :(V.p.i + Ri.n.i = 0), :(Ri.p.i + R.p.i = 0), :(R.n.i + C.p.i = 0), :(Ri.v = Ri.p.v - Ri.n.v), :(0 = Ri.p.i + Ri.n.i), :(Ri.i = Ri.p.i), :(Ri.R * Ri.i = Ri.v), :(R.v = R.p.v - R.n.v), :(0 = R.p.i + R.n.i), :(R.i = R.p.i), :(R.R * R.i = R.v), :(C.v = C.p.v - C.n.v), :(0 = C.p.i + C.n.i), :(C.i = C.p.i), :(C.C * der(C.v) = C.i), :(V.v = V.p.v - V.n.v), :(0 = V.p.i + V.n.i), :(V.i = V.p.i), :(V.v = 10), :(ground.p.v = 0)]
    unknowns = Any[:(Ri.v), :(Ri.i), :(Ri.p.i), :(Ri.p.v), :(Ri.n.i), :(Ri.n.v), :(R.v), :(R.i), :(R.p.i), :(R.p.v), :(R.n.i), :(R.n.v), :(C.v), :(der(C.v)), :(C.i), :(C.p.i), :(C.p.v), :(C.n.i), :(C.n.v), :(V.v), :(V.i), :(V.p.i), :(V.p.v), :(V.n.i), :(V.n.v), :(ground.p.i), :(ground.p.v)]
    G = [[27, 25], [23, 6], [4, 10], [12, 17], [19, 27], [26, 24, 18], [22, 5], [3, 9], [11, 16], [1, 4, 6], [3, 5], [2, 3], [2, 1], [7, 10, 12], [9, 11], [8, 9], [8, 7], [13, 17, 19], [16, 18], [15, 16], [14, 15], [20, 23, 25], [22, 24], [21, 22], [20], [27]]
    vActive = Bool[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
    substitutions = Dict{Any,Any}()
    assigned = Any[]
    assignedVar = Any[]
    blocks = [[1], [2], [3], [4], [5], [6], [7], [8], [9], [10], [11], [12], [13], [14], [15], [16], [17], [18], [19], [20], [21], [22], [23], [24], [25], [26]]
    sortVariables = false

    showIncidence(equations, unknowns, G, vActive, substitutions, assigned, assignedVar, blocks, sortVariables=sortVariables)

    # RCRL
    equations = Expr[:(ground.p.v = V.n.v), :(V.p.v = Ri.n.v), :(Ri.p.v = R.p.v), :(R.n.v = C.p.v), :(C.n.v = ground.p.v), :(ground.p.i + V.n.i + C.n.i = 0), :(V.p.i + Ri.n.i = 0), :(Ri.p.i + R.p.i = 0), :(R.n.i + C.p.i = 0), :(Ri.v = Ri.p.v - Ri.n.v), :(0 = Ri.p.i + Ri.n.i), :(Ri.i = Ri.p.i), :(Ri.R * Ri.i = Ri.v), :(R.v = R.p.v - R.n.v), :(0 = R.p.i + R.n.i), :(R.i = R.p.i), :(R.R * R.i = R.v), :(C.v = C.p.v - C.n.v), :(0 = C.p.i + C.n.i), :(C.i = C.p.i), :(C.C * der(C.v) = C.i), :(V.v = V.p.v - V.n.v), :(0 = V.p.i + V.n.i), :(V.i = V.p.i), :(V.v = 10), :(ground.p.v = 0)]
    unknowns = Any[:(Ri.v), :(Ri.i), :(Ri.p.i), :(Ri.p.v), :(Ri.n.i), :(Ri.n.v), :(R.v), :(R.i), :(R.p.i), :(R.p.v), :(R.n.i), :(R.n.v), :(C.v), :(der(C.v)), :(C.i), :(C.p.i), :(C.p.v), :(C.n.i), :(C.n.v), :(V.v), :(V.i), :(V.p.i), :(V.p.v), :(V.n.i), :(V.n.v), :(ground.p.i), :(ground.p.v)]
    G = [[27, 25], [23, 6], [4, 10], [12, 17], [19, 27], [26, 24, 18], [22, 5], [3, 9], [11, 16], [1, 4, 6], [3, 5], [2, 3], [2, 1], [7, 10, 12], [9, 11], [8, 9], [8, 7], [13, 17, 19], [16, 18], [15, 16], [14, 15], [20, 23, 25], [22, 24], [21, 22], [20], [27]]
    vActive = Bool[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
    substitutions = Dict{Any,Any}()
    assigned = [10, 13, 12, 3, 11, 2, 17, 16, 8, 14, 15, 4, 0, 21, 20, 9, 18, 19, 5, 25, 24, 7, 22, 23, 1, 6, 26]
    assignedVar = [25, 6, 4, 12, 19, 26, 22, 9, 16, 1, 5, 3, 2, 10, 11, 8, 7, 17, 18, 15, 14, 23, 24, 21, 20, 27, 0]
    blocks = [[26], [1], [25], [22], [2], [5], [18], [4], [10, 13, 12, 8, 16, 17, 14, 3], [11], [7], [23], [15], [9], [19], [6], [20], [21], [24]]
    sortVariables = true
    showIncidence(equations, unknowns, G, vActive, substitutions, assigned, assignedVar, blocks, sortVariables=sortVariables)

end

# testShowIncidence()

end


text:
using .DynamicalSystems
import .DynamicalSystems: lyapunov, lyapunovspectrum

"""
Whip up a DynamicalSystem from a Process
"""
function process2ds(P::Process)
    prob = process2problem(P)
    if getalg(P) isa FunctionMap
        d = DiscreteDynamicalSystem(prob.f, prob.u0, prob.p)
    else
        d = ContinuousDynamicalSystem(prob)
    end
end
export process2ds

"""
Calculate the largest lyapunov exponent of a process
"""
function lyapunov(P::Process)
    d = process2ds(P)
    T = gettmax(P)
    Ttr = gett0(P) - gettransient_t0(P)
    lyapunov(d, T; Ttr)
end
export lyapunov

""" Calculate the lyapunov spectrum"""
function lyapunovspectrum(P::Process, N::Integer=1000, k::Integer=length(getX0(P)))
    d = process2ds(P)
    u0 = getX0(P)
    Ttr = gett0(P) - gettransient_t0(P)
    #dt = getdt(P)
    lyapunovspectrum(d, N, k; u0, Ttr)
end
export lyapunovspectrum

"""Calculate the largest lyapunov exponent as a function of a parameter"""
function lyapunovresponse(P::Process, p, prange, N=1000, k=length(getX0(P)))
    d = process2ds(P)
    λs = Array{Float64, 2}(undef, length(prange), k)
    Threads.@threads for ip ∈ 1:length(prange)
        subP = updateparam(P, p, constantParameter, prange[ip])
        if k == 1 && !(d isa DiscreteDynamicalSystem)
            λ = lyapunov(subP)
        else
            λ = lyapunovspectrum(subP, N, k)
        end
        λs[ip, :] .= λ
    end
    return λs
end
export lyapunovresponse


text:
using Ansillary
using Ansillary.Inputs
using Documenter
using JuliaFormatter
using Test

# Style Tests!

function isformatted(path)
    if isfile(path)
        if last(splitext(path)) != ".jl" || success(`git check-ignore $path`)
            return true
        end

        original = read(path, String)
        formatted = format_text(
            original,
            always_for_in = true,
            whitespace_typedefs = true,
            whitespace_ops_in_indices = true,
            remove_extra_newlines = true,
        )

        original == formatted
    else
        isformatted(joinpath.(path, readdir(path)))
    end
end

isformatted(paths::Vector) = all(isformatted.(paths))

@testset "Formatting" begin
    let tests = dirname(@__FILE__), project = dirname(tests)
        @test isformatted(project)
    end
end

# Doc Tests!

DocMeta.setdocmeta!(Ansillary, :DocTestSetup, :(using Ansillary))
DocMeta.setdocmeta!(
    Ansillary.Cursor,
    :DocTestSetup,
    :(using Ansillary; using Ansillary.Cursor),
)
DocMeta.setdocmeta!(
    Ansillary.Inputs,
    :DocTestSetup,
    :(using Ansillary; using Ansillary.Inputs),
)
DocMeta.setdocmeta!(
    Ansillary.Screen,
    :DocTestSetup,
    :(using Ansillary; using Ansillary.Screen),
)
doctest(Ansillary)

# Unit Tests!

# TODO: Use libvterm (http://www.leonerd.org.uk/code/libvterm/) to implement tests.

@testset "Modifiers" begin
    @test Alt() + Left() == Modified(Left(), [Alt()])
    @test Ctrl() + 'c' == CTRL_C == Modified(Character('c'), [Ctrl()])
    @test Ctrl() + Alt() + Delete() == Modified(Delete(), [Ctrl(), Alt()])
    @test Shift() + Super() + 'y' == Modified(Character('y'), [Shift(), Super()])
    @test Ctrl() + Alt() + Super() + F(1) == Modified(F(1), [Alt(), Super(), Ctrl()])
    @test Ctrl() + Alt() + Shift() + Super() + Character('q') ==
          Modified(Character('q'), [Shift(), Alt(), Super(), Ctrl()])
    @test Super() + Ctrl() + Alt() + Shift() + Super() + Character('q') ==
          Modified(Character('q'), [Shift(), Alt(), Super(), Ctrl()])
end


text:
# TODO: We should probably define these in a straightforward way, rather than preserve
# generality in the definitions.
"""
    z0

`z0` is `[1, 0]`.
"""
const z0 = [1, 0]

"""
    z1

`z1` is `[0, 1]`.
"""
const z1 = [0, 1]

####
#### Single qubit gates
####

"""
    I2

The single-qubit identity operator.
"""
const I2 = [1 0; 0 1]

"""
    X

The Pauli `X` gate.
"""
const X = [0 1; 1 0]

"""
    Y

The Pauli `Y` gate.
"""
const Y = [0 -im; im 0]

"""
    Z

The Pauli `Z` gate.
"""
const Z = [1 0; 0 -1]

"""
    H

The single-qubit Hadamard gate.
"""
const H = (1/sqrt(2)) * [1 1; 1 -1]

"""
    S

The phase gate.
"""
const S = [1 0; 0 im]

"""
    T

The π/8 gate.
"""
const T = [1 0; 0 (cospi(1/4) + im * sinpi(1/4))]

# FIXME: These are easy to construct. Maybe we do not need
# to provide them. Decide.
"""
    sqrt_NOT

The square root of the `NOT` (or `X`) gate.
"""
const sqrt_NOT = sqrt(complex(X))

# FIXME: Highlighting of $cgate is broken
# Controlled X, Y, Z, H with one and two control qubits
for gate in (:X, :Y, :Z, :H)
    for n in 1:2
        cgate = Symbol(fill(:C, n)..., gate)
        gatestr = string(gate)
        cgatestr = string(cgate)
        expl = n == 2 ? "gate with two control qubits." : "gate."
        @eval begin
"""
   $($cgatestr)

The controlled-`$($gatestr)` $($expl)
"""
        const $cgate = control($gate, $n)
        end
    end
end

"""
    SWAP

The 2-qubit `SWAP` gate.
"""
const SWAP = [1 0 0 0
              0 0 1 0
              0 1 0 0
              0 0 0 1]

"""
    CSWAP

The controlled `SWAP` gate.
"""
const CSWAP = control(SWAP)

"""
    iSWAP

The 2-qubit `iSWAP` gate.
"""
const iSWAP = [1 0  0  0
               0 0  im 0
               0 im 0  0
               0 0  0  1]

"""
    Magic basis

The matrix whose columns are the magic basis.
"""
const magic_basis = (1 / sqrt(2)) *
    [1  0    0   im
     0  im   1    0
     0  im  -1    0
     1  0    0  -im]


text:
function get_parameters_from_file(fname::String)

    #Open file
    fid         = h5open( fname , "r")

    #Dissipation on?
    diss = read_attribute(fid, "Dissipation")

    if( diss == "true")

        p0 = Param(

            #Temporal variables
            t_sim_init      = read_attribute(fid, "t_min")          ,
            t_sim_final     = read_attribute(fid, "t_max")          ,
            out_every_t     = read_attribute(fid, "out_every_t")      ,
            deltat          = read_attribute(fid, "dt")             ,

            max_runtime     = "24:00:00"  ,

            #Spacial variables
            xmin           =    read_attribute(fid, "xmin")     ,
            xmax           =    read_attribute(fid, "xmax")     ,
            xnodes         =    read_attribute(fid, "xnodes")   ,

            ymin           =    read_attribute(fid, "ymin")     ,
            ymax           =    read_attribute(fid, "ymax")     ,
            ynodes         =    read_attribute(fid, "ynodes")   ,

            #Initial Conditions
            A0              =  read_attribute(fid , "A0")       ,
            σ               =  read_attribute(fid , "sigma")    ,
            r0              =  read_attribute(fid , "r0")       ,
            ω               =  read_attribute(fid , "omega")    ,
            m               =  read_attribute(fid , "m_mode")   ,

            #Klein Gordon Mass
            μ               =    read_attribute(fid , "field_mass")         ,

            #Dissipation terms
            dissipation     =  :true      ,

            R_orbit         =  read_attribute(fid , "Rorbit")       ,
            R1              =  read_attribute(fid , "R1")         ,
            R2              =  read_attribute(fid , "R2")         ,

            Ω               =  read_attribute(fid , "Omega")        ,
            alpha1          =  read_attribute(fid , "alpha1")         ,
            alpha2          =  read_attribute(fid , "alpha2")         , #Periodic BC?

            #Boundaries
            Boundaries      = Symbol(read_attribute(fid, "Boundaries"))       , # :square :periodic

            )

        close(fid)

        return p0

    else

        p0 = Param(

             #Temporal variables
             t_sim_init      = read_attribute(fid, "t_min")          ,
             t_sim_final     = read_attribute(fid, "t_max")          ,
             out_every_t     = read_attribute(fid, "out_every_t")      ,
             deltat          = read_attribute(fid, "dt")             ,

            max_runtime     = "24:00:00"  ,

            #Spacial variables
            xmin           =    read_attribute(fid, "xmin")     ,
            xmax           =    read_attribute(fid, "xmax")     ,
            xnodes         =    read_attribute(fid, "xnodes")   ,

            #Spacial variables
            ymin           =    read_attribute(fid, "ymin")     ,
            ymax           =    read_attribute(fid, "ymax")     ,
            ynodes         =    read_attribute(fid, "ynodes")   ,

            #Initial Conditions
            A0              =  read_attribute(fid , "A0")       ,
            σ               =  read_attribute(fid , "sigma")    ,
            r0              =  read_attribute(fid , "r0")       ,
            ω               =  read_attribute(fid , "omega")    ,
            m               =  read_attribute(fid , "m_mode")   ,

            #Klein Gordon Mass
            μ               =    read_attribute(fid , "field_mass")         ,

            #Dissipation terms
            dissipation     =  :false        ,

            #Periodic BC?
            Boundaries      = Symbol(read_attribute(fid, "Boundaries"))       , # :square :periodic

            )

            close(fid)

            return p0

    end

end


text:
using Random
using LinearAlgebra
using Statistics
using Optim
using DataFrames
using DataFramesMeta
using CSV
using HTTP
using GLM

# read in function to create state transitions for dynamic model
#include("create_grids.jl")
include("create_grids.jl")

#:::::::::::::::::::::::::::::::::::::::::::::::::::
# Question 1: reshaping the data
#:::::::::::::::::::::::::::::::::::::::::::::::::::
# load in the data
url = "https://raw.githubusercontent.com/OU-PhD-Econometrics/fall-2021/master/ProblemSets/PS5-ddc/busdataBeta0.csv"
df = CSV.read(HTTP.get(url).body, DataFrame)

# create bus id variable
df = @transform(df, bus_id = 1:size(df,1))

#---------------------------------------------------
# reshape from wide to long (must do this twice be-
# cause DataFrames.stack() requires doing it one
# variable at a time)
#---------------------------------------------------
# first reshape the decision variable
dfy = @select(df, :bus_id,:Y1,:Y2,:Y3,:Y4,:Y5,:Y6,:Y7,:Y8,:Y9,:Y10,:Y11,:Y12,:Y13,:Y14,:Y15,:Y16,:Y17,:Y18,:Y19,:Y20,:RouteUsage,:Branded)
dfy_long = DataFrames.stack(dfy, Not([:bus_id,:RouteUsage,:Branded]))
rename!(dfy_long, :value => :Y)
dfy_long = @transform(dfy_long, time = kron(collect([1:20]...),ones(size(df,1))))
select!(dfy_long, Not(:variable))

# next reshape the odometer variable
dfx = @select(df, :bus_id,:Odo1,:Odo2,:Odo3,:Odo4,:Odo5,:Odo6,:Odo7,:Odo8,:Odo9,:Odo10,:Odo11,:Odo12,:Odo13,:Odo14,:Odo15,:Odo16,:Odo17,:Odo18,:Odo19,:Odo20)
dfx_long = DataFrames.stack(dfx, Not([:bus_id]))
rename!(dfx_long, :value => :Odometer)
dfx_long = @transform(dfx_long, time = kron(collect([1:20]...),ones(size(df,1))))
select!(dfx_long, Not(:variable))

# join reshaped df's back together
df_long = leftjoin(dfy_long, dfx_long, on = [:bus_id,:time])
sort!(df_long,[:bus_id,:time])

#:::::::::::::::::::::::::::::::::::::::::::::::::::
# Question 2: estimate a static version of the model
#:::::::::::::::::::::::::::::::::::::::::::::::::::

res = glm(@formula(Y~Branded+Odometer),df_long, Binomial(), LogitLink())
println(res)

#:::::::::::::::::::::::::::::::::::::::::::::::::::
# Question 3a: read in data for dynamic model
#:::::::::::::::::::::::::::::::::::::::::::::::::::

# [ Load in the data ]
url = "https://raw.githubusercontent.com/OU-PhD-Econometrics/fall-2021/master/ProblemSets/PS5-ddc/busdata.csv"
df = CSV.read(HTTP.get(url).body, DataFrame)
#show(df, allcols=true)

Y = Matrix(df[:,[:Y1,:Y2,:Y3,:Y4,:Y5,:Y6,:Y7,:Y8,:Y9,:Y10,:Y11,:Y12,:Y13,:Y14,:Y15,:Y16,:Y17,:Y18,:Y19,:Y20]])

# [ convert other data frame columns to matrices]
O = Matrix(df[:,[:Odo1,:Odo2,:Odo3,:Odo4,:Odo5,:Odo6,:Odo7,:Odo8,:Odo9,:Odo10,:Odo11,:Odo12,:Odo13,:Odo14,:Odo15,:Odo16,:Odo17,:Odo18,:Odo19,:Odo20]])
X = Matrix(df[:,[:Xst1,:Xst2,:Xst3,:Xst4,:Xst5,:Xst6,:Xst7,:Xst8,:Xst9,:Xst10,:Xst11,:Xst12,:Xst13,:Xst14,:Xst15,:Xst16,:Xst17,:Xst18,:Xst19,:Xst20]])
Z = Matrix(df[:,[:Zst]])
B = Matrix(df[:,[:Branded]])

#:::::::::::::::::::::::::::::::::::::::::::::::::::
# Question 3b: generate state transition matrices
#:::::::::::::::::::::::::::::::::::::::::::::::::::

#discretize the mileage x1 transitions into 0,125-mile bins (i.e. 0.125 units of x1t ).
#specify x2 as a discrete uniform distribution ranging from 0.25 to 1.25 with 0.01 unit increments.
zval,zbin,xval,xbin,xtran = create_grids()

T = 20
FV = zeros((size(xtran,1),2,T+1))
beta = 0.9
theta = [0.1,0.2,0.3]

theta = ones(3)
for t = 1:T
    for b =0:1
        for z = 1:zbin
            for x = 1:xbin
                t = (T+1)-t
                row = x + (z-1)*xbin
                #not replace the engine
                v1t = theta[1] +theta[2]*xval[x]+theta[3]*b+ (xtran[row,:]'*FV[(z-1)*xbin+1:z*xbin,b+1,t+1])
                #replace the engine
                v0t = xtran[1+(z-1)*xbin,:]'*FV[(z-1)*xbin+1:z*xbin,b+1,t+1]
                FV[row,b+1,t] = beta*log(exp(v0t)+exp(v1t))
            end
        end
    end
end

loglik = 0
for t = 1:T
    for i = 1:size(X,1)
        t = (T+1)-t
        z = Z[i]
        x = X[i,t]
        b = B[i]
        row1 = x + (z-1)*xbin
        row0 = 1+(z-1)*xbin
        v_diff  = theta[1] +theta[2]*O[i,t]+theta[3]*b + (xtran[row1,:].-xtran[row0,:])'*FV[row0:row0+xbin-1,b+1,t+1]
        if Y[i,t] == 1
            loglik += log(exp(v_diff)/(1+exp(v_diff)))*Y[i,t]
        elseif Y[i,t] == 0
            loglik += log(1/(1+exp(v_diff)))*Y[i,t]
        end
    end
end
loglik = -loglik

@views @inbounds function myfun()

    Y = Matrix(df[:,[:Y1,:Y2,:Y3,:Y4,:Y5,:Y6,:Y7,:Y8,:Y9,:Y10,:Y11,:Y12,:Y13,:Y14,:Y15,:Y16,:Y17,:Y18,:Y19,:Y20]])

    # [ convert other data frame columns to matrices]
    O = Matrix(df[:,[:Odo1,:Odo2,:Odo3,:Odo4,:Odo5,:Odo6,:Odo7,:Odo8,:Odo9,:Odo10,:Odo11,:Odo12,:Odo13,:Odo14,:Odo15,:Odo16,:Odo17,:Odo18,:Odo19,:Odo20]])
    X = Matrix(df[:,[:Xst1,:Xst2,:Xst3,:Xst4,:Xst5,:Xst6,:Xst7,:Xst8,:Xst9,:Xst10,:Xst11,:Xst12,:Xst13,:Xst14,:Xst15,:Xst16,:Xst17,:Xst18,:Xst19,:Xst20]])
    Z = Matrix(df[:,[:Zst]])
    B = Matrix(df[:,[:Branded]])

    #discretize the mileage x1 transitions into 0,125-mile bins (i.e. 0.125 units of x1t ).
    #specify x2 as a discrete uniform distribution ranging from 0.25 to 1.25 with 0.01 unit increments.
    zval,zbin,xval,xbin,xtran = create_grids()

    T = 20
    FV = zeros((size(xtran,1),2,T+1))
    beta = 0.9

    function ddc(theta,Y,O,X,Z,B)
        #(c)
        for t = 1:T
            for b =0:1
                for z = 1:zbin
                    for x = 1:xbin
                        t = (T+1)-t
                        row = x + (z-1)*xbin
                        #not replace the engine
                        v1t = theta[1] +theta[2]*xval[x]+theta[3]*b+ (xtran[row,:]'*FV[(z-1)*xbin+1:z*xbin,b+1,t+1])
                        #replace the engine
                        v0t = xtran[1+(z-1)*xbin,:]'*FV[(z-1)*xbin+1:z*xbin,b+1,t+1]
                        FV[row,b+1,t] = beta*log(exp(v0t)+exp(v1t))
                    end
                end
            end
        end
        #(d)
        loglik = 0
        for t = 1:T
            for i = 1:size(X,1)
                t = (T+1)-t
                z = Z[i]
                x = X[i,t]
                b = B[i]
                row1 = x + (z-1)*xbin
                row0 = 1+(z-1)*xbin
                v_diff  = theta[1] +theta[2]*O[i,t]+theta[3]*b + (xtran[row1,:].-xtran[row0,:])'*FV[row0:row0+xbin-1,b+1,t+1]
                if Y[i,t] == 1
                    loglik += log(exp(v_diff)/(1+exp(v_diff)))*Y[i,t]
                elseif Y[i,t] == 0
                    loglik += log(1/(1+exp(v_diff)))*Y[i,t]
                end
            end
        end
        loglik = -loglik
        return loglik
    end

    theta0 = rand(3)
    obj = optimize(theta->ddc(theta, Y,O,X,Z,B), theta0, LBFGS(), Optim.Options(g_tol=1e-5, iterations=10, show_trace=true, show_every=50))
    theta = obj.minimizer
    println("theta is ",theta)

end

myfun()


text:
#!/usr/bin/env julia

# A sript to draw for secret santa and send pairings via e-mail using msmtp.

# Usage: julia secretsanta.jl <table> [<mail account>]

# <table> is a csv table with columns "Name", "Pair", and "Address".
# Use the "Pair" columns to prevent drawing, e.g, partners names by assigning
# partners an identical id.

using Random
using CSV
using Dates

"""
    readtable(f)

Reads a csv file f containing information of secret santa participants.
Order of columns is: names, pairs, e-mail addresses.
"""
function readtable(tablefile)
    if isfile(tablefile)
        table = CSV.read(tablefile, comment="#")
    else
        error("$tablefile is not a valid file.")
    end
    table
end

"""
    getsecretsanta(n, p)

Takes two vectors with names, n, and corresponding pairings, p,
and returns a permutated names vector for secret santa.
"""
function getsecretsanta(names, pairs)
    santaindex = shuffle(eachindex(names))
    while any(pairs.== pairs[santaindex])
        santaindex = shuffle(eachindex(names))
    end
    secretsanta = names[santaindex]
end

"""
    mailsecretsanta(t, a)

Notifies participants via e-mail using information in table `t`.
This requires a set-up msmtp. Choose your prefered mail account `a`.
"""
function mailsecretsanta(table, mailaccount = "default")
    if length(size(table)) > 1
        for i in 1:size(table, 1)
            name = table.Name[i]
            mail = table.Address[i]
            santa = table.Santa[i]
            message = "To: $mail\r\nSubject: Weihnachtswichteln $(Dates.year(Dates.today()))\r\n\r\nHallo $name,\r\n\r\nDu darfst heuer $santa bewichteln.\r\n\r\nLG\r\nDein Weihnachtswichtelskript\r\n\r\nBesuch mich auf https://github.com/lleiding/secretsanta.jl\r\n"
            run(pipeline(`echo $message`, `msmtp -a $mailaccount $mail`))
        end
    else
        name = table.Name
        mail = table.Address
        santa = table.Santa
        message = "To: $mail\r\nSubject: Weihnachtswichteln $(Dates.year(Dates.today()))\r\n\r\nHallo $name,\r\n\r\nDu darfst heuer $santa bewichteln.\r\n\r\nLG\r\nDein Weihnachtswichtelskript\r\n"
        run(pipeline(`echo $message`, `msmtp -a $mailaccount $mail`))
    end
end

"""
    runsecretsanta(f, a)

Runs the secret santa script using information in file `f`.
Mailing requires a set-up msmtp. Choose your prefered mail account `a`.
"""
function runsecretsanta(filename, mailaccount)
    table = readtable(filename)
    table.Santa = getsecretsanta(table[:,1], table[:,2])
    mailsecretsanta(table, mailaccount)
end

runsecretsanta(ARGS[1], ARGS[2])


text:
module Source
using  ..MacroDatabase
using  ..MacroExpander
using  ..ShellInterop
export source

const Included=Set{String}()
const Compiled=Set{String}()

const COMPILED_PATH = "$FXCRD_PATH/compiled"

function compiled_path(file)
  return "$(COMPILED_PATH)/$(replace(file, "/", "%%"))"
end

function outdated(file)
  return mtime(file) > mtime(compiled_path(file))
end

function source(files...; force=false)
  for file in files
    if force || !(file in Included)
      push!(Included, file)
      if !(file in Compiled) && (force || outdated(file))
        compiled = compile_io(file)
        open(compiled_path(file), "w") do cfile
          println(cfile, compiled)
        end
      end
    end
  end
end

function source_force(ARGS...; input="")
  source(ARGS..., force=true)
  return ""
end

MacroDatabase.Macros["@source"]=Macro("@source", source_force, true, true)

end

max_stars_repo_name: fcard/Fxcrd

max_issues_repo_name: fcard/Fxcrd

max_forks_repo_name: fcard/Fxcrd


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# ReverseMcCormick.jl
# A Reverse McCormick operator library in Julia
# See https://github.com/PSORLab/ReverseMcCormick.jl
#############################################################################
# src/reverse_operators/arithmetic.jl
# Contains definitions of reverse +, -, /, *, promotions, conversion, one,
# zero.
#############################################################################

"""
$(SIGNATURES)

Creates reverse McCormick contractor for `a` = `real(b)`.
"""
function real_rev(a::MC, b::MC)
    isempty(a) && (return a, a)
    b = b ∩ a
    a, b
end

"""
$(SIGNATURES)

Creates reverse McCormick contractor for `a` = `one(b)`.
"""
function one_rev(a::MC, b::MC)
    isempty(a) && (return a, a)
    1.0 ∉ a    && (return a, empty(b))
    a, b
end

"""
$(SIGNATURES)

Creates reverse McCormick contractor for `a` = `zero(b)`.
"""
function zero_rev(a::MC, b::MC)
    isempty(a) && (return a, a)
    0.0 ∉ a    && (return a, empty(b))
    a, b
end

const RAD2DEG = Interval{Float64}(pi)/Interval{Float64}(180.0)
"""
$(SIGNATURES)

Creates reverse McCormick contractor for `a` = `rad2deg(b)`. That is, it is the reverse
contractor of `a = (180/pi)*b`
"""
function rad2deg_rev(a::MC, b::MC)
    isempty(a) && (return a, a)
    b = RAD2DEG*a ∩ b
    a, b
end

const DEG2RAD = Interval{Float64}(180.0)/Interval{Float64}(pi)

"""
$(SIGNATURES)

Creates reverse McCormick contractor for `a` = `deg2rad(b)`. That is, it is the reverse
contractor of `a = (pi/180)*b`
"""
function deg2rad_rev(a::MC, b::MC)
    isempty(a) && (return a, a)
    b = DEG2RAD*a ∩ b
    a, b
end

"""
$(SIGNATURES)

Creates reverse McCormick contractor for `a` = `b` +`c`
"""
function plus_rev(a::MC, b::MC, c::MC)
    isempty(a) && (return a, a, a)
    b_new = b ∩ (a - c)
    c_new = c ∩ (a - b)
    a, b_new, c_new
end
function plus_rev(a::MC{N,T}, b::MC{N,T}, c::C) where {N, T<:RelaxTag, C<:NumberNotRelax}
    isempty(a) && (return a, a, a)
    if c ∈ a - b
        b = b ∩ (a - c)
        return a, b, MC{N,T}(c)
    end
    empty(a), empty(a), empty(a)
end
function plus_rev(a::MC{N,T}, c::C, b::MC{N,T}) where {N, T<:RelaxTag, C<:NumberNotRelax}
    a, b, c = plus_rev(a, b, c)
    a, c, b
end
function plus_rev(a::MC, b::MC)
    isempty(a) && (return a, a)
    b = b ∩ a
    a, b
end

"""
$(SIGNATURES)

Creates reverse McCormick contractor for `a` = `b`- `c`
"""
function minus_rev(a::MC, b::MC, c::MC)
    isempty(a) && (return a, a, a)
    b_new = b ∩ (a + c)
    c_new = c ∩ (b - a)
    a, b_new, c_new
end
function minus_rev(a::MC{N,T}, b::MC{N,T}, c::C) where {N, T<:RelaxTag, C<:NumberNotRelax}
    isempty(a) && (return a, a, a)
    b_new = b ∩ (a + c)
    a, b_new, MC{N,T}(c)
end

function minus_rev(a::MC{N,T}, b::C, c::MC{N,T}) where {N, T<:RelaxTag, C<:NumberNotRelax}
    isempty(a) && (return a, a, a)
    c_new = c ∩ (b - a)
    a, MC{N,T}(b), c_new
end

"""
$(SIGNATURES)

Creates reverse McCormick contractor for `a` = `-b``
"""
function minus_rev(a::MC, b::MC)
    isempty(a) && (return a, a)
    b_new = b ∩ (-a)
    return a, b_new
end

"""
$(SIGNATURES)

Creates reverse McCormick contractor for `a` = `b`*`c`
"""
function mult_rev(a::MC, b::MC, c::MC)
    isempty(a) && (return a, a, a)
    c_new = (0.0 ∉ b) ? (c ∩ (a/b)) : c
    b_new = (0.0 ∉ c) ? (b ∩ (a/c)) : b
    a, b_new, c_new
end

function mult_rev(a::MC{N,T}, b::MC{N,T}, c::C) where {N, T<:RelaxTag, C<:NumberNotRelax}
    isempty(a) && (return a, a, a)
    if !iszero(c)
        b = b ∩ (a*inv(c))
    elseif 0.0 ∉ a
        return empty(a), empty(a), empty(a)
    end
    a, b, MC{N,T}(c)
end

function mult_rev(a::MC{N,T}, c::C, b::MC{N,T}) where {N, T<:RelaxTag, C<:NumberNotRelax}
    if !iszero(c)
        b = b ∩ (a*inv(c))
    elseif 0.0 ∉ a
        return empty(a), empty(a), empty(a)
    end
    a, MC{N,T}(c), b
end

"""
$(SIGNATURES)

Creates reverse McCormick contractor for `a` = `b`/`c`
"""
function div_rev(a::MC, b::MC, c::MC)
    isempty(a) && (return a, a, a)
    b = b ∩ (a * c)
    if ~isempty(b) && 0.0 ∉ a
        c = c ∩ (b / a)
    end
    a, b, c
end

function div_rev(a::MC{N,T}, b::MC{N,T}, c::C) where {N, T<:RelaxTag, C<:NumberNotRelax}
    isempty(a) && (return a, a, a)
    b = b ∩ (a * c)
    a, b, c
end

function div_rev(a::MC{N,T}, b::C, c::MC{N,T}) where {N, T<:RelaxTag, C<:NumberNotRelax}
    isempty(a) && (return a, a, a)
    b ∉ (a * c) && (return a, a, a)
    if 0.0 ∉ a
        c = c ∩ (b / a)
    end
    a, b, c
end

"""
$(SIGNATURES)

Creates reverse McCormick contractor for `a` = `inv(b)`
"""
function inv_rev(a::MC, b::MC)
    isempty(a) && (return a, a)
    0.0 ∉ a && (b = b ∩ inv(a))
    a, b
end

"""
$(SIGNATURES)

Creates reverse McCormick contractor for `a` = `b`^`c`
"""
function power_rev(a::MC, b::MC, c::MC)
    isempty(a) && (return a, a, a)
    (isempty(b) || isempty(c)) && (return a, a, a)
    0.0 ∉ c && (b = b ∩ a^inv(c))
    if 0.0 ∉ a
        if 0.0 < lo(b) < Inf
            blog = log(b)
            if 0.0 ∉ blog
                c = c ∩ (log(a) / blog)
            end
        end
    end
    a, b, c
end
function power_rev(a::MC{N,T}, b::C, c::MC{N,T}) where {N, T<:RelaxTag, C<:NumberNotRelax}
    isempty(a) && (return a, a, a)
    if !isempty(c)
        if (0.0 ∉ a) && (b > zero(C))
            blog = log(b)
            if 0.0 ∉ blog
                c = c ∩ (log(a)/blog)
            end
        end
    else
        a = empty(a)
    end
    a, MC{N,T}(b), c
end

function int_power_rev(a::MC{N,T}, b::MC{N,T}, n::Int, c::C) where {N, T<:RelaxTag, C<:NumberNotRelax}
    if n == 2
        root = sqrt(a)
        if lo(b) > 0.0
            b = b ∩ root
        else
            b = b ∩ -root
        end
    elseif iseven(n)
        root = a^(1/n)
        if lo(b) > 0.0
            b = b ∩ root
        else
            b = b ∩ -root
        end
    elseif isodd(n)
        if lo(b) > 0.0
            b = b ∩ (a ∩ Interval(0, Inf))^(1/n)
        else
            b = b ∩ -(((-a) ∩ Interval(0, Inf))^(1/n))
        end
    end
    a, b, MC{N,T}(c)
end

function power_rev(a::MC{N,T}, b::MC{N,T}, c::C) where {N, T<:RelaxTag, C<:NumberNotRelax}
    isempty(a) && (return a, a, a)
    (isempty(b) || 0.0 > lo(b)) && (return a, a, a)
    if isone(-c)
        a, b = inv_rev(a,b)
        return a, b, MC{N,T}(c)
    elseif !iszero(c)
        if isinteger(c)
            return int_power_rev(a, b, Int(c), c)
        elseif lo(b) > 0.0
            b = b ∩ a^(one(C)/c)
        end
    elseif iszero(c) && 1.0 ∉ a
        return empty(a), empty(a), empty(a)
    end
    a, b, MC{N,T}(c)
end

"""
$(SIGNATURES)

Creates reverse McCormick contractor for `a` = `sqrt(b)`. That is
`b = b ∩ a^2`
"""
function sqrt_rev(a::MC, b::MC)
    isempty(a) && (return a, a)
    b = b ∩ (a^2)
    a, b
end
sqr_rev(f, x) = power_rev(f, x, 2)

"""
$(SIGNATURES)

Creates reverse McCormick contractor for `a` = `abs(b)`
"""
function abs_rev(y::MC, x::MC)
    isempty(y) && (return y, y)
    y, x
end

"""
$(SIGNATURES)

Creates reverse McCormick contractor for `a` = `step(b)`
"""
function step_rev(a::MC, b::MC)
    a_lo = a.Intv.lo
    a_hi = a.Intv.hi
    if isempty(a) || 0.0 > a_hi || 1.0 < a_lo ||
        (0.0 < a_lo && a_hi < 1.0)
        b = empty(b)
        return a, b
    elseif 0.0 ∈ a && 1.0 ∉ a
        b = b ∩ Interval(-Inf, 0.0)
    elseif -1.0 ∉ a && 1.0 ∈ a
        b = b ∩ Interval(0.0, Inf)
    end
    a, b
end

"""
$(SIGNATURES)

Creates reverse McCormick contractor for `a` = `sign(b)`.
"""
function sign_rev(a::MC, b::MC)
    a_lo = a.Intv.lo
    a_hi = a.Intv.hi
    if isempty(a) || -1.0 > a_hi || 1.0 < a_lo || (-1.0 < a_lo && a_hi < 1.0)
        b = empty(b)
        return a, b
    elseif -1.0 ∈ a && 1.0 ∉ a
        b = b ∩ Interval(-Inf, 0.0)
    elseif -1.0 ∉ a && 1.0 ∈ a
        b = b ∩ Interval(0.0, Inf)
    end
    a, b
end


text:
using LombScargle
using Measurements, FFTW, StableRNGs
using Test

Threads.nthreads() > 1 && (FFTW.set_num_threads(2); @info("Multi-threading enabled"))

ntimes = 1001
# Observation times
t = range(0.01, stop=10pi, length=ntimes)
# Signal
s = sin.(t)

pgram1 = @inferred(lombscargle(LombScargle.plan(t, s, fast = false, fit_mean=false)))
pgram2 = @inferred(lombscargle(LombScargle.plan(t, s, fast = false, fit_mean=true)))
pgram3 = @inferred(lombscargle(LombScargle.plan(t, s, fast = false, center_data=false, fit_mean=false)))
pgram4 = @inferred(lombscargle(LombScargle.plan(t, s, fast = false, center_data=false, fit_mean=true)))

@testset "Periodograms" begin
    @testset "Random stuff" begin
        rng = StableRNG(1)
        # Randomized times
        trandom = t .+ step(t) .* rand(rng, ntimes)
        # Randomized signal
        srandom = sinpi.(trandom) .+ cospi.(2trandom) .+ rand(rng, ntimes)
        # Frequency grid
        nfreqs = 10000
        freqs = range(0.01, stop=3, length=nfreqs)
        # Randomize frequency grid
        freqs += step(freqs) * rand(rng, nfreqs)
        # Use "freqpower" and "periodpower" just to call that function and increase code
        # coverage.  "autofrequency" function is tested below.
        prandom1 = lombscargle(trandom, srandom, frequencies=freqs, fit_mean=false)
        prandom2 = lombscargle(trandom, srandom, frequencies=freqs, fit_mean=true)
        @test freqpower(pgram1)[2] ≈ periodpower(pgram2)[2] atol = 6e-7
        @testset "Infinities" begin
            # Make sure there are no infinities in `pgram1'.  It seems to work only on
            # 64-bit systems.
            Sys.WORD_SIZE == 64 && @test(!(Inf in power(prandom1)) && !(Inf in power(prandom2)))
        end
    end

    @testset "fit_mean and center_data" begin
        @test power(pgram1) ≈ power(pgram2) atol = 5e-7
        @test power(pgram3) ≈ power(pgram4) atol = 5e-7
    end

    # Test the values in order to prevent wrong results in both algorithms
    @testset "Test values" begin
        @test power(lombscargle(t, s, frequencies=0.2:0.2:1, fit_mean=true)) ≈ [0.029886871262324886,0.0005456198989410226,1.912507742056023e-5, 4.54258409531214e-6, 1.0238342782430832e-5]
        @test power(lombscargle(t, s, frequencies=0.2:0.2:1, fit_mean=false)) ≈ [0.02988686776042212, 0.0005456197937194695,1.9125076826683257e-5,4.542583863304549e-6,1.0238340733199874e-5]
        @test power(lombscargle(t, s, frequencies=0.2:0.2:1, center_data=false, fit_mean=true)) ≈ [0.029886871262325004,0.0005456198989536703,1.9125077421448458e-5,4.5425840956285145e-6,1.023834278337881e-5]
        @test power(lombscargle(t, s, frequencies=0.2:0.2:1, center_data=false, fit_mean=false)) ≈ [0.029886868328967767,0.0005456198924872134,1.9125084251687147e-5,4.542588504467314e-6,1.0238354525870936e-5]
        @test power(lombscargle(t, s, frequencies=0.2:0.2:1, normalization=:model)) ≈ [0.030807614469885718,0.0005459177625354441,1.9125443196143085e-5,4.54260473047638e-6,1.0238447607164715e-5]
        @test power(lombscargle(t, s, frequencies=0.2:0.2:1, normalization=:log)) ≈ [0.030342586720560734,0.0005457688036440774,1.9125260307148152e-5,4.542594412890309e-6,1.0238395194654036e-5]
        @test power(lombscargle(t, s, frequencies=0.2:0.2:1, normalization=:psd)) ≈ [7.474096700871138,0.1364484040771917,0.004782791641128195,0.0011360075968541799,0.002560400630125523]
        @test power(lombscargle(t, s, frequencies=0.2:0.2:1, normalization=:Scargle)) ≈ [0.029886871262324904,0.0005456198989410194,1.912507742056126e-5,4.54258409531238e-6,1.0238342782428552e-5]
        @test power(lombscargle(t, s, frequencies=0.2:0.2:1, normalization=:HorneBaliunas)) ≈ [14.943435631162451,0.2728099494705097,0.009562538710280628,0.00227129204765619,0.005119171391214276]
        @test power(lombscargle(t, s, frequencies=0.2:0.2:1, normalization=:Cumming)) ≈ [15.372999620472974,0.2806521440709115,0.009837423440787873,0.0023365826071340815,0.005266327088140394]
        @test_throws ErrorException lombscargle(t, s, frequencies=0.2:0.2:1, normalization=:foo)
    end

    err = collect(range(0.5, stop=1.5, length=ntimes))

    @testset "Signal with uncertainties" begin
        @test power(lombscargle(t, s, err, frequencies=0.1:0.1:1, fit_mean=true)) ≈ [0.06659683848818691,0.09230959166317589,0.006625919314669043,0.0015664010997692612,0.0005085442118408477,0.00019704659245878378,9.658452525613897e-5,6.331573873913433e-5,4.903871967643573e-5,3.7948448825374025e-5]
        @test power(lombscargle(t, s, err, frequencies=0.1:0.1:1, fit_mean=false)) ≈ [0.0664002483305464,0.09219168665786254,0.006625915010614472,0.0015663421089042564,0.0005085109569237008,0.00019703233981948823,9.6577091433651e-5,6.33101344670203e-5,4.9033581990442793e-5,3.7944076990210425e-5]
        @test power(lombscargle(t, s, err, frequencies=0.1:0.1:1, fit_mean=false, center_data=false)) ≈ [0.06920814049261209,0.09360344864985352,0.006634919960009565,0.0015362072871144769,0.0004858250831632676,0.00018179850370583626,8.543727416919218e-5,5.379994730581837e-5,4.0107232867763e-5,2.9784059487535237e-5]
        @test power(lombscargle(t, s, err)) ==
            power(lombscargle(t, measurement.(s, err)))
        @test power(lombscargle(t, s, err, frequencies = 0.1:0.1:1, fast = false, normalization = :psd)) ≈ [21.851224476672318,30.287888352835566,2.1740438975167593,0.5139550589572747,0.16685947834022155,0.06465335925734642,0.031690545531213095,0.020774656147387098,0.01609019430987704,0.012451342926314715]
    end
#=
    pgram5 = lombscargle(t, s, maximum_frequency=30, fast=true)
    pgram6 = lombscargle(t, s, maximum_frequency=30, fast=false)

    @testset "Fast method" begin
        @test power(pgram5) ≈ power(pgram6) atol = 3e-4
        @test power(lombscargle(t, s, frequencies=0.2:0.2:1, fast=true, fit_mean=true, flags = FFTW.MEASURE, timelimit = 5)) ≈
            [0.029886871262325053, 0.0005447325913220627, 2.1246300058375996e-5, 4.1259517049745417e-7, 5.04610747916143e-5]
        @test power(lombscargle(t, s, frequencies=0.2:0.2:1, fast=true, fit_mean=false, flags = FFTW.MEASURE, timelimit = 5)) ≈
            [0.0298868677604221, 0.0005447324642349108, 2.1246261449681898e-5, 4.125675324400738e-7, 5.042481545188997e-5]
        @test power(lombscargle(t, s, frequencies=0.2:0.2:1, fast=true, center_data=false, fit_mean=false, flags = FFTW.MEASURE, timelimit = 5)) ≈
            [0.029886868328967798, 0.0005447325727792022, 2.1246201600761156e-5, 4.1251689656274853e-7, 5.0422982073668846e-5]
        @test power(lombscargle(t, s, err, frequencies=0.2:0.2:1, fast=true, fit_mean=true)) ≈ [0.09230959166317658, 0.0015636510410697796, 0.00019210902295493893, 8.221668511639927e-5, 0.00021747947231895047]
        @test power(lombscargle(t, s, err, frequencies=0.2:0.2:1, fast=true, fit_mean=false)) ≈ [0.09219168665786256, 0.0015635926899499144, 0.00019209076832480172, 8.2184576575774e-5, 0.00021723818031507325]
        @test power(lombscargle(t, s, err, frequencies=0.2:0.2:1, fast=true, center_data=false, fit_mean=false)) ≈ [0.09360344864985316, 0.001533862558849818, 0.00017415674032714024, 5.1945371781883936e-5, 0.00011977935090727627]
        @test power(lombscargle(t, s, err)) ==
            power(lombscargle(t, measurement.(s, err)))
    end

    # Compare result of uncertainties with both methods (fast and non-fast).
    @testset "Fast and non-fast methods" begin
        rng = StableRNG(1)
        errors = rand(rng, 0.1:1e-3:4.0, ntimes)
        @test power(lombscargle(t, s, errors)) ≈ power(lombscargle(t, s, errors, fast = false)) atol = 0.3
        @test power(lombscargle(t, s, errors, fit_mean = false)) ≈ power(lombscargle(t, s, errors, fit_mean = false, fast = false)) atol = 0.2
    end

   @testset "Float32" begin
       @test typeof(period(lombscargle(Vector{Float32}(t), Vector{Float32}(s), fast = false))) == Vector{Float32}
   end
end

@testset "Utilities" begin
    @testset "findmaxpower, findmaxfreq, findmaxperiod" begin
        @test findmaxfreq(pgram1) ≈ [31.997145470342]
        @test findmaxfreq(pgram1, 0.965) ≈ [0.15602150741832602,31.685102455505348,31.997145470342,63.52622641842902,63.838269433265665]
        @test findmaxperiod(pgram1) ≈ 1 ./ findmaxfreq(pgram1)
        @test findmaxperiod(pgram1, 0.965) ≈ 1 ./ findmaxfreq(pgram1, 0.965)
        global s = sinpi.(2t) .+ cospi.(4t)
        p = lombscargle(t, s, maximum_frequency=4)
        @test findmaxfreq(p, [0.9, 1.1]) ≈ [1.0029954048320957]
        @test findmaxfreq(p, [1.9, 2.1]) ≈ [2.002806697267899]
        @test findmaxperiod(p, [1/0.9, 1/1.1]) ≈ 1 ./ findmaxfreq(p, [0.9, 1.1])
        @test findmaxperiod(p, [1/1.9, 1/2.1]) ≈ 1 ./ findmaxfreq(p, [1.9, 2.1])
        @test findmaxpower(pgram1) ≈ 0.9695017551608017
    end

    @testset "LombScargle.autofrequency" begin
        @test LombScargle.autofrequency(t) ≈ 0.003184112396292367:0.006368224792584734:79.6824127172165
        @test LombScargle.autofrequency(t, minimum_frequency=0) ≈ 0.0:0.006368224792584734:79.6792286048202
        @test LombScargle.autofrequency(t, maximum_frequency=10) ≈ 0.003184112396292367:0.006368224792584734:9.99492881196174
        # This last test also makes sure that `freq` and `power` fields of a Periodogram
        # object can have different type.
        if Threads.nthreads() > 1 && Sys.isapple()
            # This would make Julia crash, skip it.  See
            # https://github.com/JuliaLang/julia/issues/35702
            @test_skip freq(lombscargle(1:11, big.(sin.(1:11)))) ≈ 0.01:0.02:2.75
        else
            @test freq(lombscargle(1:11, big.(sin.(1:11)))) ≈ 0.01:0.02:2.75
        end
    end

    @testset "Probabilities and FAP" begin
        global t = collect(range(0.01, stop = 10pi, length = 101))
        global s = sin.(t)
        for norm in (:standard, :Scargle, :HorneBaliunas, :Cumming)
            P = lombscargle(t, s, normalization = norm)
            for z_0 in (0.1, 0.5, 0.9)
                @test prob(P, probinv(P, z_0)) ≈ z_0
                @test fap(P,  fapinv(P, z_0)) ≈ z_0
            end
        end
        P = lombscargle(t, s, normalization = :log)
        @test_throws ErrorException prob(P, 0.5)
        @test_throws ErrorException probinv(P, 0.5)
    end

    @testset "LombScargle.model" begin
        @test s ≈ LombScargle.model(t, s, 1/2pi, center_data=false, fit_mean=false)
        global s = sinpi.(t) .+ pi .* cospi.(t) .+ ℯ
        @test s ≈ LombScargle.model(t, measurement.(s, fill(1, size(s))), 0.5)
    end

    @testset "LombScargle.add_at!" begin
        a = ones(Int, 3)
        LombScargle.add_at!(a, [3, 1, 3, 1, 2], 1:5)
        @test a ≈ [7, 6, 5]
    end

    @testset "LombScargle.extirpolate!" begin
        x = collect(range(0, stop = 10, length = 50))
        y = sin.(x)
        vec13 = Vector{Complex{Float64}}(undef, 13)
        LombScargle.extirpolate!(vec13, x, y, 13)
        @test vec13 ≈ [0.39537718210649553,3.979484140636793,4.833090108345013,0.506805556164743,-3.828112427525919,-4.748341359084166,-1.3022050566901917,3.3367666084342256,5.070478111668922,1.291245296032218,-0.8264466821981216,0.0,0.0]
        x = collect(range(0, stop = 10, length = 50))
        y = sin.(x)
        @test LombScargle.extirpolate!(Vector{Complex{Float64}}(undef, 11), x, y, 11) ≈ vec13[1:11]
    end

    x = collect(range(0, stop = 10, length = 50))
    y = sin.(x)

    @testset "LombScargle.trig_sum!" begin
        N = 10
        Nfft = nextpow(2, 5N)
        fftgrid = Vector{Complex{Float64}}(undef, N)
        bfft_vec = Vector{Complex{Float64}}(undef, Nfft)
        p = plan_bfft(bfft_vec)
        grid = similar(bfft_vec)
        C, S = LombScargle.trig_sum!(grid, fftgrid, bfft_vec, p, x, y, 1, N, Nfft, minimum(x))
        @test S ≈ [0.0,0.3753570125888358,0.08163980192703546,-0.10139634351774979,-0.4334223744905633,-2.7843373311769777,0.32405810159838055,0.05729663600471602,-0.13191736591325876,-0.5295781583202946]
        @test C ≈ [8.708141477890015,-0.5402668064176129,-0.37460815054027985,-0.3793457539084364,-0.5972351546196192,14.612204307982497,-0.5020253140297526,-0.37724493022381034,-0.394096831764578,-0.6828241623474718]
    end

    @testset "Bootstrap" begin
        rng = StableRNG(1)
        @test LombScargle.bootstrap(rng, 5, x, y).p ≈
            [0.2583163570869385, 0.24972129609003385, 0.23092196417031927, 0.18502993723773883, 0.1789661587851332]
        rng = StableRNG(1)
        @test LombScargle.bootstrap(rng, 5, x, y, fit_mean = false).p ≈
            [0.2580212594813987, 0.24897444742007394, 0.23090538831280463, 0.18492853793841382, 0.17895144706266247]
        err = collect(range(0.5, stop = 1.5, length = 50))
        rng = StableRNG(1)
        b = LombScargle.bootstrap(rng, 50, x, measurement.(y, err), fast = true)
        @test fap(b, fapinv(b, 0.02)) ≈ 0.02
        rng = StableRNG(1)
        b = LombScargle.bootstrap(rng, 50, x, measurement.(y, err), fast = false)
        @test fap(b, fapinv(b, 0.02)) ≈ 0.02
        rng = StableRNG(1)
        @test fapinv(LombScargle.bootstrap(rng, 50, x, measurement.(y, err), fast = false, fit_mean = false),
                     0.2) ≈ 0.23917691901908134
        rng = StableRNG(1)
        @test fapinv(LombScargle.bootstrap(rng, 1000, x, y, fast=false, fit_mean=false), 0.2) ≈
            0.2157617143004672
    end
=#
end


text:
module DiscreteHO
    include("Krav.jl")
end


text:
using Reactive
using Test
# Stop the runner task

@testset "Queue runner" begin
    @test (istaskdone(Reactive.runner_task[])) == (false)
    Reactive.stop()
    @test (istaskdone(Reactive.runner_task[])) == (true)
end

step() = Reactive.run(1)
queue_size() = Base.n_avail(Reactive._messages)
number() = round(Int, rand()*1000)

include("basics.jl")
include("push_to_non_input.jl")
# include("gc.jl")
include("node_order.jl")
include("call_count.jl")
include("flatten.jl")
include("time.jl")
include("async.jl")
include("queue_runner.jl")


text:
function print_block(fun::Function,out=stdout)
    io = IOBuffer()
    fun(io)
    output = String(take!(io))
    isempty(output) && return
    data = split(output, "\n")
    isempty(data[end]) && (data = data[1:end-1])
    if length(data) == 1
        println(out, "[ $(data[1])")
    elseif length(data) > 1
        for (p,dl) in zip(vcat("╭", repeat(["│"], length(data)-2), "╰"),data)
            println(out, "$(p) $(dl)")
        end
    end
    println(out)
end


text:
import Distributions

"""
    perfectNormal(n::Int, mean::Number=0, sd::Number=1)

Generate an almost-perfect normal distribution of size `n`.

# Arguments
- `n::Int`: Length of the vector.
- `mean::Number`: Mean of the vector.
- `sd::Number`: SD of the vector.

# Examples
```julia
perfectNormal(10, 0, 1)
```
"""
function perfectNormal(n::Int, mean::Number=0, sd::Number=1)
    x = Distributions.quantile.([Distributions.Normal(mean, sd)],
                                range(1/n, stop=1-1/n, length=n))
    return x
end
perfectNormal(; n::Int, mean::Number=0, sd::Number=1) = perfectNormal(n, mean, sd)


text:
module FourierTools

using PaddedViews, ShiftedArrays
using FFTW
FFTW.set_num_threads(4)
using LinearAlgebra
using IndexFunArrays
using ChainRulesCore
using NDTools

include("utils.jl")
include("resampling.jl")
include("custom_fourier_types.jl")
include("fourier_resizing.jl")
include("fourier_shifting.jl")
include("fourier_resample_1D_based.jl")
include("fourier_rotate.jl")
include("fourier_shear.jl")
include("fftshift_alternatives.jl")
include("fft_helpers.jl")
include("convolutions.jl")
include("damping.jl")
include("czt.jl")

end # module


text:
function decode_msg(line::AbstractString)
    pos = our_findfirst(':', line)
    pos2 = our_findnext(':', line, pos + 1)

    msg_id = line[1:pos - 1]
    msg_cmd = line[pos + 1:pos2 - 1]
    msg_body_encoded = line[pos2 + 1:end]
    msg_body = String(Base64.base64decode(msg_body_encoded))
    return msg_id, msg_cmd, msg_body
end

function send_notification(conn, msg_cmd::AbstractString, msg_body::AbstractString = "")
    encoded_msg_body = Base64.base64encode(msg_body)
    println(conn, msg_cmd, ":notification:", encoded_msg_body)
end

function send_response(conn, msg_id::AbstractString, msg_body::AbstractString)
    encoded_msg_body = Base64.base64encode(msg_body)
    println(conn, "RESPONSE:", msg_id, ':', encoded_msg_body)
end

max_stars_repo_name: cloudgiant/julia-vscode

max_issues_repo_name: cloudgiant/julia-vscode

max_forks_repo_name: cloudgiant/julia-vscode