text:
abstract type ODESolver end

abstract type RK4Solver{nsteps} <: ODESolver  end
abstract type OutOfPlaceRK4Solver{nsteps} <: ODESolver  end

@doc doc"""

    RK4Solver(F!::Function, y₀, t₀, t₁, nsteps)

Solve for $y(t_1)$ with 4th order Runge-Kutta assuming $dy/dt = F(t,y)$ and $y(t_0)$ = $y_0$.

Arguments:
* `F!` — a function `F!(v,t,y)`` which sets `v=F(t,y)`
"""
function RK4Solver(F!::Function, y₀, t₀, t₁, nsteps)
    h, h½, h⅙ = (t₁-t₀)/nsteps ./ (1,2,6)
    y = copy(y₀)
    k₁, k₂, k₃, k₄, y′ = @repeated(similar(y₀),5)
    for t in range(t₀,t₁,length=nsteps+1)[1:end-1]
        @! k₁ = F(t, y)
        @! k₂ = F(t + h½, (@. y′ = y + h½*k₁))
        @! k₃ = F(t + h½, (@. y′ = y + h½*k₂))
        @! k₄ = F(t + h,  (@. y′ = y + h*k₃))

        # due to https://github.com/JuliaLang/julia/issues/27988, if this were
        # written the natural way as:
        #    @. y .+= h*(k₁ + 2k₂ + 2k₃ + k₄)/6
        # it has god-awful performance for FieldTuples (although is fine for
        # FlatS0s). until a solution for that issue comes around, a workaround
        # is to write out the broadcasting kernel by hand:
        broadcast!((y,h,k₁,k₂,k₃,k₄)->(y+h⅙*(k₁+2*(k₂+k₃)+k₄)), y, (y,h,k₁,k₂,k₃,k₄)...)
    end
    return y
end

function OutOfPlaceRK4Solver(F::Function, y₀, t₀, t₁, nsteps)
    h, h½, h⅙ = (t₁-t₀)/nsteps ./ (1,2,6)
    y = copy(y₀)
    for i in 0:nsteps-1
        t = i*h
        k₁ = F(t, y)
        k₂ = F(t + h½, y + h½*k₁)
        k₃ = F(t + h½, y + h½*k₂)
        k₄ = F(t + h,  y + h*k₃)
        y += @. h*(k₁ + 2k₂ + 2k₃ + k₄)/6
    end
    return y
end

odesolve(::Type{RK4Solver{N}},F!,y₀,t₀,t₁) where {N} = RK4Solver(F!,y₀,t₀,t₁,N)
odesolve(::Type{OutOfPlaceRK4Solver{N}},F,y₀,t₀,t₁) where {N} = OutOfPlaceRK4Solver(F,y₀,t₀,t₁,N)

@doc doc"""
    conjugate_gradient(M, A, b, x=M\b; nsteps=length(b), tol=sqrt(eps()), progress=false, callback=nothing, hist=nothing, histmod=1)

Compute `x=A\b` (where `A` is positive definite) by conjugate gradient. `M` is the
preconditioner and should be `M≈A`, and `M\x` should be fast.

The solver will stop either after `nsteps` iterations or when `dot(r,r)<tol`
(where `r=A*x-b` is the residual  at that step), whichever occurs first.

Info from the iterations of the solver can be returned if `hist` is specified.
`hist` can be one or a tuple of:

* `:i` — current iteration number
* `:x` — current solution
* `:r` — current residual `r=A*x-b`
* `:res` — the norm of `r`
* `:t` — the time elapsed (in seconds) since the start of the algorithm

`histmod` can be used to include every N-th iteration only in `hist`.
"""
function conjugate_gradient(M, A, b, x=0*b; nsteps=length(b), tol=sqrt(eps()), progress=false, callback=nothing, hist=nothing, histmod=1)
    gethist() = hist == nothing ? nothing : NamedTuple{hist}(getindex.(Ref(@dict(i,x,p,r,res,t)),hist))
    t₀ = time()
    i = 1
    r = b - A*x
    z = M \ r
    p = z
    bestres = res = res₀ = dot(r,z)
    @assert !isnan(res)
    bestx = x
    t    = time() - t₀
    _hist = [gethist()]

    prog = Progress(100, (progress!=false ? progress : Inf), "Conjugate Gradient: ")
    for outer i = 2:nsteps
        Ap   = A * p
        α    = res / dot(p,Ap)
        x    = x + α * p
        r    = r - α * Ap
        z    = M \ r
        res′ = dot(r,z)
        p    = z + (res′ / res) * p
        res  = res′
        t    = time() - t₀

        if all(res<bestres)
            bestres,bestx = res,x
        end
        if callback!=nothing
            callback(i,x,res)
        end
        if hist!=nothing && (i%histmod)==0
            push!(_hist, gethist())
        end
        if all(res<tol)
            break
        end

        # update progress bar to whichever we've made the most progress on,
        # logarithmically reaching the toleranace limit or doing the maximum
        # number of steps
        if progress
            progress_nsteps = round(Int,100*(i-1)/(nsteps-1))
            progress_tol = res isa AbstractFloat ? round(Int,100^min(1, (log10(res/res₀)) / log10(tol/res₀))) : 0
            ProgressMeter.update!(prog, max(progress_nsteps,progress_tol))
        end
    end
    ProgressMeter.finish!(prog)
    hist == nothing ? bestx : (bestx, _hist)
end

@doc doc"""

    itp = LinearInterpolation(xdat::AbstractVector, ydat::AbstractVector; extrapolation_bc=NaN)
    itp(x) # interpolate at x

A simple 1D linear interpolation code which is fully Zygote differentiable in
either `xdat`, `ydat`, or the evaluation point `x`.
"""
function LinearInterpolation(xdat::AbstractVector, ydat::AbstractVector; extrapolation_bc=NaN)

    @assert issorted(xdat)
    @assert extrapolation_bc isa Number || extrapolation_bc == :line

    m = diff(ydat) ./ diff(xdat)

    function (x::Number)

        if x<xdat[1] || x>xdat[end]
            if extrapolation_bc isa Number
                return extrapolation_bc
            elseif extrapolation_bc == :line
                if x<xdat[1]
                    i = 1
                elseif x>xdat[end]
                    i = length(m)
                end
            end
        else
            # sets i such that x is between xdat[i] and xdat[i+1]
            i = max(1, searchsortedfirst(xdat,x) - 1)
        end

        # do interpolation
        ydat[i] + m[i]*(x-xdat[i])

    end

end

@doc doc"""

    gmres(A, b; maxiter, Pl=I)

Solve `A \ b` with `maxiter` iterations of the [generalized minimal
residual](https://en.wikipedia.org/wiki/Generalized_minimal_residual_method)
algorithm. `Pl` is a left-preconditioner which should approximate `inv(A)`.

Note: the implemenation is memory inefficient and uses O(n * maxiter) memory, where
`n,n=size(A)` (may not be a big deal for small `maxiter`), although is totally generic
and works with CPU or GPU and dense or sparse matrices, unlike IterativeSolver's
`gmres`.
"""
function gmres(A, b; Pl=I, maxiter)

    n = maxiter
    T = promote_op(matprod, eltype(A), eltype(b))

    # build Krylov matrix K = [(Pl*A)*b (Pl*A)²*b ...]
    K = similar(b, T, length(b), n+1)
    mul!(view(K, :, 1), Pl, b)
    for i=2:n+1
        mul!(view(K,:,i), A, view(K,:,i-1))
        if Pl !== I
            mul!(view(K,:,i), Pl, K[:,i]) # copy needed here
        end
    end

    # solve least-squares problem |Pl * A * K * α - Pl * b|²
    α = qr(view(K, :, 2:n+1)) \ view(K, :, 1)

    # return solution, K * α
    view(K, :, 1:n) * α

end


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

## generate repl-docs ##

function generate(io, command)
    cmd_nospace = replace(command, " " => "-")
    println(io, """
    ```@raw html
    <article class="docstring">
        <header>
            <a class="docstring-binding" id="repl-$(cmd_nospace)" href="#repl-$(cmd_nospace)">
                <code>$(command)</code>
            </a>
            —
            <span class="docstring-category">REPL command</span>
        </header>
        <section>
    ```
    ```@eval
    using Pkg
    Dict(Pkg.REPLMode.canonical_names())["$(command)"].help
    ```
    ```@raw html
        </section>
    </article>
    ```
    """)
end
function generate()
    io = IOBuffer()
    println(io, """
        # [**11.** REPL Mode Reference](@id REPL-Mode-Reference)

        This section describes available commands in the Pkg REPL.
        The Pkg REPL mode is mostly meant for interactive use,
        and for non-interactive use it is recommended to use the
        functional API, see [API Reference](@ref API-Reference).
        """)
    # list commands
    println(io, "## `package` commands")
    foreach(command -> generate(io, command), ["add", "build", "compat", "develop", "free", "generate", "pin", "remove", "test", "update"])
    println(io, "## `registry` commands")
    foreach(command -> generate(io, command), ["registry add", "registry remove", "registry status", "registry update"])
    println(io, "## Other commands")
    foreach(command -> generate(io, command), ["activate", "gc", "help", "instantiate", "precompile", "resolve", "status"])
    # write to file
    write(joinpath(@__DIR__, "src", "repl.md"), seekstart(io))
    return
end

generate()

ext: jl


text:
@testset "knight.jl" begin
    @printf("> knight\n")

    #corners
    @test Bb.KNIGHT[p2u["a1"]] == reduce(|, p2u[s] for s = ["b3", "c2"])
    @test Bb.KNIGHT[p2u["h1"]] == reduce(|, p2u[s] for s = ["f2", "g3"])
    @test Bb.KNIGHT[p2u["h8"]] == reduce(|, p2u[s] for s = ["f7", "g6"])
    @test Bb.KNIGHT[p2u["a8"]] == reduce(|, p2u[s] for s = ["b6", "c7"])

    #sides
    @test Bb.KNIGHT[p2u["e1"]] == reduce(|, p2u[s] for s = ["d3", "c2", "g2", "f3"])
    @test Bb.KNIGHT[p2u["d8"]] == reduce(|, p2u[s] for s = ["b7", "c6", "f7", "e6"])
    @test Bb.KNIGHT[p2u["b4"]] == reduce(|, p2u[s] for s = ["a6", "c6", "a2",
                                                            "c2", "d5", "d3"])

    #center
    @test Bb.KNIGHT[p2u["d4"]] == reduce(|, p2u[s] for s =  ["c2", "b3", "b5", "c6",
                                                             "e6", "f5", "f3", "e2"])
end


text:
# test phase shifting

T = Float32
fs = 100.  # sampling rate
N = 2^14   # number of points
starttime = d2u(DateTime(Date(now())))

@testset "phase shift" begin

C = SeisChannel()
C.x = rand(T,N)
C.fs = fs
C.t = [1 starttime * 1e6;length(C.x) 0]

# test no shift
Cnew = phase_shift(C,ϕshift=true)
@test Cnew == C

# test phase shift + time shift
C.t = [1 (starttime + 0.007) * 1e6;length(C.x) 0]
Cnew = phase_shift(C,ϕshift=true)
@test Cnew != C
@test Cnew.t != C.t
@test Cnew.x != C.x

# test just time shift
Cnew = phase_shift(C,ϕshift=false)
@test Cnew != C
@test Cnew.t != C.t
@test Cnew.x == C.x

# test in-place
Cnew = phase_shift(C,ϕshift=true)
phase_shift!(C)
@test Cnew == C
@test Cnew.t == C.t
@test Cnew.x == C.x

## test SeisData
S = SeisData(2)
S.x[1] = rand(T,N)
S.fs[1] = fs
# unchanged starttime
S.t[1] = [1 starttime * 1e6;length(C.x) 0]
S.x[2] = rand(T,N)
S.fs[2] = fs
# offset starttime
S.t[2] = [1 (starttime + 0.007) * 1e6;length(C.x) 0]

# test phase shift + time shift
Snew = phase_shift(S,ϕshift=true)
@test Snew != S
@test Snew[1] == S[1]
@test Snew[2] != S[2]

# test just time shift
Snew = phase_shift(S,ϕshift=false)
@test Snew != S
@test Snew.t[1] == S.t[1]
@test Snew.x[1] == S.x[1]
@test Snew.t[2] != S.t[2]
@test Snew.x[2] == S.x[2]

# test in-place
Snew = phase_shift(S,ϕshift=true)
phase_shift!(S)
@test Snew == S
@test Snew.t[1] == S.t[1]
@test Snew.x[1] == S.x[1]
@test Snew.t[2] == S.t[2]
@test Snew.x[2] == S.x[2]
end


text:
export BrickletIndustrialCounterSignalData
struct BrickletIndustrialCounterSignalData
    duty_cycle::Integer
    period::Integer
    frequency::Integer
    value::Bool
end

export BrickletIndustrialCounterAllSignalData
struct BrickletIndustrialCounterAllSignalData
    duty_cycle::Vector{Integer}
    period::Vector{Integer}
    frequency::Vector{Integer}
    value::Vector{Bool}
end

export BrickletIndustrialCounterCounterConfiguration
struct BrickletIndustrialCounterCounterConfiguration
    count_edge::Integer
    count_direction::Integer
    duty_cycle_prescaler::Integer
    frequency_integration_time::Integer
end

export BrickletIndustrialCounterAllCounterCallbackConfiguration
struct BrickletIndustrialCounterAllCounterCallbackConfiguration
    period::Integer
    value_has_to_change::Bool
end

export BrickletIndustrialCounterAllSignalDataCallbackConfiguration
struct BrickletIndustrialCounterAllSignalDataCallbackConfiguration
    period::Integer
    value_has_to_change::Bool
end

export BrickletIndustrialCounterSPITFPErrorCount
struct BrickletIndustrialCounterSPITFPErrorCount
    error_count_ack_checksum::Integer
    error_count_message_checksum::Integer
    error_count_frame::Integer
    error_count_overflow::Integer
end

export BrickletIndustrialCounterIdentity
struct BrickletIndustrialCounterIdentity
    uid::String
    connected_uid::String
    position::Char
    hardware_version::Vector{Integer}
    firmware_version::Vector{Integer}
    device_identifier::Integer
end

export BrickletIndustrialCounter
"""
4 channel counter up to 4MHz
"""
mutable struct BrickletIndustrialCounter <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

    """
    Creates an object with the unique device ID *uid* and adds it to
    the IP Connection *ipcon*.
    """
    function BrickletIndustrialCounter(uid::String, ipcon::IPConnection)
        package = pyimport("tinkerforge.bricklet_industrial_counter")
        deviceInternal = package.BrickletIndustrialCounter(uid, ipcon.ipconInternal)

        return new(ipcon, deviceInternal)
    end
end

export get_counter
"""
    $(SIGNATURES)

Returns the current counter value for the given channel.
"""
function get_counter(device::BrickletIndustrialCounter, channel)
    return device.deviceInternal.get_counter(channel)
end

export get_all_counter
"""
    $(SIGNATURES)

Returns the current counter values for all four channels.
"""
function get_all_counter(device::BrickletIndustrialCounter)
    return device.deviceInternal.get_all_counter()
end

export set_counter
"""
    $(SIGNATURES)

Sets the counter value for the given channel.

The default value for the counters on startup is 0.
"""
function set_counter(device::BrickletIndustrialCounter, channel, counter)
    device.deviceInternal.set_counter(channel, counter)
end

export set_all_counter
"""
    $(SIGNATURES)

Sets the counter values for all four channels.

The default value for the counters on startup is 0.
"""
function set_all_counter(device::BrickletIndustrialCounter, counter)
    device.deviceInternal.set_all_counter(counter)
end

export get_signal_data
"""
    $(SIGNATURES)

Returns the signal data (duty cycle, period, frequency and value) for the
given channel.
"""
function get_signal_data(device::BrickletIndustrialCounter, channel)
    return device.deviceInternal.get_signal_data(channel)
end

export get_all_signal_data
"""
    $(SIGNATURES)

Returns the signal data (duty cycle, period, frequency and value) for all four
channels.
"""
function get_all_signal_data(device::BrickletIndustrialCounter)
    return device.deviceInternal.get_all_signal_data()
end

export set_counter_active
"""
    $(SIGNATURES)

Activates/deactivates the counter of the given channel.

true = activate, false = deactivate.

By default all channels are activated.
"""
function set_counter_active(device::BrickletIndustrialCounter, channel, active)
    device.deviceInternal.set_counter_active(channel, active)
end

export set_all_counter_active
"""
    $(SIGNATURES)

Activates/deactivates the counter of all four channels.

true = activate, false = deactivate.

By default all channels are activated.
"""
function set_all_counter_active(device::BrickletIndustrialCounter, active)
    device.deviceInternal.set_all_counter_active(active)
end

export get_counter_active
"""
    $(SIGNATURES)

Returns the activation state of the given channel.

true = activated, false = deactivated.
"""
function get_counter_active(device::BrickletIndustrialCounter, channel)
    return device.deviceInternal.get_counter_active(channel)
end

export get_all_counter_active
"""
    $(SIGNATURES)

Returns the activation state of all four channels.

true = activated, false = deactivated.
"""
function get_all_counter_active(device::BrickletIndustrialCounter)
    return device.deviceInternal.get_all_counter_active()
end

export set_counter_configuration
"""
    $(SIGNATURES)

Sets the counter configuration for the given channel.

* Count Edge: Counter can count on rising, falling or both edges.
* Count Direction: Counter can count up or down. You can also use
  another channel as direction input, see
  `here <https://www.tinkerforge.com/en/doc/Hardware/Bricklets/Industrial_Counter.html#external-count-direction>`__
  for details.
* Duty Cycle Prescaler: Sets a divider for the internal clock. See
  `here <https://www.tinkerforge.com/en/doc/Hardware/Bricklets/Industrial_Counter.html#duty-cycle-prescaler-and-frequency-integration-time>`__
  for details.
* Frequency Integration Time: Sets the integration time for the
  frequency measurement. See
  `here <https://www.tinkerforge.com/en/doc/Hardware/Bricklets/Industrial_Counter.html#duty-cycle-prescaler-and-frequency-integration-time>`__
  for details.
"""
function set_counter_configuration(device::BrickletIndustrialCounter, channel, count_edge, count_direction, duty_cycle_prescaler, frequency_integration_time)
    device.deviceInternal.set_counter_configuration(channel, count_edge, count_direction, duty_cycle_prescaler, frequency_integration_time)
end

export get_counter_configuration
"""
    $(SIGNATURES)

Returns the counter configuration as set by :func:`Set Counter Configuration`.
"""
function get_counter_configuration(device::BrickletIndustrialCounter, channel)
    return device.deviceInternal.get_counter_configuration(channel)
end

export set_all_counter_callback_configuration
"""
    $(SIGNATURES)

The period is the period with which the :cb:`All Counter`
callback is triggered periodically. A value of 0 turns the callback off.

If the `value has to change`-parameter is set to true, the callback is only
triggered after the value has changed. If the value didn't change within the
period, the callback is triggered immediately on change.

If it is set to false, the callback is continuously triggered with the period,
independent of the value.
"""
function set_all_counter_callback_configuration(device::BrickletIndustrialCounter, period, value_has_to_change)
    device.deviceInternal.set_all_counter_callback_configuration(period, value_has_to_change)
end

export get_all_counter_callback_configuration
"""
    $(SIGNATURES)

Returns the callback configuration as set by
:func:`Set All Counter Callback Configuration`.
"""
function get_all_counter_callback_configuration(device::BrickletIndustrialCounter)
    return device.deviceInternal.get_all_counter_callback_configuration()
end

export set_all_signal_data_callback_configuration
"""
    $(SIGNATURES)

The period is the period with which the :cb:`All Signal Data`
callback is triggered periodically. A value of 0 turns the callback off.

If the `value has to change`-parameter is set to true, the callback is only
triggered after the value has changed. If the value didn't change within the
period, the callback is triggered immediately on change.

If it is set to false, the callback is continuously triggered with the period,
independent of the value.
"""
function set_all_signal_data_callback_configuration(device::BrickletIndustrialCounter, period, value_has_to_change)
    device.deviceInternal.set_all_signal_data_callback_configuration(period, value_has_to_change)
end

export get_all_signal_data_callback_configuration
"""
    $(SIGNATURES)

Returns the callback configuration as set by
:func:`Set All Signal Data Callback Configuration`.
"""
function get_all_signal_data_callback_configuration(device::BrickletIndustrialCounter)
    return device.deviceInternal.get_all_signal_data_callback_configuration()
end

export set_channel_led_config
"""
    $(SIGNATURES)

Each channel has a corresponding LED. You can turn the LED off, on or show a
heartbeat. You can also set the LED to "Channel Status". In this mode the
LED is on if the channel is high and off otherwise.
"""
function set_channel_led_config(device::BrickletIndustrialCounter, channel, config)
    device.deviceInternal.set_channel_led_config(channel, config)
end

export get_channel_led_config
"""
    $(SIGNATURES)

Returns the channel LED configuration as set by :func:`Set Channel LED Config`
"""
function get_channel_led_config(device::BrickletIndustrialCounter, channel)
    return device.deviceInternal.get_channel_led_config(channel)
end

export get_spitfp_error_count
"""
    $(SIGNATURES)

Returns the error count for the communication between Brick and Bricklet.

The errors are divided into

* ACK checksum errors,
* message checksum errors,
* framing errors and
* overflow errors.

The errors counts are for errors that occur on the Bricklet side. All
Bricks have a similar function that returns the errors on the Brick side.
"""
function get_spitfp_error_count(device::BrickletIndustrialCounter)
    return device.deviceInternal.get_spitfp_error_count()
end

export set_bootloader_mode
"""
    $(SIGNATURES)

Sets the bootloader mode and returns the status after the requested
mode change was instigated.

You can change from bootloader mode to firmware mode and vice versa. A change
from bootloader mode to firmware mode will only take place if the entry function,
device identifier and CRC are present and correct.

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function set_bootloader_mode(device::BrickletIndustrialCounter, mode)
    return device.deviceInternal.set_bootloader_mode(mode)
end

export get_bootloader_mode
"""
    $(SIGNATURES)

Returns the current bootloader mode, see :func:`Set Bootloader Mode`.
"""
function get_bootloader_mode(device::BrickletIndustrialCounter)
    return device.deviceInternal.get_bootloader_mode()
end

export set_write_firmware_pointer
"""
    $(SIGNATURES)

Sets the firmware pointer for :func:`Write Firmware`. The pointer has
to be increased by chunks of size 64. The data is written to flash
every 4 chunks (which equals to one page of size 256).

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function set_write_firmware_pointer(device::BrickletIndustrialCounter, pointer)
    device.deviceInternal.set_write_firmware_pointer(pointer)
end

export write_firmware
"""
    $(SIGNATURES)

Writes 64 Bytes of firmware at the position as written by
:func:`Set Write Firmware Pointer` before. The firmware is written
to flash every 4 chunks.

You can only write firmware in bootloader mode.

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function write_firmware(device::BrickletIndustrialCounter, data)
    return device.deviceInternal.write_firmware(data)
end

export set_status_led_config
"""
    $(SIGNATURES)

Sets the status LED configuration. By default the LED shows
communication traffic between Brick and Bricklet, it flickers once
for every 10 received data packets.

You can also turn the LED permanently on/off or show a heartbeat.

If the Bricklet is in bootloader mode, the LED is will show heartbeat by default.
"""
function set_status_led_config(device::BrickletIndustrialCounter, config)
    device.deviceInternal.set_status_led_config(config)
end

export get_status_led_config
"""
    $(SIGNATURES)

Returns the configuration as set by :func:`Set Status LED Config`
"""
function get_status_led_config(device::BrickletIndustrialCounter)
    return device.deviceInternal.get_status_led_config()
end

export get_chip_temperature
"""
    $(SIGNATURES)

Returns the temperature as measured inside the microcontroller. The
value returned is not the ambient temperature!

The temperature is only proportional to the real temperature and it has bad
accuracy. Practically it is only useful as an indicator for
temperature changes.
"""
function get_chip_temperature(device::BrickletIndustrialCounter)
    return device.deviceInternal.get_chip_temperature()
end

export reset
"""
    $(SIGNATURES)

Calling this function will reset the Bricklet. All configurations
will be lost.

After a reset you have to create new device objects,
calling functions on the existing ones will result in
undefined behavior!
"""
function reset(device::BrickletIndustrialCounter)
    device.deviceInternal.reset()
end

export write_uid
"""
    $(SIGNATURES)

Writes a new UID into flash. If you want to set a new UID
you have to decode the Base58 encoded UID string into an
integer first.

We recommend that you use Brick Viewer to change the UID.
"""
function write_uid(device::BrickletIndustrialCounter, uid)
    device.deviceInternal.write_uid(uid)
end

export read_uid
"""
    $(SIGNATURES)

Returns the current UID as an integer. Encode as
Base58 to get the usual string version.
"""
function read_uid(device::BrickletIndustrialCounter)
    return device.deviceInternal.read_uid()
end

export get_identity
"""
    $(SIGNATURES)

Returns the UID, the UID where the Bricklet is connected to,
the position, the hardware and firmware version as well as the
device identifier.

The position can be 'a', 'b', 'c', 'd', 'e', 'f', 'g' or 'h' (Bricklet Port).
A Bricklet connected to an :ref:`Isolator Bricklet <isolator_bricklet>` is always at
position 'z'.

The device identifier numbers can be found :ref:`here <device_identifier>`.
|device_identifier_constant|
"""
function get_identity(device::BrickletIndustrialCounter)
    return device.deviceInternal.get_identity()
end

export register_callback
"""
Registers the given *function* with the given *callback_id*.
"""
function register_callback(device::BrickletIndustrialCounter, callback_id, function_)
    if isnothing(function_)
        device.registered_callbacks.pop(callback_id, None)
    else
        device.registered_callbacks[callback_id] = function_
    end
end


text:
using HypothesisTests
using Test
using StatsBase
using Statistics
using DelimitedFiles

@testset "Equality of Variances" begin

    # https://en.wikipedia.org/wiki/One-way_analysis_of_variance#Example
    groups = [
        [6, 8, 4, 5, 3, 4],
        [8, 12, 9, 11, 6, 8],
        [13, 9, 11, 8, 7, 12]
    ]
    t = OneWayANOVATest(groups...)
    @test nobs(t) == fill(6, 3)
    @test dof(t) == (2,15)
    @test pvalue(t) ≈ 0.002 atol=1e-3
    @test occursin("reject h_0", sprint(show, t))

    show(IOContext(IOBuffer(), :table => true), t)
    show(IOBuffer(), t)

    # http://www.real-statistics.com/one-way-analysis-of-variance-anova/confidence-interval-anova/
    groups = [
        [51, 87, 50, 48, 79, 61, 53],
        [82, 91, 92, 80, 52, 79, 73, 74],
        [79, 84, 74, 98, 63, 83, 85, 58],
        [85, 80, 65, 71, 67, 51],
    ]
    t = OneWayANOVATest(groups...)
    @test nobs(t) == [7, 8, 8, 6]
    @test dof(t) == (3,25)
    @test pvalue(t) ≈ 0.072 atol=1e-3
    @test occursin("reject h_0", sprint(show, t))

    # http://www.real-statistics.com/one-way-analysis-of-variance-anova/homogeneity-variances/levenes-test/
    groups = [
        [51, 87, 50, 48, 79, 61, 53, 54],
        [82, 91, 92, 80, 52, 79, 73, 74],
        [79, 84, 74, 98, 63, 83, 85, 58],
        [85, 80, 65, 71, 67, 51, 63, 93],
    ]
    # with means
    l = LeveneTest(groups...; statistic=mean)
    @test nobs(l) == fill(8, 4)
    @test dof(l) == (3,28)
    @test pvalue(l) ≈ 0.90357 atol=1e-4
    # with medians
    l = LeveneTest(groups...; statistic=median)
    @test pvalue(l) ≈ 0.97971 atol=1e-4
    # with 10% trimmed means
    l = LeveneTest(groups...; statistic=v->mean(trim(v, prop=0.1)))
    @test pvalue(l) ≈ 0.90357 atol=1e-4

    # Fligner-Killeen Test
    t = FlignerKilleenTest(groups...)
    @test nobs(t) == fill(8, 4)
    @test dof(t) == 3
    @test pvalue(t) ≈ 0.9878 atol=1e-4
    @test HypothesisTests.teststatistic(t) ≈ 0.1311 atol=1e-5
    @test occursin("reject h_0", sprint(show, t))

    # https://www.itl.nist.gov/div898/handbook/eda/section3/eda35a.htm

    # Columns are gear diameter and batch number
    gear = readdlm(joinpath(@__DIR__, "data", "gear.txt"))
    samples = reshape(gear[:,1],:,10)
    groups2 = tuple((samples[:,i] for i in 1:size(samples,1))...)

    # Brown-Forsythe Test
    l = BrownForsytheTest(groups2...)
    @test nobs(l) == fill(10, 10)
    @test dof(l) == (9,90)
    @test HypothesisTests.teststatistic(l) ≈ 1.705910 atol=1e-5
    @test pvalue(l) ≈ 0.0991 atol=1e-4
    @test occursin("reject h_0", sprint(show, l))
end


text:
@testset "buffers" begin
    include("circular_array_buffer.jl")
    include("circular_turn_buffer.jl")
    include("episode_turn_buffer.jl")
end


text:
using Test, DynamicEnergyBudgets, Unitful
using DynamicEnergyBudgets: STATELEN, Records, define_organs, split_state, sum_flux!, offset_apply!

# @testset "build records" begin
    # time=0u"hr":1u"hr":10u"hr"
    # a = build_record(build_flux(1.0u"mol/hr", [:one,:two], [:A,:B]), time)
    # @test axisnames(a) == (:time,)
    # @test axisvalues(a) == (0u"hr":1u"hr":10u"hr",)
    # @test axisnames(a[1]) == (:state, :transformations)
    # @test axisvalues(a[1]) == (([:one, :two]), ([:A, :B]))
    # @test a[1][:one,:A] == 0.0u"mol/hr"
    # @test a[2][:one,:A] == 0.0u"mol/hr"
    # a[1][:one,:A] = 10.0u"mol/hr"
    # @test a[0.0u"hr"][:one,:A] == 10.0u"mol/hr"
    # @test a[1.0u"hr"][:one,:A] == 0.0u"mol/hr"
    # o = Organ()
    # Records(o, time)
# end

@testset "sum_flux" begin
    du = fill(0.0u"mol/hr", 12)
    o = Organism();
    or1, or2 = define_organs(o, 1);
    or1.J .= [1.0,2.0,3.0,4.0,5.0,6.0]oneunit(or1.J[1,1])
    or2.J .= [1.0,2.0,3.0,4.0,5.0,6.0]oneunit(or2.J[1,1]) / 2
    sum_flux!(du, (or1, or2))
    @test du == [7.0,14.0,21.0,28.0,35.0,42.0,3.5,7.0,10.5,14.0,17.5,21.0] .* u"mol/hr"
end

@testset "split_state" begin
    o = Organism();
    os = define_organs(o, 1);
    u = [1.0:12.0...]u"mol"
    us = split_state(os, u)
    @test us[1][1] == 1.0u"mol"
    @test us[2][1] == 7.0u"mol"
    @test sum(us[1]) == sum(1:6) .* u"mol"
    @test sum(us[2]) == sum(7:12) .* u"mol"
end


text:
module A
struct Null end
end

StructTypes.StructType(::Type{A.Null}) = StructTypes.NullType()
StructTypes.construct(::Type{A.Null}, ::Nothing) = nothing
Base.convert(::Type{Union{A.Null, Snowflake}}, s::AbstractString) = Snowflake(s)
Base.convert(::Type{Union{A.Null, Snowflake}}, n::Integer) = Snowflake(n)

function json(x)
    data = Dict{Symbol, Any}()
    for name in propertynames(x)
        val = getproperty(x, name)
        if !ismissing(val)
            data[name] = val
        end
    end
    return JSON3.write(data)
end

macro discord_object(typedef)
    # Make the type mutable.
    typedef.args[1] = true

    # Make all the types nullable/optional with a default value of `missing`.
    fields = typedef.args[3].args
    for (i, field) in enumerate(fields)
        if field isa Expr
            type = field.args[2]
            field.args[2] = :(Union{$type, A.Null, Nothing, Missing})
            fields[i] = :($field = missing)
        end
    end
    typedef_withkw = esc(:(@with_kw $typedef))

    # Enable JSON IO for the type.
    name = esc(typedef.args[2])
    json_methods = quote
        StructTypes.StructType(::Type{$name}) = StructTypes.Mutable()
        JSON3.write(x::$name) = json(x)
    end

    return quote
        export $name
        $typedef_withkw
        $json_methods
    end
end

macro discord_enum(name, block, kwargs...)
    # TODO: Why do I need to eval here? I want to just insert a macrocall expression.
    @eval @se $name $block

    exp = :(export $(esc(name)))
    name = Expr(:., name, QuoteNode(Symbol(name, :Enum)))

    extras = Expr(:block)
    for ex in kwargs
        k, v = ex.args
        if v === true
            if k === :or
                or = :(Base.:(|)(a::$(esc(name)), b::$(esc(name))) = Int(a) | Int(b))
                push!(extras.args, or)
            end
        end
    end

    return quote
        StructTypes.StructType(::Type{$(esc(name))}) = StructTypes.NumberType()
        StructTypes.numbertype(::Type{$(esc(name))}) = Int
        $exp
        $extras
    end
end

include(joinpath("objects", "enums.jl"))
include(joinpath("objects", "user.jl"))
include(joinpath("objects", "integration_account.jl"))
include(joinpath("objects", "integration.jl"))
include(joinpath("objects", "connection.jl"))
include(joinpath("objects", "overwrite.jl"))
include(joinpath("objects", "channel.jl"))
include(joinpath("objects", "channel_mention.jl"))
include(joinpath("objects", "guild_member.jl"))
include(joinpath("objects", "attachment.jl"))
include(joinpath("objects", "embed_footer.jl"))
include(joinpath("objects", "embed_image.jl"))
include(joinpath("objects", "embed_thumbnail.jl"))
include(joinpath("objects", "embed_video.jl"))
include(joinpath("objects", "embed_provider.jl"))
include(joinpath("objects", "embed_author.jl"))
include(joinpath("objects", "embed_field.jl"))
include(joinpath("objects", "embed.jl"))
include(joinpath("objects", "emoji.jl"))
include(joinpath("objects", "reaction.jl"))
include(joinpath("objects", "message_activity.jl"))
include(joinpath("objects", "message_application.jl"))
include(joinpath("objects", "message_reference.jl"))
include(joinpath("objects", "message.jl"))
include(joinpath("objects", "role.jl"))
include(joinpath("objects", "voice_state.jl"))
include(joinpath("objects", "activity_timestamps.jl"))
include(joinpath("objects", "activity_emoji.jl"))
include(joinpath("objects", "activity_party.jl"))
include(joinpath("objects", "activity_assets.jl"))
include(joinpath("objects", "activity_secrets.jl"))
include(joinpath("objects", "activity.jl"))
include(joinpath("objects", "client_status.jl"))
include(joinpath("objects", "presence.jl"))
include(joinpath("objects", "guild.jl"))
include(joinpath("objects", "guild_widget.jl"))
include(joinpath("objects", "webhook.jl"))
include(joinpath("objects", "audit_log_change.jl"))
include(joinpath("objects", "audit_log_info.jl"))
include(joinpath("objects", "audit_log_entry.jl"))
include(joinpath("objects", "audit_log.jl"))
include(joinpath("objects", "allowed_mentions.jl"))
include(joinpath("objects", "ban.jl"))
include(joinpath("objects", "invite.jl"))
include(joinpath("objects", "voice_region.jl"))
include(joinpath("objects", "prune_count.jl"))
include(joinpath("objects", "user_nick_change.jl"))


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

using Markdown
using InteractiveUtils

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

# ╔═╡ 79b40eaa-f6b6-11ea-05a0-d3a166f4fd39
using OrdinaryDiffEq, LightGraphs, GraphPlot, Random, LinearAlgebra, Plots, PlutoUI

# ╔═╡ 2109877c-f6b5-11ea-1941-a7f99c8bab86
md"# Graph Matching using Gradient Flow in ``\cal{O}_n``"

# ╔═╡ ed92164a-8b01-11eb-031d-7997c875255a
md"""
This demonstrates the findings of [Zavlanos and Pappas, 2008](https://doi.org/10.1016/j.automatica.2008.04.009) on using *Analog Computing* in the form of gradient flows to solve graph matching problems.
"""

# ╔═╡ f0fb46aa-f6d8-11ea-2467-1b2fa447a503
md"## Original and permuted graph"

# ╔═╡ 999b3430-8b02-11eb-1076-354db3137ad1
md"""
We start by constructing a simple graph. Just chose one of the built-in examples.
"""

# ╔═╡ 1c300130-fa76-11ea-3a12-47d1dffd179b
md"""**Select example graph:** $(@bind example_graph Select(
	[ "chain5" => "chain of 5"
	, "chain20" => "chain of 20"
	, "tutte" => "tutte"
	, "housex" => "house"
	]))
"""

# ╔═╡ 12ec0c50-8b04-11eb-0c2d-714d2edc1f4b
md"""
And here is a visualisation of that graph.
(The layouting is non-deterministic, so try rerunning the cell if it looks strange.)
"""

# ╔═╡ 40af3fc4-8b04-11eb-2b1c-9b6658331a25
md"""
We consider `G1` to be the "original" graph.
Let's construct another graph, `G2`, by permuting all vertex labels in `G1`:
"""

# ╔═╡ ca129dcc-8b04-11eb-2f89-7f6855899278
md"""
Our **goal** will be to find a permuation of the vertices in `G2` such that we obtain `G1` again.
"""

# ╔═╡ ba799a4e-fa79-11ea-07fb-4fa60e4829d5
md"### Adjacency matrices"

# ╔═╡ 9eb927fe-8b04-11eb-0699-47c6e029292c
md"""
We will use some linear algebra to solve the graph matching.
To do this, we need canonical matrix representation of a graph: its *adjacency matrix*.

Recall that for an undirected graph ``G = (V, E)``, its adjacency matrix ``A \in \{0, 1\}^{|V| \times |V|}`` has entries
```math
a_{ij} = \begin{cases}
	1, & \text{if } \{i, j\} \in E \\
	0, & \text{else.}
\end{cases}
```
"""

# ╔═╡ be035310-8b08-11eb-26af-2d8fe12b1a29
md"""
## Our Strategy
We briefly restate the strategy that Zavlanos and Pappas suggest:

We want to permute the vertices in `G2`, which is equivalent to permuting the rows and columns in `A2`.
More specifically, we want to find a permuatation matrix ``P`` such that ``A_1`` and ``P^\top A_2 P`` (or equivalently ``P A_1`` and ``A_2 P``) are as similar as possible.
To quantify the approximation error, we pick the squared *Frobenius*-norm as our notion of similarity, which just means that we sum the elementwise squared differences of both matrices:
"""

# ╔═╡ 3a7dd4d2-8b0f-11eb-0658-f3d3779b33f3
md"### Minimise approximation error"

# ╔═╡ c748b142-8b0b-11eb-0199-c5924c5d19cf
md"""
The authors then explain that the permutation matrices are exactly the orthogonal matrices with non-negative entries.
We exploit that by first relaxing our problem to ``P`` being any orthogonal matrix and performing a **gradient flow**.

This is the ODE problem and its solution that describes the gradient flow:
"""

# ╔═╡ 262ef806-f6c4-11ea-1b96-8776b89f8071
function min_approx_err!(dP, P, params, t)
	A1, A2 = params
	dP .= P * (P' * A2 * P * A1 - A1 * P' * A2 * P)
end

# ╔═╡ 9d320016-8b0d-11eb-0d83-45212b00c532
md"Define for how long we want to follow the gradient flow:"

# ╔═╡ 6092784c-f741-11ea-1a0b-3b40c5a88a67
t1 = 20.0

# ╔═╡ 1a9da0ce-f6da-11ea-24a1-398e86eea280
md"""
### Penalising negative entries in ``P``

We now need to make the orthogonal matrix into a permutation matrix, which means we have to get rid of the negative entries.

This function assesses how negative the entries in ``P`` are:
"""

# ╔═╡ ca89c966-f6cf-11ea-05d8-717423c34f3c
negativity(P) = 2/3 * tr(P' * (P .- P .* P))

# ╔═╡ f161ae44-8b0f-11eb-07fc-5364c9e5aa55
md"""
And here we define the corresponding gradient flow:
"""

# ╔═╡ 222dbe04-f6d0-11ea-2997-b5d642fddf4c
function penalise_negative!(dP, P, params, t)
	PP = P .* P
	dP .= P * (P' * PP .- PP' * P)
end

# ╔═╡ 5d6b088e-8b0f-11eb-26f6-dfa9cacb669e
md"""
It turns out that it makes sense to actually combine the gradient flows for the actual objective and for penalising negative entries with some weighting factor `k`.
"""

# ╔═╡ cb658266-fa79-11ea-0db1-7567dea3d069
k = 0.5

# ╔═╡ 490ce0c0-f72f-11ea-1466-3b1aeb47fc11
function combined_gradient!(dP, P, params, t)
	A1, A2, k = params
	dP1 = similar(dP)
	dP2 = similar(dP)
	min_approx_err!(dP1, P, (A1, A2), t)
	penalise_negative!(dP2, P, nothing, t)
	dP .= (1-k) * dP1 .+ k * dP2
end

# ╔═╡ 6767b2d2-8b10-11eb-116c-1972bef114fe
md"""
And this is the overall ODE problem describing this second gradient flow, which starts at the last ``P`` of the first system and evolves until `t2`.
"""

# ╔═╡ cabbd10e-fa79-11ea-30f2-bd1c70613a3c
t2 = 60.0

# ╔═╡ e63643ac-f739-11ea-32ac-d124ec648ff8
md"""
We can summarise the strategy:

**Minimise approximation error until time $t1, then additionally make entries in ``P`` non-negative until time $t2 with $(round(Int, 100k)) % importance.**
"""

# ╔═╡ ba3571ac-8b10-11eb-3654-f3861db16d26
md"""
## Results

This combines both parts of the gradient flow into one time-dependent function:
"""

# ╔═╡ 05f036a2-fa7c-11ea-1902-7592154ca34d
md"### Plotting"

# ╔═╡ a730b14e-f74e-11ea-3ec9-ed8b3a479bdc
trange = 0 : 1.0 : t2;

# ╔═╡ 02cf1af8-8b11-11eb-00a1-2b2e298b4a85
md"""
This is the final solution for ``P`` that the algorithm came up with:
"""

# ╔═╡ 17f5d944-8b11-11eb-0f62-c5e233214748
md"""
Sanity check, that the final solution is still orthogonal:
"""

# ╔═╡ 369c1c8a-f6da-11ea-2069-3fde5a57b953
md"""
## Reconstructed graph

Finally, we can perfom the permutation and find the reconstructed adjacency matrix:
"""

# ╔═╡ 59b29f46-8b11-11eb-2b3c-277d36a0debb
md"""
And this is the corresponding graph.
Does it look at least somewhat like `G1`?
I hope so!
Try playing with the `k` parameter to see how that affects the solution quality.
"""

# ╔═╡ 64e538b2-f6bf-11ea-375e-a5f3e77d6cd0
function permute_vertices(g::AbstractGraph)
	perm_g = SimpleGraph(nv(g))
	perm_v = shuffle(vertices(g))

	for e in edges(g)
		perm_s = perm_v[src(e)]
		perm_d = perm_v[dst(e)]
		add_edge!(perm_g, perm_s, perm_d)
	end

	perm_g
end

# ╔═╡ 12c0e9ee-fa77-11ea-3e61-ed34d1ea848e
function construct_example(id)
	if id == "tutte"
		smallgraph(:tutte)
	elseif id == "housex"
		smallgraph(:housex)
	elseif id == "chain5"
		path_graph(5)
	elseif id == "chain20"
		path_graph(20)
	else
		path_graph(1)
	end
end

# ╔═╡ c548c0b2-f6bc-11ea-1ba0-f3de59ef9987
G1 = construct_example(example_graph)

# ╔═╡ a9ac43ac-f6c2-11ea-2114-338b468c054f
gplot(G1, nodelabel=collect(vertices(G1)))

# ╔═╡ 14eeb4fe-f6c0-11ea-1211-972d5d4f6a4b
G2 = permute_vertices(G1)

# ╔═╡ 78bb17c8-f6c2-11ea-04df-69c5d723ec3c
gplot(G2, nodelabel=collect(vertices(G2)))

# ╔═╡ d8cbb096-f6c2-11ea-1f32-ed95b80b4d8f
A2 = adjacency_matrix(G2);

# ╔═╡ cc0052ea-f6c2-11ea-2710-e175f7ec3f8a
A1 = adjacency_matrix(G1);

# ╔═╡ 637607f8-f6cc-11ea-2a32-a15f4d93c7e7
approx_err(P) = .5 * sum((P*A1 - A2*P) .^ 2)

# ╔═╡ afba6f40-f6c6-11ea-20a3-e16e48da9005
function rand_SO(s::Tuple{Int,Int})
	Q, _ = rand(s...) |> qr
	Q = Matrix(Q)
	if det(Q) < 0
		Q[1,:] .*= -1
	end
	Q
end

# ╔═╡ f4fd59b4-f6c6-11ea-165b-eb798b8fbff5
min_approx_err_solution = let
	P0 = rand_SO(size(A1))
	t0 = 0.0
	min_approx_err_problem = ODEProblem(min_approx_err!, P0, (t0, t1), (A1, A2))

	solve(min_approx_err_problem, Tsit5())
end;

# ╔═╡ 19739754-fa7a-11ea-2b0e-037d4266a5ae
P_min_approx_err = min_approx_err_solution(t1);

# ╔═╡ 7519397a-f6d0-11ea-19ad-f5aaa7f8f228
non_negative_solution = let
	non_negative_problem = ODEProblem(combined_gradient!, P_min_approx_err, (t1, t2), (A1, A2, k))
	solve(non_negative_problem, Tsit5())
end;

# ╔═╡ f6580a7e-f6d1-11ea-2200-f9c80bc69b31
P(t) = if t < t1
		min_approx_err_solution(t)
	else
		non_negative_solution(t)
	end

# ╔═╡ 3358b806-f6d2-11ea-203e-d53ae160ffd6
anim = map(trange) do t
	l = heatmap(P(t) * A1 - A2 * P(t) .|> abs,
		clim=(0,2),
		title="t = $(round(Int, t))",
		xlabel="elem. wise deviation")
	r = heatmap(P(t),
		clim=(-1,1),
		xlabel="permutation matrix")
	plot(l, r, layout=(1,2), size=(600,300), aspect_ratio=:equal)
end;

# ╔═╡ ce19a088-fa7b-11ea-30b3-7b80bfbbdf1f
@bind t_anim Slider(eachindex(anim))

# ╔═╡ ec1ba20c-fa7b-11ea-1238-41e99f6ded6b
anim[t_anim]

# ╔═╡ 507ccbd4-f6d2-11ea-3a8d-1311265cdc54
begin
	plot(trange,
		[[approx_err(P(t)) for t in trange],
		[negativity(P(t)) for t in trange]],
		label=["approximation error"  "negativity of entries"],
		xlabel="time")
	vline!([t1]; linestyle=:dash, label="switching systems")
	vline!([trange[t_anim]]; linestyle=:solid, label="currently displayed")
end

# ╔═╡ 422d4c82-f6d4-11ea-3e87-d7f548a73eb7
final_P = P(t2);

# ╔═╡ e527660c-fc60-11ea-309e-c7372eea3e23
round.(final_P, digits=3)

# ╔═╡ 41f9aea8-f6d5-11ea-0330-af1c77272674
round.(final_P' * final_P, digits=3)

# ╔═╡ a184f726-f6d8-11ea-2c93-8fb47ac719a3
A3 = round.(Int, final_P' * A2 * final_P)

# ╔═╡ 2735a812-f6d8-11ea-16c4-dbbc3b090adb
G3 = SimpleGraph(A3)

# ╔═╡ 7649f264-f6d8-11ea-1712-5dac97376cdc
gplot(G3, nodelabel=collect(vertices(G3)))

# ╔═╡ 97c742fc-f6d5-11ea-1b58-6be1318afa25
begin
	plot(trange,
		[minimum(P(t)) for t in trange],
		ylim=(-1,0.1),
		label="min. entry in P",
		xlabel="time")
	vline!([t1]; linestyle=:dash, label="switching systems", legend=:bottomright)
end

# ╔═╡ Cell order:
# ╟─2109877c-f6b5-11ea-1941-a7f99c8bab86
# ╟─ed92164a-8b01-11eb-031d-7997c875255a
# ╟─f0fb46aa-f6d8-11ea-2467-1b2fa447a503
# ╟─999b3430-8b02-11eb-1076-354db3137ad1
# ╟─1c300130-fa76-11ea-3a12-47d1dffd179b
# ╠═c548c0b2-f6bc-11ea-1ba0-f3de59ef9987
# ╟─12ec0c50-8b04-11eb-0c2d-714d2edc1f4b
# ╠═a9ac43ac-f6c2-11ea-2114-338b468c054f
# ╟─40af3fc4-8b04-11eb-2b1c-9b6658331a25
# ╠═14eeb4fe-f6c0-11ea-1211-972d5d4f6a4b
# ╠═78bb17c8-f6c2-11ea-04df-69c5d723ec3c
# ╟─ca129dcc-8b04-11eb-2f89-7f6855899278
# ╟─ba799a4e-fa79-11ea-07fb-4fa60e4829d5
# ╟─9eb927fe-8b04-11eb-0699-47c6e029292c
# ╠═cc0052ea-f6c2-11ea-2710-e175f7ec3f8a
# ╠═d8cbb096-f6c2-11ea-1f32-ed95b80b4d8f
# ╟─be035310-8b08-11eb-26af-2d8fe12b1a29
# ╠═637607f8-f6cc-11ea-2a32-a15f4d93c7e7
# ╟─3a7dd4d2-8b0f-11eb-0658-f3d3779b33f3
# ╟─c748b142-8b0b-11eb-0199-c5924c5d19cf
# ╠═262ef806-f6c4-11ea-1b96-8776b89f8071
# ╟─9d320016-8b0d-11eb-0d83-45212b00c532
# ╠═6092784c-f741-11ea-1a0b-3b40c5a88a67
# ╠═f4fd59b4-f6c6-11ea-165b-eb798b8fbff5
# ╠═19739754-fa7a-11ea-2b0e-037d4266a5ae
# ╟─1a9da0ce-f6da-11ea-24a1-398e86eea280
# ╠═ca89c966-f6cf-11ea-05d8-717423c34f3c
# ╟─f161ae44-8b0f-11eb-07fc-5364c9e5aa55
# ╠═222dbe04-f6d0-11ea-2997-b5d642fddf4c
# ╟─5d6b088e-8b0f-11eb-26f6-dfa9cacb669e
# ╠═cb658266-fa79-11ea-0db1-7567dea3d069
# ╠═490ce0c0-f72f-11ea-1466-3b1aeb47fc11
# ╟─6767b2d2-8b10-11eb-116c-1972bef114fe
# ╠═cabbd10e-fa79-11ea-30f2-bd1c70613a3c
# ╠═7519397a-f6d0-11ea-19ad-f5aaa7f8f228
# ╟─e63643ac-f739-11ea-32ac-d124ec648ff8
# ╟─ba3571ac-8b10-11eb-3654-f3861db16d26
# ╠═f6580a7e-f6d1-11ea-2200-f9c80bc69b31
# ╟─05f036a2-fa7c-11ea-1902-7592154ca34d
# ╠═a730b14e-f74e-11ea-3ec9-ed8b3a479bdc
# ╟─3358b806-f6d2-11ea-203e-d53ae160ffd6
# ╟─ce19a088-fa7b-11ea-30b3-7b80bfbbdf1f
# ╟─ec1ba20c-fa7b-11ea-1238-41e99f6ded6b
# ╟─507ccbd4-f6d2-11ea-3a8d-1311265cdc54
# ╟─02cf1af8-8b11-11eb-00a1-2b2e298b4a85
# ╠═422d4c82-f6d4-11ea-3e87-d7f548a73eb7
# ╠═e527660c-fc60-11ea-309e-c7372eea3e23
# ╟─17f5d944-8b11-11eb-0f62-c5e233214748
# ╠═41f9aea8-f6d5-11ea-0330-af1c77272674
# ╟─97c742fc-f6d5-11ea-1b58-6be1318afa25
# ╟─369c1c8a-f6da-11ea-2069-3fde5a57b953
# ╠═a184f726-f6d8-11ea-2c93-8fb47ac719a3
# ╟─59b29f46-8b11-11eb-2b3c-277d36a0debb
# ╠═2735a812-f6d8-11ea-16c4-dbbc3b090adb
# ╠═7649f264-f6d8-11ea-1712-5dac97376cdc
# ╟─64e538b2-f6bf-11ea-375e-a5f3e77d6cd0
# ╟─12c0e9ee-fa77-11ea-3e61-ed34d1ea848e
# ╟─afba6f40-f6c6-11ea-20a3-e16e48da9005
# ╠═79b40eaa-f6b6-11ea-05a0-d3a166f4fd39

max_stars_repo_name: andreasKroepelin/graphmatching

max_issues_repo_name: andreasKroepelin/graphmatching

max_forks_repo_name: andreasKroepelin/graphmatching


text:
let
    method = ME_Z_EXPECTED_G_MOVING

    # Feed jobs
    Ch = Channel(1) do ch
        cGLCs = Hd.val("cGLC")
        for (exp, cGLC)  in enumerate(cGLCs)
            put!(ch, (exp, cGLC))
        end
    end

    # @threads
    for thid in 1:nthreads()
        for (exp, cGLC) in Ch

            ## -------------------------------------------------------------------
            # handle cache
            datfile = dat_file(;method, exp)
            check_cache(;method, exp) && continue

            ## -------------------------------------------------------------------
            # SetUp
            model_cache_id = (:MODEL0_CACHE, exp)
            # UJL.delete_cache(model_cache_id) # uncomment to reset
            model =  UJL.load_cache(model_cache_id; verbose = true) do
                BASE_MODELS = ChU.load_data(iJR.BASE_MODELS_FILE; verbose = false);
                model_dict = BASE_MODELS["fva_models"][exp]
                model0 = ChU.MetNet(;model_dict...) |> ChU.uncompressed_model
                ChU.scale_polytope!(model0, 0.5)
                return ChLP.fva_preprocess(model0;
                    verbose = false, check_obj = iJR.BIOMASS_IDER
                )
            end
            objidx = ChU.rxnindex(model, iJR.BIOMASS_IDER)
            M, N = size(model)
            exp_growth = Hd.val("D", exp)
            biom_lb, biom_ub = ChU.bounds(model, iJR.BIOMASS_IDER)
            if biom_ub < exp_growth
                lock(WLOCK) do
                    INDEX[method, :DFILE, exp] = :unfeasible
                    @info("Not feasible (skipping)",
                        biom_ub ,exp_growth,
                        thid
                    ); println()
                end
                continue
            end
            cgD_X = -Hd.cval(:GLC, exp) * Hd.val(:D, exp) / Hd.val(:DCW, exp)
            exglc_L = ChU.lb(model, iJR.EX_GLC_IDER)
            exglc_qta = abs(exglc_L * 0.005)
            expβ = 0.0
            epouts = Dict()
            epout = nothing

            for movround in 1:1000

                empty!(epouts)

                ## -------------------------------------------------------------------
                # GRAD DESCEND
                x0 = expβ
                x1 = 10.0
                maxΔ = max(expβ * 0.05, 5e3)
                gd_th = 1e-3
                target = exp_growth
                beta_vec = zeros(size(model, 2))

                upfrec_time = 10 # secunds
                last_uptime = time()
                gd_it = 1

                ## -------------------------------------------------------------------
                function upfun(beta)

                    beta_vec[objidx] = beta
                    epout = ChEP.maxent_ep(model;
                        beta_vec,
                        alpha = Inf,
                        maxiter = 5000,
                        epsconv = 1e-3,
                        verbose = false,
                        solution = epout
                    )
                    epouts[beta] = epout

                    # ep_growth = ChU.av(epout)[objidx]
                    ep_growth = ChU.av(model, epout, iJR.BIOMASS_IDER)

                    update = gd_it == 1 || abs(last_uptime - time()) > upfrec_time ||
                        epout.status != ChEP.CONVERGED_STATUS

                    update && lock(WLOCK) do
                        diff = abs.(exp_growth - ep_growth)
                        @info(
                            "Grad descent... ",
                            exp, gd_it,
                            epout.status, epout.iter,
                            ep_growth, exp_growth, diff,
                            beta, thid
                        ); println()
                        gd_it += 1
                        last_uptime = time()

                    end
                    return ep_growth
                end

                ## -------------------------------------------------------------------
                # FIND BETA
                expβ = UJL.grad_desc(upfun; x0, x1, th = gd_th, maxΔ,
                    target, maxiters = 2000, verbose = false
                )

                ## -------------------------------------------------------------------
                # MOVE V_UB
                Δstep = 0.5
                exglc_lb = ChU.lb(model, iJR.EX_GLC_IDER)
                exglc_ub = ChU.ub(model, iJR.EX_GLC_IDER)
                vg_avPME = ChU.av(model, epout, iJR.EX_GLC_IDER)

                # lb is the uptake limit
                dist = cgD_X - vg_avPME
                Δexglc_lb = sign(dist) * max(exglc_qta, abs(dist * Δstep))
                exglc_lb = min(exglc_ub,
                    max(exglc_L, exglc_lb + Δexglc_lb)
                )
                ChU.lb!(model, iJR.EX_GLC_IDER, exglc_lb)

                ## -------------------------------------------------------------------
                # INFO AND CONV
                ep_growth = ChU.av(model, epout, iJR.BIOMASS_IDER)
                gd_err = abs(exp_growth - ep_growth) / exp_growth
                conv = cgD_X <= vg_avPME && epout.status == ChEP.CONVERGED_STATUS && gd_err < gd_th

                lock(WLOCK) do
                    @info("Round Done",
                        movround, conv,
                        gd_err, exglc_qta, Δexglc_lb,
                        (vg_avPME, cgD_X),
                        exglc_ub, exglc_lb,  exglc_L,
                        thid
                    ); println()
                end
                conv && break

            end #  for movround in 1:1000

            ## -------------------------------------------------------------------
            lock(WLOCK) do
                # Storing
                dat = Dict()
                dat[:exp_beta] = expβ
                dat[:epouts] = epouts
                dat[:model] = model |> ChU.compressed_model

                # caching
                serialize(datfile, dat)

                ep_growth = ChU.av(epouts[expβ])[objidx]
                diff = abs.(exp_growth - ep_growth)
                @info("Finished ",
                    exp, expβ,
                    length(epouts),
                    ep_growth, exp_growth, diff,
                    thid
                ); println()
            end

        end # for (exp, cGLC) in Ch
    end # for thid in 1:nthreads()
end

## -------------------------------------------------------------------
# Further convergence
let
    method = ME_Z_EXPECTED_G_MOVING
    objider = iJR.BIOMASS_IDER

    iterator = Hd.val("cGLC") |> enumerate |> collect
    @threads for (exp, cGLC) in iterator

        datfile = dat_file(;method, exp)
        datfile == :unfeasible && continue
        dat = deserialize(datfile)
        model, epouts = ChU.uncompressed_model(dat[:model]) , dat[:epouts]

        exp_growth = Hd.val(:D, exp)
        exp_beta = dat[:exp_beta]
        exp_epout = epouts[exp_beta]

        lock(WLOCK) do
            @info("Converging...",
                exp, method,
                exp_beta, exp_epout.status,
                threadid()
            ); println()
        end
        converg_status = get!(dat, :converg_status, :undone)
        converg_status == :done && continue

        model = ChLP.fva_preprocess(model; verbose = false,
            check_obj = iJR.BIOMASS_IDER
        )

        new_epout = nothing
        try
            objidx = ChU.rxnindex(model, objider)
            beta_vec = zeros(size(model, 2));
            beta_vec[objidx] = exp_beta
            new_epout = ChEP.maxent_ep(model;
                beta_vec, alpha = Inf,
                epsconv = 1e-5, verbose = false,
                solution = exp_epout, maxiter = 5000
            )
        catch err; @warn("ERROR", err); println() end

        ep_growth = isnothing(new_epout) ? 0.0 : ChU.av(model, new_epout, objider)
        fail = isnan(ep_growth) || ep_growth == 0.0
        epouts[exp_beta] = fail ? exp_epout : new_epout

        # Saving
        lock(WLOCK) do
            @info("Saving...",
                exp, method,
                exp_beta,
                new_epout.status,
                new_epout.iter,
                threadid()
            ); println()
        end
        dat[:model] = ChU.compressed_model(model)
        dat[:converg_status] = :done
        serialize(datfile, dat)
    end
end


text:
module jlpkg

# --compile=min segfaults Julia 1.0 and the fix can't trivially be backported.
# --color should be removed for the new LTS due to better defaults in Julia >1.5
@static if VERSION < v"1.1"
    const default_julia_flags = ["--color=yes", "--startup-file=no", "-q", "-O0"]
else
    const default_julia_flags = ["--color=yes", "--startup-file=no", "-q", "--compile=min", "-O0"]
end

"""
    jlpkg.install(; kwargs...)

Install the command line interface.

*Keyword arguments:*
 - `command`: name of the executable command, defaults to `jlpkg`.
 - `julia`: path to julia executable, defaults to the path of the
   current running julia.
 - `destdir`: writable directory (available in PATH) for the executable,
   defaults to `~/.julia/bin`.
 - `julia_flags`: vector with command line flags for the julia executable,
   defaults to `["--color=yes", "--startup-file=no", "-q", "--compile=min", "-O0"]`.
 - `force`: boolean used to overwrite any existing commands.
"""
function install(; julia::String=joinpath(Sys.BINDIR, Base.julia_exename()),
                   command::String="jlpkg",
                   destdir::String=joinpath(DEPOT_PATH[1], "bin"),
                   julia_flags::Vector{String}=default_julia_flags,
                   force::Bool=false)
    Sys.iswindows() && (command *= ".cmd")
    destdir = abspath(expanduser(destdir))
    exec = joinpath(destdir, command)
    if ispath(exec) && !force
        error("file `$(exec)` already exists; use `jlpkg.install(force=true)` to overwrite.")
    end
    mkpath(destdir)
    open(exec, "w") do f
        if Sys.iswindows()
            # TODO: Find a way to embed the script in the file
            print(f, """
                @ECHO OFF
                $(julia) $(join(julia_flags, ' ')) $(abspath(@__DIR__, "cli.jl")) %*
                """)
        else # unix
            print(f, """
                #!/usr/bin/env bash
                #=
                exec $(julia) $(join(julia_flags, ' ')) "\${BASH_SOURCE[0]}" "\$@"
                =#
                """)
            open(abspath(@__DIR__, "cli.jl"), "r") do cli
                write(f, cli)
            end
        end
    end
    chmod(exec, 0o0100775) # equivalent to -rwxrwxr-x (chmod +x exec)
    @info "Installed jlpkg to `$(Base.contractuser(exec))`."
    if Sys.which(command) === nothing
        @warn """
              `Sys.which(\"$(command)\")` returns `nothing`, meaning that `command`
              can not be found in PATH. Either make sure that `$(destdir)` is in PATH,
              or manually add a symlink from a directory in PATH to the installed
              program file.
              Path to installed program:
                  $(exec)
              """
    elseif realpath(Sys.which(command)) !== realpath(exec)
        @warn """
              `Sys.which(\"$(command)\")` points to a different program than the one
              just installed. Please check your PATH.
              `Sys.which(\"$(command)\")`:
                  $(Sys.which(command))
              Path to installed program:
                  $(exec)
              """
    end
    return
end

"""
    jlpkg.install_completion(; kwargs...)

Install shell completion for jlpkg. Currently only bash completion is supported.

*Keyword arguments:*
 - `command`: name of the executable command to be completed, defaults to `jlpkg`.
 - `destdir`: writable directory to place the completion file in,
   defaults to `~/.bash_completion.d`.
 - `rcfile`: shell startup file to source the completion file in,
   defaults to `~/.bashrc`. If you want to handle sourcing yourself, use `rcfile=nothing`.
 - `force`: boolean used to overwrite an existing completion file.
"""
function install_completion(; command::String="jlpkg",
                              destdir::String=joinpath(homedir(), ".bash_completion.d"),
                              rcfile::Union{String,Nothing}=joinpath(homedir(), ".bashrc"),
                              shell::String="bash",
                              force::Bool=false)
    if shell != "bash"
        throw(ArgumentError("only bash completion is currently supported"))
    end

    # Install script into destdir
    fname = "jlpkg-completion.bash"
    mkpath(destdir)
    dst = joinpath(destdir, fname)
    dstc = Base.contractuser(dst)
    src = joinpath(@__DIR__, fname)
    if isfile(dst) && !force
        throw(ArgumentError("destination file `$(dstc)` already exist; use force=true to overwrite."))
    end
    src_str = read(src, String)
    src_str = replace(src_str, r"^complete -F _jlpkg jlpkg$"m => "complete -F _jlpkg $(command)")
    open(iod -> write(iod, src_str), dst, "w")
    @info "Installed jlpkg bash completion to `$(dstc)`."

    # Source the file in the rc file
    if rcfile !== nothing
        rc = isfile(rcfile) ? read(rcfile, String) : ""
        rcc = Base.contractuser(rcfile)
        header = """

        # Bash completion for jlpkg
        # See https://github.com/fredrikekre/jlpkg for details."""
        if occursin(header, rc)
            @info "`$(rcc)` already appears to source `$(dstc)`, skipping modification."
        else
            open(rcfile, "a") do io
                print(io, """
                $(header)
                if [[ -f $(dstc) ]]; then
                    . $(dstc)
                fi
                """)
            end
            @info "Modified `$(rcc)` to source `$(dstc)`."
        end
    end
    return dst
end

end # module

max_stars_repo_name: fredrikekre/jlpkg

max_issues_repo_name: fredrikekre/jlpkg

max_forks_repo_name: fredrikekre/jlpkg


text:
import StatsBase: countmap

"extend the countmap in StatsBase to support general iterator"
function countmap(iter)
    res = Dict{eltype(iter),Int}()
    for x in iter
        if haskey(res, x)
            res[x] += 1
        else
            res[x] = 1
        end
    end
    res
end


text:
"""
     rouwenhorst(μ,ρ,σ,n)
Generates a finite state Markov representation of the AR process defined by the
supplied parameters. Returns the states and transition matrix.
"""
function rouwenhorst(μ, ρ, σ, n)
    n == 1 && (return [μ], [1.0])

    mu_eps = 0
    q = (ρ + 1) / 2
    nu = ((n - 1) / (1 - ρ^2))^(1 / 2) * σ
    P = [q 1 - q; 1 - q q]
    for i = 2:n - 1
        P = q * [P zeros(i, 1);zeros(1, i + 1)] + (1 - q) * [zeros(i, 1) P; zeros(1, i + 1)] + (1 - q) * [zeros(1, i + 1); P zeros(i, 1)] + q * [zeros(1, i + 1); zeros(i, 1) P]
        P[2:i, :] = P[2:i, :] / 2
    end
    x = [linspace(mu_eps / (1 - ρ) .- nu, mu_eps / (1 - ρ) .+ nu, n);] .+ μ
    return x, P
end

"""
     tauchen(μ,ρ,σ,n)
Generates a finite state Markov representation of the AR process defined by the
supplied parameters. Returns the states and transition matrix.
"""
function tauchen(μ, ρ, σ, N, m = 3)
    cdf_normal(x) = 0.5 * erfc(-x / sqrt(2))
    if N == 1
        return [μ], [1.0]
    end
    Z     = zeros(N)
    Zprob = zeros(N, N)
    a     = (1 - ρ) * μ

    Z[N]  = m * sqrt(σ^2 / (1 - ρ^2))
    Z[1]  = -Z[N]
    zstep = (Z[N] - Z[1]) / (N - 1)

    for i = 2:(N - 1)
        Z[i] = Z[1] + zstep * (i - 1)
    end

    Z = Z .+ a / (1 - ρ)

    for j = 1:N
        for k = 1:N
            if k == 1
                Zprob[j,k] = cdf_normal((Z[1] - a - ρ * Z[j] + zstep / 2) / σ)
            elseif k == N
                Zprob[j,k] = 1 - cdf_normal((Z[N] - a - ρ * Z[j] - zstep / 2) / σ)
            else
                Zprob[j,k] = cdf_normal((Z[k] - a - ρ * Z[j] + zstep / 2) / σ) - cdf_normal((Z[k] - a - ρ * Z[j] - zstep / 2) / σ)
            end
        end
    end
    return Z, Zprob
end

max_stars_repo_name: ZacLN/HADSGE

max_issues_repo_name: ZacLN/HADSGE

max_forks_repo_name: ZacLN/HADSGE


text:
if  pwd() ∉ LOAD_PATH
  push!(LOAD_PATH,pwd())
end

if "$(pwd())\\data" ∉ LOAD_PATH
  push!(LOAD_PATH,"$(pwd())\\data")
end

using CSV, DataFrames, CT, Random

const OUTPUT_PATH = pwd() * "\\output" #path where we store the data
const FOOTER_NAME = "footer.tex" #these are just headers for the tex tables
const HEADER_NAME = "header.tex"

function mwe()
  local df::DataFrame
  local b::IOBuffer

  iStream = open("mwe.csv")
  df = CSV.read(iStream)
  println("This is using open: $(df[:NAME])")
  close(iStream)

  iStream = open("mwe.csv")
  b = IOBuffer(read(iStream))
  df = CSV.read(b)
  println("This uses an IOBuffer: $(df[:NAME])")
  close(iStream)

  df = CSV.read("mwe.csv")
  println("This is direct: $(df[:NAME])")
end

#mwe()

function toyRegression()

  Random.seed!(11)
  years::Vector{Int} = collect(2001:2010) #10 years of data
  post::Vector{Int} = [zeros(Int, 5); ones(Int, 5)]
  groups::Vector{Int} = collect(1:100) #100 groups

  yearGroupFE::Vector{Float64} = # a fixed effect for each group and year with a time trend
    rand(length(groups) * length(years)) .*collect(1:(length(groups) * length(years))) ./ 100

  treated::Vector{Int} = [zeros(Int, 60); ones(Int, 40)]

  individualsPerGroup = 10
  individualFE::Vector{Float64} = rand(individualsPerGroup*length(groups))
  treatmentEffect = 2.0

  observationsN::Int = individualsPerGroup*length(groups)*length(years)
  individuals::Vector{Int} = collect(1:(individualsPerGroup*length(groups)))
  individualsN::Int = length(individuals)
  confounder::Vector{Float64} = 1.0 .* ( ones(observationsN) .* collect(1:observationsN) ./ 1000).^2 #time trend confounder
  ε::Vector{Float64} = rand(observationsN)

  irrelevantEffect::Vector{Float64} =  1.0 .* (-1.0 .* rand(observationsN) .* collect(1:observationsN) ./ 2000)
  irrelevant::Vector{Int} = [fill((i -> (i≤4)).(1:individualsPerGroup), length(groups)*length(years))...;]
  irrelevantEffect .*= irrelevant

  #this will hold the information
  df = DataFrame(N = collect(1:observationsN),
    ε = ε,
    outcome=Vector{MFloat64}(undef, observationsN),
    confounder=confounder,
    years=Vector{MSymbol}(undef, observationsN),
    groups=Vector{MSymbol}(undef, observationsN),
    individuals=Vector{MSymbol}(undef, observationsN),
    post=Vector{MInt}(undef, observationsN),
    irrelevant=irrelevant,
    treated=Vector{MInt}(undef, observationsN),
    keep = trues(observationsN)
    )

  ctr::Int = 0 #now construct the simulated sample
  yearGroupCtr::Int = 0
  for y::Int ∈ 1:length(years)
    for g::Int ∈ 1:length(groups)
      yearGroupCtr += 1
      for i::Int ∈ 1:individualsPerGroup
        individual::Int = (g-1)*individualsPerGroup + i #this is the index of the individual
        ctr+=1
        df[ctr, :years] = Symbol(years[y])
        df[ctr, :groups] = Symbol(groups[g])
        df[ctr, :individuals] = Symbol(individuals[individual])
        df[ctr, :post] = post[y]
        df[ctr, :treated] = treated[g]

        #determine if these are the folks we need to exclude in the post sample
        df[ctr, :keep] = !((irrelevant[ctr] == 1) && (post[y] == 1))

        #the actual process
        df[ctr, :outcome] = individualFE[individual] + yearGroupFE[yearGroupCtr] + confounder[ctr] +
          treated[g]*post[y]*treatmentEffect + irrelevantEffect[ctr] + ε[ctr]
      end
    end
  end

  df = df[df[:keep] .== true,:]
  df[:postXtreated]  = df[:treated] .* df[:post] #create interaction terms
  df[:yearXgroup] = ((y::Symbol,g::Symbol)->Symbol(y, "_", g)).(df[:years],df[:groups]) #make an interacted variable

  categorical!(df, :years) #make everything categorical
  categorical!(df, :groups)
  categorical!(df, :individuals)
  categorical!(df, :yearXgroup)

  models::Vector{CTLM} = Vector{CTLM}()
  #XNames::Vector{Vector{Symbol}} = Vector{Vector{Symbol}}() # holds the names of the x variables

  focal::Vector{Symbol} = [:intercept, :post, :postXtreated]
  push!(models, CTLM(df, Meta.parse("post + postXtreated + individuals"), :outcome,
    XNames=focal, YName=:outcome))

  push!(models, CTLM(df[df[:irrelevant] .≠ 1,:], Meta.parse("post + postXtreated + individuals"), :outcome,
    XNames=focal, YName=:outcome))

  focal = [:intercept, :post, :postXtreated, :confounder]
  push!(models, CTLM(df, Meta.parse("post + postXtreated + confounder + individuals"), :outcome,
    XNames=focal, YName=:outcome))

  push!(models, CTLM(df[df[:irrelevant] .≠ 1,:], Meta.parse("post + postXtreated + confounder + individuals"), :outcome,
    XNames=focal, YName=:outcome))

  displayed::Vector{Symbol} = [:post, :postXtreated, :confounder]

  tableInTex::String = tableText::String = texTable(models,
    getModWhiteΣ!,
    displayed,
    titleCaption = "Comparison",
    colNames = [["Just FE", "Dropped", "Just FE + confounder", "Dropped + confounder"]],
    contentRowNames = (string).(displayed),
    #descRowNames = descRowNames,
    #descContent = descContent,
    decimalDigits = 4,
    #stars=true,
    #starStrings = OVERRIDE_STAR_STRINGS,
    #clearMem = USE_AGGRESSIVE_GC,
    caption = "to be written")

  writeTables2File([tableInTex],
      HEADER_NAME, FOOTER_NAME, path=OUTPUT_PATH,
      outName = "comparisonTable.tex")

  #push!(models, CTLM(longDF, "",  :outcome, XNames=XNames[i], YName = YSpecs[i]))
end

@time toyRegression()

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


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

@test Base.Cartesian.exprresolve(:(1 + 3)) == 4
ex = Base.Cartesian.exprresolve(:(if 5 > 4; :x; else :y; end))
@test ex.args[2] == QuoteNode(:x)

@test Base.Cartesian.lreplace!("val_col", Base.Cartesian.LReplace{String}(:col, "col", 1)) == "val_1"

max_stars_repo_name: visr/julia

max_issues_repo_name: visr/julia

max_forks_repo_name: visr/julia


text:
"""
    x=cvp(y,R::UpperTriangular)

Solve the problem `argmin_x ||y-Rx||` for integer x using the technique from
the paper below. The input vector `y` is of length `n`, with upper
triangular `R` of dimension `n` by `n`, and the returned vector `x` of
length `n`.

    x=cvp(y,R,infinite=Val(true),Umin=-1,Umax=-Umin,nxMax=Int(ceil(log2(n)*1e6)))

If `infinite==Val(true)` then we search the (infinite) lattice, otherwise we
search integers in `[Umin,Umax]`. When `nxMax` is included, it gives a
maximum number of steps to take, otherwise a heuristic value is used.
Note that `cvp` does not handle complex numbers.

Follows "Faster Recursions in Sphere Decoding" Arash Ghasemmehdi, Erik
Agrell, IEEE Transactions on Information Theory, vol 57, issue 6 , June 2011.

# Examples
```jldoctest
julia> H=[1 2; 3 4]; Q,R=qr(H); uhat = cvp(Q'*[0,2],UpperTriangular(R))
2-element Vector{Float64}:
  2.0
 -1.0

julia> uhat = cvp(Q'*[0,2],UpperTriangular(R),Val(false),0,100)
2-element Vector{Float64}:
 1.0
 0.0

julia> n=100;H=randn(n,n); Q,Rt=qr(H); R=UpperTriangular(Rt);

julia> u=Int.(rand(0:1e10,n));y=H*u+rand(n)/100;

julia> uhat=cvp(Q'*y,R); sum(abs.(u-uhat))
0.0

julia> n=500;H=randn(n,n); Q,Rt=qr(H); R=UpperTriangular(Rt);

julia> u=Int.(rand([-1,1],n));y=H*u+rand(n)/100;

julia> uhat=cvp(Q'*y,R,Val(false),-1,1); sum(abs.(u-uhat))
0.0
```
"""
function cvp(r::AbstractArray{Td,1},G::UpperTriangular{Td,Array{Td,2}},
             infinite=Val(true),Umin=-1,Umax=-Umin,
             nxMax=typemax(Int)) where {Td<:Number}

    roundGA(x) = roundFinite(x,Td(Umin),Td(Umax))
    n = length(r)
    C = Inf
    i = n+1
    d = ones(Int,n)*n
    λ = zeros(n+1)
    F = zeros(Td,n,n)
    F[:,n] = r
    p = zeros(n)
    u = NaN*zeros(Int,n)
    û = NaN*zeros(Int,n)
    Δ = zeros(Int,n)

    if nxMax == typemax(Int)
        nxMax=Int(ceil(log2(n)*1e6)) # heuristic
    end
    nx=1

    begin
        @label LOOP

        nx+=1
        while λ[i]<C
            if i != 1
                i-=1
                for j=d[i]:-1:i+1
                    F[i,j-1] = F[i,j] - G[i,j]*u[j]
                end
                p[i] = F[i,i]/G[i,i]
                if infinite==Val(true)
                    u[i] = round(p[i])
                else
                    u[i] = roundGA(p[i])
                end
                y = (p[i]-u[i])*G[i,i]
                Δ[i] = sgns(y)
                λ[i] = λ[i+1]+y*y
            else
                û .= u
                C = λ[1]
            end
        end
        m = i

        while λ[i]≥C
            if i==n || nx>nxMax
                if nx>nxMax
                    @warn("terminated search after $(nx) iterations.")
                end
                return û
            else
                i+=1
                if infinite ≠ Val(true)
                    y = typemax(Td)
                end
                u[i]+=Δ[i]
                Δ[i]=-Δ[i]-sgns(Δ[i])
                if infinite ≠ Val(true)
                    if Umin ≤ u[i] ≤ Umax
                        y = (p[i]-u[i])*G[i,i]
                    else
                        u[i]+=Δ[i]
                        Δ[i]=-Δ[i]-sgns(Δ[i])
                        if Umin ≤ u[i] ≤ Umax
                            y = (p[i]-u[i])*G[i,i]
                        end
                    end
                else
                    y = (p[i]-u[i])*G[i,i]
                end
                λ[i] = λ[i+1]+y*y
            end
        end
        d[m:i-1] .= i
        for j=m-1:-1:1
            if d[j]<i
                d[j]=i
            else
                @goto LOOP
            end
        end
        @goto LOOP
    end
end

sgns(x) = x<=0 ? -one(x) : one(x)

function roundFinite(x,Umin,Umax)
    y=round(x)
    if x<Umin
        return Umin
    elseif x>Umax
        return Umax
    else
        return y
    end
end

"""
    x=svp(B)

Find the shortest basis vector `b` for the lattice formed by the matrix
`B`. This solves the 'shortest vector problem' (SVP). Note that `svp`
does not handle complex numbers.

Follows "Algorithm SHORTESTVECTOR(G)" from "Closest Point Search in
Lattices" by Erik Agrell, Thomas Eriksson, Alexander Vardy, and
Kenneth Zeger in IEEE Transactions on Information Theory, vol. 48, no. 8,
August 2002.

# Examples
```jldoctest
julia> H=[1 2; 3 4]; svp(H)
2-element Vector{Int64}:
 -1
 -1

julia> H = [1 0 0 0;   0 1 0 0;   208 175 663 0;     651 479 0  663];

julia> svp(H)
4-element Vector{Int64}:
  16
 -19
   3
 -11

```
"""
function svp(G::AbstractArray{Td,2}) where {Td<:Number}

    m,n = size(G)
    if n==1 return G[:,1];    end

    G2,_ = lll(G)
    # if qr gave positive values on diagonal we'd use  Q,G3 = qr(G2)
    G3 = cholesky(Hermitian(G2'*G2)).U
    H3 = inv(G3)
    #Q = G2*H3
    uhat = decodeSVPAgrell(H3)
    yhat = G2*uhat # the "shortest vector"
    return yhat
end

"""
    z=svpu(B)

Find the vector of integers z such that B*z is the shortest vector
in the lattice with basis B.

# Examples
```jldoctest
julia> H=[1 2; 3 4]; s=svp(H); u=LLLplus.svpu(H); [u s H*u]
2×3 Matrix{Int64}:
  1  -1  -1
 -1  -1  -1

```
"""
function svpu(G::AbstractArray{Td,2}) where {Td<:Number}
    m,n = size(G)
    if n==1 return [1]; end
    G2,T = lll(G)
    # if qr gave positive values on diagonal we'd use  Q,G3 = qr(G2)
    G3 = cholesky(Hermitian(G2'*G2)).U
    H3 = inv(G3)
    uhat = decodeSVPAgrell(H3)
    return T*uhat
end

"""
    x=decodeSVPAgrell(x,H)

Internal function.

Follows "Algorithm DECODE(G)" with "SHORTESTVECTOR" changes from the paper
"Closest Point Search in Lattices" by Erik Agrell, Thomas Eriksson,
Alexander Vardy, and Kenneth Zeger in IEEE Transactions on Information
Theory, vol. 48, no. 8, August 2002.

# Examples
```jldoctest
julia> H = inv([1.0 2; 0 4]);

julia> uhat = LLLplus.decodeSVPAgrell(H)
2-element Vector{Int64}:
 -1
  0

```
"""
function decodeSVPAgrell(H::AbstractArray{Td,2}) where {Td<:Number}

    n = size(H,1)
    bestdist=Inf
    k=n
    dist = zeros(n)
    e = zeros(n,n)
#    e[:,k] = H*x
    u = zeros(Int,n)
    û = zeros(Int,n)
    u[k] = round(e[k,k])
    y = (e[k,k]-u[k])/H[k,k]
    Δ = NaN*zeros(n)
    Δ[k] = sgns(y)
    begin
        @label LOOP
        newdist = dist[k]+y*y
        if newdist<bestdist
            if k!=1
                e[1:k-1,k-1]=e[1:k-1,k]-y*H[1:k-1,k]
                k -=1
                dist[k]=newdist
                u[k] = round(e[k,k])
                y = (e[k,k]-u[k])/H[k,k]
                Δ[k] = sgns(y)
            else
                if !(newdist≈0.0)
                    û[:] = u[:]
                    bestdist=newdist
                    k+=1
                end
                # û[:] = u[:]
                # bestdist=newdist
                # k+=1
                u[k] += Δ[k]
                y = (e[k,k]-u[k])/H[k,k]
                Δ[k] = -Δ[k] - sgns(Δ[k])
            end
        else
            if k==n
                return û
            else
                k+=1
                u[k]+= Δ[k]
                y = (e[k,k]-u[k])/H[k,k]
                Δ[k] = -Δ[k] - sgns(Δ[k])
            end
        end
        @goto LOOP
    end
end

max_stars_repo_name: christianpeel/LLLplus

max_issues_repo_name: christianpeel/LLLplus

max_forks_repo_name: christianpeel/LLLplus


text:
using DataDeps

register(DataDep(
  "communities-and-crime-2",
  "https://archive.ics.uci.edu/ml/datasets/Communities+and+Crime+Unnormalized",
  "https://archive.ics.uci.edu/ml/machine-learning-databases/00211/CommViolPredUnnormalizedData.txt",
  "383a9530b2802b4457986095b13af8f02a802a2504feaef45783f9c129e003f1",
  post_fetch_method=(path -> begin
    UCIData.process_dataset(path,
        target_index=146,
        feature_indices=6:129,
        id_indices=1,
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
@testset "169.majority-element.jl" begin
    @test majority_element([3, 2, 3]) == 3
    @test majority_element([2, 2, 1, 1, 1, 2, 2]) == 2
end


text:
module Oceanostics

#++++ TKEBudgetTerms exports
export TurbulentKineticEnergy, KineticEnergy
export IsotropicViscousDissipationRate, IsotropicPseudoViscousDissipationRate
export AnisotropicPseudoViscousDissipationRate
export XPressureRedistribution, YPressureRedistribution, ZPressureRedistribution
export XShearProduction, YShearProduction, ZShearProduction
#----

#++++ FlowDiagnostics exports
export RichardsonNumber, RossbyNumber
export ErtelPotentialVorticity, ThermalWindPotentialVorticity
export IsotropicBuoyancyMixingRate, AnisotropicBuoyancyMixingRate
export IsotropicTracerVarianceDissipationRate, AnisotropicTracerVarianceDissipationRate
#----

include("TKEBudgetTerms.jl")
include("FlowDiagnostics.jl")
include("progress_messengers.jl")

using .TKEBudgetTerms

end # module


text:
if "./src" in LOAD_PATH
    # print("in src")
else
    dir = "./src"
    push!(LOAD_PATH,dir)
    # print("not in src")
end

module bifurcations

struct bifurcation
    # Stokes Number Properties
    ro_d :: Float64     # Droplet density (ro_d) SI units: kg m-3
    d_d :: Float64      # Droplet diameter (dd) SI units: m
    d_0 :: Float64      # Tube inner diamater of parent branch (d_0)
    d_1 :: Float64      # Tube inner diameter of the daughter branch (d_1)
    mu :: Float64       # Carrier fluid dynamic viscosity (mu) SI units: N s m-2
    Q :: Float64        # Carrier fluid flow rate (Q) SI units: m3 s-1

    U :: Float64        # Fluid flow mean velocity (U) SI units: m s-1

    alpha :: Float64    # Angle between axes of parent and daughter branch
    phi :: Float64      # Bifurcation geometry SI units: degrees

    Stk :: Float64      # Stoke's Number (Stk)
    P :: Float64        # Efficiency of inertial deposition in an airway bifurcation

    function bifurcation(ro_d, d_d, d_0, d_1, mu, Q, alpha)
        U=4*Q/(pi*d_0^2)
        Stk = (ro_d*(d_d^2)*U)/(18*mu*d_1)

        phi = 2*(d_0/d_1)^2*sin(alpha)

        # print((sqrt( (4/3)*phi*Stk )))
        # print("\n")

        Deposition = (16/3*pi)*phi*Stk*(2-(sqrt( (4/3)*phi*Stk )))

        # print(Deposition)
        # print("\n")

        P = 1-Deposition
        # print(U)
        # print("\n")

        # print(Stk)
        # print("\n")

        # print(P)
        # print("\n")
    end
end

dataset1 = [5.4444E-07, 5.8472E-07, 6.2656E-07, 6.6727E-07, 7.2656E-07, 7.6727E-07, 8.3492E-07, 9.0000E-07, 9.6627E-07, 1.0557E-06, 1.1621E-06, 1.2616E-06, 1.3686E-06, 1.4746E-06, 1.5741E-06, 1.6736E-06, 1.7828E-06, 1.8863E-06, 1.9858E-06, 2.1567E-06, 2.3432E-06, 2.5136E-06, 2.6939E-06, 2.8643E-06, 3.0672E-06, 3.3113E-06, 3.5555E-06, 3.8230E-06]

# # bifurcation(ro_d, d_d, d_0, d_1, mu, Q, alpha)
# ro_d = 1000 # Droplet density of water
# # d_d = 6e-6 # droplet diameter of inhalation
# d_d = 3.8230E-06 # droplet diameter of inhalation
# d_0 = 0.018 # Diameter of trachea (https://www.google.com/url?sa=i&url=https%3A%2F%2Fwww.slideshare.net%2Fjinojustinj%2Ftracheal-pathologies&psig=AOvVaw38sp20tGz3PGVizR_R8F_N&ust=1605869019853000&source=images&cd=vfe&ved=0CAIQjRxqFwoTCICtjbi2ju0CFQAAAAAdAAAAABAT)
# d_1 = 0.0122 # Diameter of trachea (https://www.google.com/url?sa=i&url=https%3A%2F%2Fwww.slideshare.net%2Fjinojustinj%2Ftracheal-pathologies&psig=AOvVaw38sp20tGz3PGVizR_R8F_N&ust=1605869019853000&source=images&cd=vfe&ved=0CAIQjRxqFwoTCICtjbi2ju0CFQAAAAAdAAAAABAT)
# mu = 1.895e-5 # Dynamic viscosity of air (https://www.engineersedge.com/physics/viscosity_of_air_dynamic_and_kinematic_14483.htm#:~:text=The%20viscosity%20of%20air%20depends,10%2D5%20Pa%C2%B7s%20.&text=At%2025%20%C2%B0C%2C%20the,the%20kinematic%20viscosity%2015.7%20cSt.)
# # Q = 0.001 # Based on 60L/min flowrate of ventilators
# Q = 1e-4 # Based on 6L/min of average adult male (Guyton and Hall Textbook of Medical Physiology 12th Edition)
# # alpha = 1.06116 # SCA, all ages. Haskin Goodman et. al. "Normal Tracheal Bifurcation Angle"
# alpha = 1.040216 # SCA, all ages. Haskin Goodman et. al. "Normal Tracheal Bifurcation Angle"

# bifurcation(ro_d, d_d, d_0, d_1, mu, Q, alpha)
Q = 0.00050038/2 #Chapter 3
ro_d = 1000 # Droplet density of water
d_0 = 0.018 # Diameter of trachea (https://www.google.com/url?sa=i&url=https%3A%2F%2Fwww.slideshare.net%2Fjinojustinj%2Ftracheal-pathologies&psig=AOvVaw38sp20tGz3PGVizR_R8F_N&ust=1605869019853000&source=images&cd=vfe&ved=0CAIQjRxqFwoTCICtjbi2ju0CFQAAAAAdAAAAABAT)
d_1 = 0.018 # Diameter of daughter (https://www.google.com/url?sa=i&url=https%3A%2F%2Fwww.slideshare.net%2Fjinojustinj%2Ftracheal-pathologies&psig=AOvVaw38sp20tGz3PGVizR_R8F_N&ust=1605869019853000&source=images&cd=vfe&ved=0CAIQjRxqFwoTCICtjbi2ju0CFQAAAAAdAAAAABAT)
mu = 1.895e-5 # Dynamic viscosity of air (https://www.engineersedge.com/physics/viscosity_of_air_dynamic_and_kinematic_14483.htm#:~:text=The%20viscosity%20of%20air%20depends,10%2D5%20Pa%C2%B7s%20.&text=At%2025%20%C2%B0C%2C%20the,the%20kinematic%20viscosity%2015.7%20cSt.)
alpha = 1.040216 # IBA, all ages. Haskin Goodman et. al. "Normal Tracheal Bifurcation Angle"
# d_d = 3.8230E-06 # droplet diameter of inhalation

# bifurcation(ro_d, d_d, d_0, d_1, mu, Q, alpha)

post_bifurc1 = []

for d_d = dataset1
    bifurc_penetration = bifurcations.bifurcation(ro_d, d_d, d_0, d_1, mu, Q, alpha)
    push!(post_bifurc1, bifurc_penetration)
end

print(post_bifurc1)

end


text:
using RandomNumbers.Xorshifts: Xoshiro256StarStar

const _GLOBAL_RNG = Xoshiro256StarStar()
global_rng() = _GLOBAL_RNG

macro napi_call(env, sym)
    has_return = sym.head == :(::)
    status = gensym()
    func, result = if has_return
        sym.args[1], gensym()
    else
        sym, status
    end
    define_vars = if has_return
        return_type = sym.args[2]
        push!(func.args, :($result::Ptr{$return_type}))
        quote
          $result = Ref{$return_type}()
        end
    else
        quote end
    end
    insert!(func.args, 2, :($env::NodeCall.NapiEnv))
    fname = func.args[1]
    func.args[1] = :(:libjlnode.x)
    func.args[1].args[2] = QuoteNode(fname)
    esc(quote
        $define_vars
        $status = @ccall $func::NodeCall.NapiStatus
        if $status != NodeCall.NapiTypes.napi_ok
            @debug $status
            NodeCall.throw_error($env)
        end
        $result[]
    end)
end
macro napi_call(sym)
    env = gensym()
    esc(quote
        $env = node_env()
        @napi_call $env $sym
    end)
end

function throw_error(env::NapiEnv = node_env())
    info = @napi_call env napi_get_last_error_info()::Ptr{NapiExtendedErrorInfo}
    is_exception_pending = @napi_call env napi_is_exception_pending()::Bool
    err = if is_exception_pending
        @napi_call env napi_get_and_clear_last_exception()::NapiValue
    else
        error_info = unsafe_load(info)
        message = NapiValue(
            error_info.error_message == C_NULL ? "Error in native callback" : error_info.error_message
        )
        if error_info.error_code ∈ (
            NapiTypes.napi_object_expected,
            NapiTypes.napi_string_expected,
            NapiTypes.napi_boolean_expected,
            NapiTypes.napi_number_expected
        )
            @napi_call env napi_create_type_error(C_NULL::NapiValue, message::NapiValue)::NapiValue
        else
            @napi_call env napi_create_error(C_NULL::NapiValue, message::NapiValue)::NapiValue
        end
    end
    throw(NodeError(err))
end

const _SCOPE_STACK = Set{NapiHandleScope}()

function open_scope()
    if length(_SCOPE_STACK) < 2
        scope = @napi_call napi_open_handle_scope()::NapiHandleScope
        push!(_SCOPE_STACK, scope)
        scope
    end
end

close_scope(scope::Union{Nothing, NapiHandleScope}) = if !isnothing(scope)
    @napi_call napi_close_handle_scope(scope::NapiHandleScope)
    pop!(_SCOPE_STACK, scope)
end

function open_scope(f)
    scope = open_scope()
    try
        return f(scope)
    finally
        close_scope(scope)
    end
end
macro with_scope(code)
    esc(:(open_scope() do _
        $code
    end))
end

string_pointer(ptr::Ptr) = string("0x", string(UInt64(ptr), base = 16))
string_pointer(o) = string_pointer(pointer_from_objref(o))


text:
using Flux, NNlib, XLA, Test

xresult(f, args...) = collect(xla(f)(args...))

x, y = randn(2), randn(1, 2)
@test collect(xla(*)(x, y)) ≈ x*y

f(x) = gradient(x -> sum(x.*x), x)

@test collect(xla(f)([1, 2, 3])[1]) == [2, 4, 6]

f(x) = gradient(x -> sum(x.+x), x)

@test collect(xla(f)([1, 2, 3])[1]) == [2, 2, 2]

W = rand(2, 3)

f(x) = gradient(x -> sum(W*x), x)

@test collect(xla(f)([1.0, 2, 3])[1]) == f([1.0, 2, 3])[1]

@test xla(sin'')(0.5) == -sin(0.5)

x = rand(10)
m = Chain(Dense(10, 5), Dense(5, 2)) |> f64

f(x) = gradient(m -> sum(m(x)), m)[1]

xf = xla(f)

@test collect(xf(x).layers[1].W) ≈ f(x).layers[1].W

x = rand(10, 10, 5, 1)
w = rand(2, 2, 5, 3)
y = rand(size(conv(x, w))...)

@test xresult(conv, x, w) ≈ conv(x, w)

@test ∇conv_data(y, w, DenseConvDims(x, w)) ≈ xresult(y, w) do y, w
  ∇conv_data(y, w, DenseConvDims(x, w))
end

@test ∇conv_filter(x, y, DenseConvDims(x, w)) ≈ xresult(x, y) do x, y
  ∇conv_filter(x, y, DenseConvDims(x, w))
end

@test xresult(x -> maxpool(x, (2, 2)), x) ≈ maxpool(x, (2, 2))

y = maxpool(x, (2, 2))

@test xresult(x -> ∇maxpool(y, y, x, PoolDims(x, (2, 2))), x) ≈ ∇maxpool(y, y, x, PoolDims(x, (2, 2)))


text:
@doc raw"""
    SquaredEuclidean

The squared Euclidean function is defined by:

```math
f(\mathbf{x}, \mathbf{y}) = (\mathbf{x} - \mathbf{y})^{\intercal}(\mathbf{x} - \mathbf{y})
```
"""
struct SquaredEuclidean <: Metric end

@inline base_aggregate(::SquaredEuclidean, s::T, x::T, y::T) where {T} = s + (x-y)^2

@inline isstationary(::SquaredEuclidean) = true
@inline isisotropic(::SquaredEuclidean)  = true


text:
module GlobalEnergyGIS

export GISwind, GISsolar, GIShydro, GIStemp, GISturbines, makedistances, annualwindindex

using MAT, HDF5, ProgressMeter, Random, Interpolations, BenchmarkTools, Images,
    Statistics, DelimitedFiles, Dates, NCDatasets, JLD, Parameters, ImageSegmentation,
    StatsBase, CSV, Distances, Printf, TimeZones, DataFrames

include("rasterize_shapefiles.jl")
include("make_auxiliary_datasets.jl")
include("era5download.jl")
include("helperfunctions.jl")
include("make_regions.jl")
include("regiondefinitions.jl")
include("solarposition.jl")
include("makewindera5.jl")
include("makesolarera5.jl")
include("maketempera5.jl")
include("makedistances.jl")
include("GISwind.jl")
include("GISsolar.jl")
include("GIShydro.jl")
include("GIStemp.jl")
include("GISturbines.jl")
include("downloaddatasets.jl")
include("sspregiondefinitions.jl")
include("mapping.jl")
include("syntheticdemand_inputdata.jl")
include("syntheticdemand_training.jl")

function GISsequence(regionname, regionmat)
    saveregions(regionname, regionmat)
    makedistances(regionname)
    createmaps(regionname)
    GISsolar(gisregion=regionname, plotmasks=true)
    GISwind(gisregion=regionname, plotmasks=true)
    GIShydro(gisregion=regionname)
    predictdemand(gisregion=regionname, sspscenario="ssp2-26", sspyear=2050, era_year=2018)
end

end

max_stars_repo_name: henrikhodel/GlobalEnergyGIS

max_issues_repo_name: henrikhodel/GlobalEnergyGIS

max_forks_repo_name: henrikhodel/GlobalEnergyGIS