text:
"""
A generally useful event scheduler class.
Each instance of this class manages its own queue.
No multi-threading is implied; you are supposed to hack that
yourself, or use a single instance per application.
Each instance is parametrized with two functions, one that is
supposed to return the current time, one that is supposed to
implement a delay.    You can implement real-time scheduling by
substituting time and sleep from built-in module time, or you can
implement simulated time by writing your own functions.    This can
also be used to integrate scheduling with STDWIN events; the delay
function is allowed to modify the queue.    Time can be expressed as
integers or floating point numbers, as long as it is consistent.
Events are specified by tuples (time, priority, action, argument, kwargs).
As in UNIX, lower priority numbers mean higher priority; in this
way the queue can be maintained as a priority queue.    Execution of the
event means calling the action function, passing it the argument
sequence in "argument" (remember that in Python, multiple function
arguments are be packed in a sequence) and keyword parameters in "kwargs".
The action function may be an instance method so it
has another way to reference private data (besides global variables).
"""
module Sched

    export Scheduler, enter, enterabs, cancel, queue
    export FloatTimeFunc, UTCDateTimeFunc
    # Base exports: run, isempty
    import Base: run, isempty

    using DataStructures: PriorityQueue, peek, dequeue!, dequeue_pair!
    using Dates

    """
    Abstract type for struct that returns real-time or simulated time
    when called (functor)
    """
    abstract type TimeFunc end

    """
        UTCDateTimeFuncStruct()

    Functor that return real-time as DateTime (UTC) when called
    """
    struct UTCDateTimeFuncStruct <: TimeFunc
        func::Function
        args

        UTCDateTimeFuncStruct() = new(now, [Dates.UTC])
    end
    function (timefunc::UTCDateTimeFuncStruct)()
        timefunc.func(timefunc.args...)
    end
    UTCDateTimeFunc = UTCDateTimeFuncStruct()

    """
        FloatTimeFuncStruct()

    Functor that return real-time as Float when called
    """
    struct FloatTimeFuncStruct <: TimeFunc
        func::Function

        FloatTimeFuncStruct() = new(time)
    end
    function (timefunc::FloatTimeFuncStruct)()
        timefunc.func()
    end
    FloatTimeFunc = FloatTimeFuncStruct()

    """
    Default time function
    """
    global _time = UTCDateTimeFunc

    # Time as Float64
    #_time = FloatTimeFunc

    # Time as DateTime (UTC)
    #_time = UTCDateTimeFunc

    """
        Event(time_, priority, action, args...; kwargs...)

    Event structure

     - `time_`: Numeric type compatible with the return value of the timefunc function passed to the constructor.'
     - `priority`: Events scheduled for the same time will be executed in the order of their priority.
     - `action`: Executing the event means executing action(args...; kwargs...)
     - `args`: args is a sequence holding the positional arguments for the action.
     - `kwargs`: kwargs is a dictionary holding the keyword arguments for the action.

    """
    struct Event
        time_
        priority
        action
        args
        kwargs

        Event(time_, priority, action, args...; kwargs...) = new(time_, priority, action, args, kwargs)
    end
    run(event::Event) = event.action(event.args...; event.kwargs...)

    """
        Scheduler(; timefunc=_time, delayfunc=sleep)

    Initialize a new Scheduler instance, passing optionaly
    the time and delay functions

    The scheduler struct defines a generic interface to scheduling events.
    It needs two functions to actually deal with the “outside world”

    - The `timefunc` should be callable without arguments, and return a number (the “time”, in any units whatsoever). `timefunc` default is `UTCDateTimeFunc`.

    - The `delayfunc` function should be callable with one argument, compatible  with the output of `timefunc`, and should delay that many time units. delayfunc  will also be called with the argument 0 after each event is run to allow other threads an opportunity to run in multi-threaded applications.
    """
    struct Scheduler
        timefunc::TimeFunc
        delayfunc::Function

            _queue::PriorityQueue
            _lock::ReentrantLock

            function Scheduler(; timefunc=_time, delayfunc=sleep)
                q = PriorityQueue{Event, Priority}(Base.Order.Reverse)
                l = ReentrantLock()
                new(timefunc, delayfunc, q, l)
            end
    end

    """
        Priority(time_, priority)

    Priority of events
    """
    struct Priority
        time_
        priority
    end
    function Base.isless(p1::Sched.Priority, p2::Sched.Priority)
        if p1.time_ > p2.time_
            true
        elseif p1.time_ < p2.time_
            false
        else
            p1.priority > p2.priority
        end
    end

    """
        enterabs(sched, time_, priority, action, args...; kwargs...)

    Enter a new event in the queue at an absolute time.
    Returns an ID for the event which can be used to remove it,
    if necessary.
    """
    function enterabs(sched::Scheduler, time_, priority, action, args...; kwargs...)
        #println("enterabs $sched $time_ $priority $action with args=$args kwargs=$kwargs")
        event = Event(time_, priority, action, args...; kwargs...)
        l = sched._lock
        lock(l)
        try
            sched._queue[event] = Priority(time_, priority)
        finally
            unlock(sched._lock)
        end
        event # The ID
    end

    """
        enter(sched, delay, priority, action, args...; kwargs...)

    Enter a new event in the queue at a relative time.
    A variant of enterabs that specifies the time as a relative time.
    This is actually the more commonly used interface.
    """
    function enter(sched::Scheduler, delay, priority, action, args...; kwargs...)
        next_time = sched.timefunc() + delay
        event = enterabs(sched, next_time, priority, action, args...; kwargs...)
        event
    end

    """
        run(sched; blocking=true)

    Execute events until the queue is empty.
    If blocking is False executes the scheduled events due to
    expire soonest (if any) and then return the deadline of the
    next scheduled call in the scheduler.
    When there is a positive delay until the first event, the
    delay function is called and the event is left in the queue;
    otherwise, the event is removed from the queue and executed
    (its action function is called, passing it the argument).  If
    the delay function returns prematurely, it is simply
    restarted.
    It is legal for both the delay function and the action
    function to modify the queue or to raise an exception;
    exceptions are not caught but the scheduler's state remains
    well-defined so run() may be called again.
    A questionable hack is added to allow other threads to run:
    just after an event is executed, a delay of 0 is executed, to
    avoid monopolizing the CPU when other threads are also
    runnable.
    """
    function run(sched::Scheduler; blocking=true)
        l = sched._lock
        q = sched._queue
        delayfunc = sched.delayfunc
        timefunc = sched.timefunc
        while(true)
            lock(l)
            if length(q) == 0
                break
            end
            next_event, priority = peek(q)
            now_ = timefunc()
            if next_event.time_ > now_
                delay = true
            else
                delay = false
                event = dequeue!(q)
            end
            unlock(l)
            if delay
                if !blocking
                    return next_event.time_ - now_
                end
                delayfunc(next_event.time_ - now_)
            else
                run(event)
                delayfunc(0)     # Let other threads run
            end
        end
    end

    """
        cancel(sched, event)

    Remove an event from the queue.
    This must be presented the ID as returned by enter().
    If the event is not in the queue, this raises ValueError.
    """
    function cancel(sched::Scheduler, event::Event)
        lock(sched._lock)
        dequeue_pair!(sched._queue, event)
        unlock(sched._lock)
    end

    """
        isempty(sched) -> Bool

    Check whether the queue is empty.
    """
    function isempty(sched::Scheduler)
        lock(sched._lock)
        _isempty = isempty(sched._queue)
        unlock(sched._lock)
        _isempty
    end

    """
        queue(sched)

    Return an ordered list of upcoming events.
    """
    function queue(sched::Scheduler)
        q = sched._queue
        q2 = deepcopy(q)
        [dequeue!(q) for i in 1:length(q)]
    end

end # module


text:
module RulesNormalMeanPrecisionOutTest

using Test
using ReactiveMP
using Random

import ReactiveMP: @test_rules

@testset "rules:NormalMeanPrecision:out" begin

    @testset "Belief Propagation: (m_μ::PointMass, m_τ::PointMass)" begin

        @test_rules [ with_float_conversions = true ] NormalMeanPrecision(:out, Marginalisation) [
            (input = (m_μ = PointMass(-1.0), m_τ = PointMass(2.0)), output = NormalMeanPrecision(-1.0, 2.0)),
            (input = (m_μ = PointMass(1.0), m_τ = PointMass(2.0)),  output = NormalMeanPrecision(1.0, 2.0)),
            (input = (m_μ = PointMass(2.0), m_τ = PointMass(1.0)),  output = NormalMeanPrecision(2.0, 1.0))
        ]

    end

    @testset "Belief Propagation: (m_μ::UnivariateNormalDistributionsFamily, m_τ::PointMass)" begin

        @test_rules [ with_float_conversions = true ] NormalMeanPrecision(:out, Marginalisation) [
            (input = (m_μ = NormalMeanPrecision(0.0, 1.0),  m_τ = PointMass(2.0)), output = NormalMeanPrecision(0.0, 2.0 / 3.0)),
            (input = (m_μ = NormalMeanPrecision(-1.0, 1.0), m_τ = PointMass(1.5)), output = NormalMeanPrecision(-1.0, 0.6)),
            (input = (m_μ = NormalMeanPrecision(2.0, 0.5),  m_τ = PointMass(1.0)), output = NormalMeanPrecision(2.0, 1.0 / 3.0)),
        ]

        @test_rules [ with_float_conversions = true ] NormalMeanPrecision(:out, Marginalisation) [
            (input = (m_μ = NormalMeanVariance(0.0, 1.0),  m_τ = PointMass(2.0)), output = NormalMeanPrecision(0.0, 2.0 / 3.0)),
            (input = (m_μ = NormalMeanVariance(-1.0, 1.0), m_τ = PointMass(1.5)), output = NormalMeanPrecision(-1.0, 0.6)),
            (input = (m_μ = NormalMeanVariance(2.0, 0.5),  m_τ = PointMass(1.0)), output = NormalMeanPrecision(2.0, 2.0 / 3.0)),
        ]

        @test_rules [ with_float_conversions = true ] NormalMeanPrecision(:out, Marginalisation) [
            (input = (m_μ = NormalWeightedMeanPrecision(0.0, 1.0),  m_τ = PointMass(2.0)), output = NormalMeanPrecision(0.0, 2.0 / 3.0)),
            (input = (m_μ = NormalWeightedMeanPrecision(-1.0, 1.0), m_τ = PointMass(1.5)), output = NormalMeanPrecision(-1.0, 0.6)),
            (input = (m_μ = NormalWeightedMeanPrecision(2.0, 0.5),  m_τ = PointMass(1.0)), output = NormalMeanPrecision(4.0, 1.0 / 3.0)),
        ]

    end

    @testset "Variational: (m_μ::UnivariateNormalDistributionsFamily, q_τ::Any)" begin

        @test_rules [ with_float_conversions = true ] NormalMeanPrecision(:out, Marginalisation) [
            (input = (m_μ = NormalMeanPrecision(0.0, 1.0),  q_τ = GammaShapeRate(2.0, 1.0)), output = NormalMeanPrecision(0.0, 2.0 / 3.0)),
            (input = (m_μ = NormalMeanPrecision(-1.0, 1.0), q_τ = GammaShapeRate(3.0, 2.0)), output = NormalMeanPrecision(-1.0, 0.6)),
            (input = (m_μ = NormalMeanPrecision(2.0, 0.5),  q_τ = GammaShapeScale(10.0, 0.1)), output = NormalMeanPrecision(2.0, 1.0 / 3.0)),
        ]

        @test_rules [ with_float_conversions = true ] NormalMeanPrecision(:out, Marginalisation) [
            (input = (m_μ = NormalMeanVariance(0.0, 1.0),  q_τ = GammaShapeRate(2.0, 1.0)), output = NormalMeanPrecision(0.0, 2.0 / 3.0)),
            (input = (m_μ = NormalMeanVariance(-1.0, 1.0), q_τ = GammaShapeRate(3.0, 2.0)), output = NormalMeanPrecision(-1.0, 0.6)),
            (input = (m_μ = NormalMeanVariance(2.0, 0.5),  q_τ = GammaShapeScale(10.0, 0.1)), output = NormalMeanPrecision(2.0, 2.0 / 3.0)),
        ]

        @test_rules [ with_float_conversions = true ] NormalMeanPrecision(:out, Marginalisation) [
            (input = (m_μ = NormalWeightedMeanPrecision(0.0, 1.0),  q_τ = GammaShapeRate(2.0, 1.0)), output = NormalMeanPrecision(0.0, 2.0 / 3.0)),
            (input = (m_μ = NormalWeightedMeanPrecision(-1.0, 1.0), q_τ = GammaShapeRate(3.0, 2.0)), output = NormalMeanPrecision(-1.0, 0.6)),
            (input = (m_μ = NormalWeightedMeanPrecision(2.0, 0.5),  q_τ = GammaShapeScale(10.0, 0.1)), output = NormalMeanPrecision(4.0, 1.0 / 3.0)),
        ]

    end

    @testset "Variational: (q_μ::Any, q_τ::Any)" begin

        @test_rules [ with_float_conversions = true ] NormalMeanPrecision(:out, Marginalisation) [
            (input = (q_μ = PointMass(-1.0), q_τ = PointMass(2.0)), output = NormalMeanPrecision(-1.0, 2.0)),
            (input = (q_μ = PointMass(1.0), q_τ = PointMass(2.0)),  output = NormalMeanPrecision(1.0, 2.0)),
            (input = (q_μ = PointMass(2.0), q_τ = PointMass(1.0)),  output = NormalMeanPrecision(2.0, 1.0)),
        ]

        @test_rules [ with_float_conversions = true ] NormalMeanPrecision(:out, Marginalisation) [
            (input = (q_μ = NormalMeanVariance(-1.0, 2.0), q_τ = PointMass(2.0)), output = NormalMeanPrecision(-1.0, 2.0)),
            (input = (q_μ = NormalMeanPrecision(1.0, 4.0), q_τ = PointMass(3.0)),  output = NormalMeanPrecision(1.0, 3.0)),
            (input = (q_μ = NormalWeightedMeanPrecision(2.0, 4.0), q_τ = PointMass(1.0)),  output = NormalMeanPrecision(0.5, 1.0))
        ]

        @test_rules [ with_float_conversions = true ] NormalMeanPrecision(:out, Marginalisation) [
            (input = (q_μ = PointMass(-1.0), q_τ = Gamma(2.0, 1.0)), output = NormalMeanPrecision(-1.0, 2.0)),
            (input = (q_μ = PointMass(1.0), q_τ = Gamma(4.0, 2.0)),  output = NormalMeanPrecision(1.0, 8.0)),
            (input = (q_μ = PointMass(2.0), q_τ = Gamma(4.0, 6.0)),  output = NormalMeanPrecision(2.0, 24.0))
        ]

    end

end

end

ext: jl


text:
module BitBasis

include("utils.jl")
include("bit_str.jl")
include("operations.jl")
include("reorder.jl")
include("iterate_control.jl")

include("deprecations.jl")

end # module


text:
# ---
# title: 1515. Best Position for a Service Centre
# id: problem1515
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Geometry
# link: <https://leetcode.com/problems/best-position-for-a-service-centre/description/>
# hidden: true
# ---
#
# A delivery company wants to build a new service centre in a new city. The
# company knows the positions of all the customers in this city on a 2D-Map and
# wants to build the new centre in a position such that **the sum of the
# euclidean distances to all customers is minimum**.
#
# Given an array `positions` where `positions[i] = [xi, yi]` is the position of
# the `ith` customer on the map, return _the minimum sum of the euclidean
# distances_ to all customers.
#
# In other words, you need to choose the position of the service centre
# `[xcentre, ycentre]` such that the following formula is minimized:
#
# ![](https://assets.leetcode.com/uploads/2020/06/25/q4_edited.jpg)
#
# Answers within `10^-5` of the actual value will be accepted.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/06/25/q4_e1.jpg)
#
#
#
#     Input: positions = [[0,1],[1,0],[1,2],[2,1]]
#     Output: 4.00000
#     Explanation: As shown, you can see that choosing [xcentre, ycentre] = [1, 1] will make the distance to each customer = 1, the sum of all distances is 4 which is the minimum possible we can achieve.
#
#
# **Example 2:**
#
# ![](https://assets.leetcode.com/uploads/2020/06/25/q4_e3.jpg)
#
#
#
#     Input: positions = [[1,1],[3,3]]
#     Output: 2.82843
#     Explanation: The minimum possible sum of distances = sqrt(2) + sqrt(2) = 2.82843
#
#
# **Example 3:**
#
#
#
#     Input: positions = [[1,1]]
#     Output: 0.00000
#
#
# **Example 4:**
#
#
#
#     Input: positions = [[1,1],[0,0],[2,0]]
#     Output: 2.73205
#     Explanation: At the first glance, you may think that locating the centre at [1, 0] will achieve the minimum sum, but locating it at [1, 0] will make the sum of distances = 3.
#     Try to locate the centre at [1.0, 0.5773502711] you will see that the sum of distances is 2.73205.
#     Be careful with the precision!
#
#
# **Example 5:**
#
#
#
#     Input: positions = [[0,1],[3,2],[4,5],[7,6],[8,9],[11,1],[2,12]]
#     Output: 32.94036
#     Explanation: You can use [4.3460852395, 4.9813795505] as the position of the centre.
#
#
#
#
# **Constraints:**
#
#   * `1 <= positions.length <= 50`
#   * `positions[i].length == 2`
#   * `0 <= positions[i][0], positions[i][1] <= 100`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
# # SIMPLE CONTROLS

const EXTRACTOR_GIVEN_STR = Dict(
    "iterations"    => :(ic_model.n_iterations),
    "evaluation"    => :(ic_model.evaluation),
    "fitted_params" => :(fitted_params(ic_model.machine)),
    "report"        => :(report(ic_model.machine)),
    "machine"       => :(ic_model.machine),
    "model"         => :(ic_model.machine.model))

# maps "fitted_params" to ":WithFittedParamsDo":
_control_name(str) = string("With",
                            join(uppercasefirst.(split(str, "_"))),
                            "Do") |> Symbol

const NAME_GIVEN_STR =
    Dict([str=>_control_name(str) for str in keys(EXTRACTOR_GIVEN_STR)]...)

const DOC_GIVEN_STR = Dict(
    "iterations" =>
    "`x` is the current number of model iterations "*
    "(generally more than the number of control cycles)",

    "evaluation" =>
    "`x` is the latest performance evaluation, as returned by "*
    "`evaluate!(train_mach, resampling=..., ...)`. Not valid if "*
    "`resampling=nothing`",

    "fitted_params" =>
    "`x = fitted_params(mach)` is the fitted parameters "*
    "of the training machine, `mach`, in its current state",

    "report"       =>
    "`x = report(mach)` is the report associated with the training "*
    "machine, `mach`,  in its current state",

    "machine"      =>
    "`x` is the training machine in its current state",

    "model" =>
    "`x` is the model associated with the "*
    "training machine; `f` may mutate `x`, as in "*
    "`f(x) = (x.learning_rate *= 0.9)`")

for str in keys(EXTRACTOR_GIVEN_STR) # eg, "fitted_params"
    sym = Symbol(str)                # eg, :fitted_params
    C = NAME_GIVEN_STR[str]          # eg, :WithFittedParamsDo
    C_str = string(C)
    doc = DOC_GIVEN_STR[str]         # eg, "`x` is the training machine..."
    extractor = EXTRACTOR_GIVEN_STR[str]

    quote
        struct $C{F<:Function}
            f::F
            stop_if_true::Bool
            stop_message::Union{String,Nothing}
        end

        # constructor:
        $C(f::Function;
           stop_if_true=false,
           stop_message=nothing) = $C(f, stop_if_true, stop_message)
        $C(; f=x->@info("$($str): $x"), kwargs...) =
            $C(f, kwargs...)

        IterationControl.@create_docs(
        $C,
        header="$($C_str)(f=x->@info(\"$($str): \$x\"), "*
            "stop_if_true=false, "*
            "stop_message=nothing)",
        example="$($C_str)(x->put!(my_channel, x))",
        body="Call `f(x)`, where $($doc). "*
            "If `stop_if_true` is `true`, then trigger an early stop "*
            "if the value returned by `f` is `true`, logging the "*
            "`stop_message` if specified. ")

        function IterationControl.update!(c::$C,
                                          ic_model,
                                          verbosity,
                                          n,
                                          state...)
            x = $extractor
            r = c.f(x)
            done = (c.stop_if_true && r isa Bool && r) ? true : false
            return (done = done, $sym = x)
        end

        IterationControl.done(c::$C, state) = state.done

        function IterationControl.takedown(c::$C, verbosity, state)
            if state.done
                message = c.stop_message === nothing ?
                    "Stop triggered by a `$($C_str)` control. " :
                    c.stop_message
                verbosity > 0 && @info message
                return (done = true, log = message)
            else
                return (done = false, log = "")
            end
        end
    end |> eval
end

# # CYCLE LEARNING RATE

struct CycleLearningRate{F<:AbstractFloat}
    stepsize::Int
    lower::F
    upper::F
    learning_rate_parameter::Union{Symbol,Expr}
end
CycleLearningRate(; stepsize=4,
                  lower::T=0.001,
                  upper::T=0.05,
                  learning_rate_parameter=:(optimiser.η)) where T<:Real =
                      CycleLearningRate(stepsize,
                                        float(lower),
                                        float(upper),
                                        learning_rate_parameter)

IterationControl.@create_docs(
    CycleLearningRate,
    header="CycleLearningRate(stepsize=4, lower=0.001, "*
    "upper=0.05, learning_rate_parameter=:(optimiser.η))",
    example="CycleLearningRate(learning_rate_parameter=:learning_rate)",
    body="Mutate the specified `learning_rate_parameter` "*
    "(defaulting to the field appropriate for MLJFlux models) "*
    "to the next value "*
    "in a triangular cyclic learning rate policy, using the specified "*
    "`upper` and `lower` learning rate bounds. "*
    "See [L. N. Smith (2019)](https://ieeexplore.ieee.org/document/7926641): "*
    "\"Cyclical Learning Rates for Training "*
    "Neural Networks,\" 2017 IEEE Winter Conference on Applications "*
    "of Computer Vision (WACV), Santa Rosa, CA, USA, pp. 464-472.\n\n"*
    "Here `stepsize` is half the period of the mutation cycle, measured "*
    "in iterations of the underlying model. "*
    "For example, `stepsize=4` means the learning rate undergoes one "*
    "cycle of mutation every 8 epochs in an MLJFlux model.\n\n"*
    "**Note.** Since \"one iteration\" is the same as \"one epoch\" "*
    "in MLJFlux models, this means learning rate updates can be applied "*
    "once per epoch, at most, "*
    "rather than \"continuously\" throughout an epoch as in "*
    "the cited paper. ")

# return one cycle of learning rate values:
function one_cycle(c::CycleLearningRate)
    rise = range(c.lower, c.upper, length=c.stepsize + 1)
    fall = reverse(rise)
    return vcat(rise[1:end - 1], fall[1:end - 1])
end

function IterationControl.update!(control::CycleLearningRate,
                                  wrapper,
                                  verbosity,
                                  ncycles,
                                  state = (n = 0,))
    n = state.n
    rates = n == 0 ? one_cycle(control) : state.learning_rates
    index = mod(n, length(rates)) + 1
    r = rates[index]
    verbosity > 1 && @info "learning rate: $r"
    MLJBase.recursive_setproperty!(wrapper.model,
                                   control.learning_rate_parameter,
                                   r)
    return (n = n + 1, learning_rates = rates)
end


text:
using PyPlot
frequencies = range(16, 20000, length=10000)
fs = 96e3
function show_grid()
	minorticks_on()
	grid(which="major", color="gray", linestyle="-")
	grid(which="minor", color="lightgray", linestyle="--")
	xlabel("f in Hz")
end

function magplot(myfilter, frequencies=frequencies)
	plot(frequencies, 20*log10.(abs.(freqz(myfilter, frequencies, fs))), label="Level/dB")
	semilogx()
	show_grid()
end

function phaseplot(myfilter, frequencies=frequencies)
	plot(frequencies, angle.(freqz(myfilter, frequencies, fs)), label="Phase/rad")
	semilogx()
	show_grid()
end

function gdplot(myfilter, frequencies=frequencies)
	deltaf = frequencies.step.hi
	gdfrequencies = range(frequencies[1]+deltaf/2, frequencies[end]-deltaf/2, length=frequencies.len-1)
	phaseresponse = angle.(freqz(myfilter, frequencies, fs))|>unwrap
	gd = -diff(phaseresponse)/deltaf
	plot(gdfrequencies, gd)
	semilogx()
	show_grid()
end

function csvGDplot(filename)
	input = readdlm(filename, ';', skipstart=1)
	deltaf = [input[i+1,1]-input[i,1] for i=1:size(input,1)-1]
	gd = -diff(unwrap(input[:,3]/180*pi))./deltaf
	fAxis = [(input[i+1,1]+input[i,1])/2 for i=1:size(input,1)-1]
	plot(fAxis, gd)
	semilogx()
	show_grid()
end

"""
testplot(filename, myfilter)

plot csv (sep ';', skipstart=1, f/Hz, Level/dB, Phase/deg) level and phase
also calls magplot and phaseplot on given filter
"""
function testplot(filename, myfilter)
	input = readdlm(filename, ';', skipstart=1)
	plot(input[:,1], input[:,2], label="ref Level/dB")
	plot(input[:,1], input[:,3]/180*pi, label="ref Phase/rad")
	magplot(myfilter, input[:,1])
	phaseplot(myfilter, input[:,1])
	legend()
end


text:
export MinStepsGoal

"Goal specification where all steps have equal cost."
struct MinStepsGoal <: Goal
    terms::Vector{Term} # Goal terms to be satisfied
end

MinStepsGoal(problem::Problem) = MinStepsGoal(flatten_conjs(problem.goal))
MinStepsGoal(goal::Term) = MinStepsGoal(flatten_conjs(goal))

Base.hash(spec::MinStepsGoal, h::UInt) =
    hash(Set(spec.terms), h)
Base.isequal(s1::MinStepsGoal, s2::MinStepsGoal) =
    Set(s1.terms) == Set(s2.terms)

is_goal(spec::MinStepsGoal, domain::Domain, state::State) =
    satisfy(spec.terms, state, domain)[1]
is_violated(spec::MinStepsGoal, domain::Domain, state::State) = false
get_cost(spec::MinStepsGoal, ::Domain, ::State, ::Term, ::State) = 1
get_reward(spec::MinStepsGoal, ::Domain, ::State, ::Term, ::State) = -1
get_goal_terms(spec::MinStepsGoal) = spec.terms


text:
# Scheduler
# =========

struct Scheduler
    # global configuration
    config::Config

    # runtime statistics
    stats::Stats

    # trace for profiling
    trace::Trace

    # immediately available workers
    availworkers::Vector{Int}

    # waiting remote calls (key: worker, value: set of remote calls)
    waitcalls::Dict{Int,Set{Task}}

    # futures associated with remote calls
    futures::WeakKeyDict{Task,Future}

    # remote functions
    remotefuncs::Dict{Function,Dict{Int,RemoteChannel}}

    function Scheduler(config::Config)
        stats = Stats()
        workers = Int[]
        waitcalls = Dict{Int,Set{Task}}()
        for worker in config.workers
            for _ in 1:config.multiplicity
                push!(workers, worker)
            end
            waitcalls[worker] = Set{Task}()
        end
        scheduler = new(
            config, stats, Trace(),
            workers, waitcalls,
            WeakKeyDict{Task,Future}(), Dict())
        initstats!(stats, config.workers)
        return scheduler
    end
end

function Scheduler(workers::Vector{Int}=Base.workers())
    return Scheduler(Config(workers))
end

function teardown(scheduler)
    @tic
    for waitcalls in values(scheduler.waitcalls)
        while !isempty(waitcalls)
            wait(pop!(waitcalls))
        end
    end
    @toc scheduler.stats.teardowntime_ns
    finstats!(scheduler.stats)
    if scheduler.config.show_stats
        println(scheduler.config.output, sprint(show, scheduler.stats))
    end
    if scheduler.config.show_trace
        println(scheduler.config.output, scheduler.trace)
    end
    return nothing
end

function scoreworker(worker, scheduler, args)
    score = 0
    for arg in args
        if arg isa Future && arg.where == worker
            # data locality
            score += 1
        end
    end
    score -= length(scheduler.waitcalls[worker])
    return score
end

function selectworker(scheduler, args)
    @tic
    worker = 0
    while worker == 0
        while isempty(scheduler.availworkers)
            yield()
        end
        bestscore = typemin(Int)
        for w in scheduler.availworkers
            s = scoreworker(w, scheduler, args)
            if s > bestscore
                worker = w
                bestscore = s
            end
        end
    end
    @toc scheduler.stats.selecttime_ns
    return worker
end

function waitargs(args, stats)
    @tic
    newargs = map(a -> a isa RemoteCall ? wait(a) : a, args)
    @toc stats.waitargstime_ns
    return newargs
end

function waitarg(arg)
    if arg isa RemoteCall
        # Future
        return wait(arg)
    else
        return arg
    end
end

function call(scheduler, func, args...)
    @tic
    if elapsedtime(scheduler.trace) ≥ snapshot_interval(scheduler.config)
        takesnapshot!(scheduler.trace, scheduler.stats)
    end
    remotefunc = args -> func(map(fetch, args)...)
    if !haskey(scheduler.remotefuncs, remotefunc)
        scheduler.remotefuncs[remotefunc] = Dict{Int,RemoteChannel}()
    end
    while isempty(scheduler.availworkers)
        yield()
    end
    call = @task begin
        args = waitargs(args, scheduler.stats)
        worker = selectworker(scheduler, args)
        deleteat!(scheduler.availworkers, findfirst(scheduler.availworkers, worker))
        remotefuncs = scheduler.remotefuncs[remotefunc]
        cache = get!(remotefuncs, worker) do
            cache = RemoteChannel(worker)
            remotefuncs[worker] = cache
            put!(cache, remotefunc)
            cache
        end
        scheduler.stats.workerstats[worker].ncalls += 1
        call = current_task()
        waitcalls = scheduler.waitcalls[worker]
        push!(waitcalls, call)
        @tic
        future = remotecall_wait(args -> fetch(cache)(args), worker, args)
        @toc scheduler.stats.remotecalltime_ns
        scheduler.futures[call] = future
        delete!(waitcalls, call)
        push!(scheduler.availworkers, worker)
        future  # result
    end
    yield(call)
    @toc scheduler.stats.calltime_ns
    return RemoteCall(call)
end

function pmap(scheduler, func, iter)
    return [call(scheduler, func, item) for item in iter]
end


text:
using NamedArrays
using Compat
using Base.Test
using DataStructures

include("test.jl")


text:
context("HyperSphere") do
    @fact Sphere(Point(0,0,0), 4) --> HyperSphere{3,Int}(Point{3,Int}((0,0,0)),4)
    c = Circle(Point(0,0), 1)
    @fact c --> HyperSphere(Point(0,0), 1)
    @fact origin(c) --> Point(0,0)
    @fact isinside(c, 1, 2) --> false
    @fact isinside(c, 0.5, 0.5) --> true

    centered_rect = centered(HyperSphere)
    @fact centered_rect --> HyperSphere{3,Float32}(Point3f0(0), 0.5f0)
    centered_rect = centered(HyperSphere{2})
    @fact centered_rect --> HyperSphere{2,Float32}(Point2f0(0),0.5f0)

    centered_rect = centered(HyperSphere{2, Float64})
    @fact centered_rect --> HyperSphere{2,Float64}(Point(0.,0.), 0.5)

    centered_rect = centered(HyperSphere{3, Float32})
    @fact centered_rect --> HyperSphere{3,Float32}(Point3f0(0), 0.5f0)

    @fact widths(centered_rect) --> Vec3f0(1)
    @fact radius(centered_rect) --> 0.5f0
    @fact maximum(centered_rect) --> Vec3f0(0.5f0)
    @fact minimum(centered_rect) --> Vec3f0(-0.5f0)
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct ConnectionMonitorTestConfiguration <: SwaggerModel
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    testFrequencySec::Any # spec type: Union{ Nothing, Int32 } # spec name: testFrequencySec
    protocol::Any # spec type: Union{ Nothing, String } # spec name: protocol
    preferredIPVersion::Any # spec type: Union{ Nothing, String } # spec name: preferredIPVersion
    httpConfiguration::Any # spec type: Union{ Nothing, ConnectionMonitorHttpConfiguration } # spec name: httpConfiguration
    tcpConfiguration::Any # spec type: Union{ Nothing, ConnectionMonitorTcpConfiguration } # spec name: tcpConfiguration
    icmpConfiguration::Any # spec type: Union{ Nothing, ConnectionMonitorIcmpConfiguration } # spec name: icmpConfiguration
    successThreshold::Any # spec type: Union{ Nothing, ConnectionMonitorSuccessThreshold } # spec name: successThreshold

    function ConnectionMonitorTestConfiguration(;name=nothing, testFrequencySec=nothing, protocol=nothing, preferredIPVersion=nothing, httpConfiguration=nothing, tcpConfiguration=nothing, icmpConfiguration=nothing, successThreshold=nothing)
        o = new()
        validate_property(ConnectionMonitorTestConfiguration, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(ConnectionMonitorTestConfiguration, Symbol("testFrequencySec"), testFrequencySec)
        setfield!(o, Symbol("testFrequencySec"), testFrequencySec)
        validate_property(ConnectionMonitorTestConfiguration, Symbol("protocol"), protocol)
        setfield!(o, Symbol("protocol"), protocol)
        validate_property(ConnectionMonitorTestConfiguration, Symbol("preferredIPVersion"), preferredIPVersion)
        setfield!(o, Symbol("preferredIPVersion"), preferredIPVersion)
        validate_property(ConnectionMonitorTestConfiguration, Symbol("httpConfiguration"), httpConfiguration)
        setfield!(o, Symbol("httpConfiguration"), httpConfiguration)
        validate_property(ConnectionMonitorTestConfiguration, Symbol("tcpConfiguration"), tcpConfiguration)
        setfield!(o, Symbol("tcpConfiguration"), tcpConfiguration)
        validate_property(ConnectionMonitorTestConfiguration, Symbol("icmpConfiguration"), icmpConfiguration)
        setfield!(o, Symbol("icmpConfiguration"), icmpConfiguration)
        validate_property(ConnectionMonitorTestConfiguration, Symbol("successThreshold"), successThreshold)
        setfield!(o, Symbol("successThreshold"), successThreshold)
        o
    end
end # type ConnectionMonitorTestConfiguration

const _property_map_ConnectionMonitorTestConfiguration = Dict{Symbol,Symbol}(Symbol("name")=>Symbol("name"), Symbol("testFrequencySec")=>Symbol("testFrequencySec"), Symbol("protocol")=>Symbol("protocol"), Symbol("preferredIPVersion")=>Symbol("preferredIPVersion"), Symbol("httpConfiguration")=>Symbol("httpConfiguration"), Symbol("tcpConfiguration")=>Symbol("tcpConfiguration"), Symbol("icmpConfiguration")=>Symbol("icmpConfiguration"), Symbol("successThreshold")=>Symbol("successThreshold"))
const _property_types_ConnectionMonitorTestConfiguration = Dict{Symbol,String}(Symbol("name")=>"String", Symbol("testFrequencySec")=>"Int32", Symbol("protocol")=>"String", Symbol("preferredIPVersion")=>"String", Symbol("httpConfiguration")=>"ConnectionMonitorHttpConfiguration", Symbol("tcpConfiguration")=>"ConnectionMonitorTcpConfiguration", Symbol("icmpConfiguration")=>"ConnectionMonitorIcmpConfiguration", Symbol("successThreshold")=>"ConnectionMonitorSuccessThreshold")
Base.propertynames(::Type{ ConnectionMonitorTestConfiguration }) = collect(keys(_property_map_ConnectionMonitorTestConfiguration))
Swagger.property_type(::Type{ ConnectionMonitorTestConfiguration }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ConnectionMonitorTestConfiguration[name]))}
Swagger.field_name(::Type{ ConnectionMonitorTestConfiguration }, property_name::Symbol) =  _property_map_ConnectionMonitorTestConfiguration[property_name]

const _allowed_ConnectionMonitorTestConfiguration_protocol = ["Tcp", "Http", "Icmp"]

const _allowed_ConnectionMonitorTestConfiguration_preferredIPVersion = ["IPv4", "IPv6"]

function check_required(o::ConnectionMonitorTestConfiguration)
    (getproperty(o, Symbol("name")) === nothing) && (return false)
    true
end

function validate_property(::Type{ ConnectionMonitorTestConfiguration }, name::Symbol, val)
    if name === Symbol("protocol")
        Swagger.validate_param(name, "ConnectionMonitorTestConfiguration", :enum, val, _allowed_ConnectionMonitorTestConfiguration_protocol)
    end
    if name === Symbol("preferredIPVersion")
        Swagger.validate_param(name, "ConnectionMonitorTestConfiguration", :enum, val, _allowed_ConnectionMonitorTestConfiguration_preferredIPVersion)
    end
end


text:
module NFFT

using Printf
using Base.Cartesian
using FFTW
using Distributed
using SparseArrays
using LinearAlgebra
using FLoops
import Base.size
using Reexport
using SpecialFunctions: besseli, besselj

#using StaticArrays

@reexport using AbstractNFFTs

export NDFTPlan, NDCTPlan, NDSTPlan, NNDFTPlan,
       NFFTPlan, NFFTParams

#########################
# utility functions
#########################
include("utils.jl")
include("windowFunctions.jl")

#################################
# implementations
#################################

include("implementation.jl")
include("direct.jl")
include("precomputation.jl")

#################
# factory methods
#################

"""
        plan_nfft(x::Matrix{T}, N::NTuple{D,Int}, rest...;  kargs...)

compute a plan for the NFFT of a size-`N` array at the nodes contained in `x`.
"""
function AbstractNFFTs.plan_nfft(::Type{<:Array}, x::Matrix{T}, N::NTuple{D,Int}, rest...;
                   timing::Union{Nothing,TimingStats} = nothing, kargs...) where {T,D}
  t = @elapsed begin
    p = NFFTPlan(x, N, rest...; kargs...)
  end
  if timing != nothing
    timing.pre = t
  end
  return p
end

include("directional.jl")
include("multidimensional.jl")

function __init__()
  NFFT._use_threads[] = (Threads.nthreads() > 1)
  FFTW.set_num_threads(Threads.nthreads())
end

end


text:
using DataStructures

function parseinput(file)
  Rules = Dict{String, String}()
  lines = readlines(file)
  for line in lines[3:end]
    sp = split(line, " -> ")
    a, b = String(sp[1]), String(sp[2])
    Rules[a] = b
  end
  Rules, lines[1]
end

function day14(Rules::Dict{String, String}, S::String, steps::Int)
# function day14(Rules, S, steps)
  Pairs = DefaultDict{String, Int}(0)
  for i in 1:length(S)-1
    Pairs[S[i]*S[i+1]] += 1
  end
  for _ in 1:steps
    N = DefaultDict{String, Int}(0)
    for k in keys(Pairs)
      N[k[1]*Rules[k]] += Pairs[k]
      N[Rules[k]*k[2]] += Pairs[k]
    end
    Pairs = N
  end
  chars = DefaultDict{Char,Int}(0)
  for k in keys(Pairs)
    chars[k[1]] += Pairs[k]
  end
  chars[S[end]] += 1
  v = values(chars)
  maximum(v) - minimum(v)
end

Rules, S = parseinput("inputs/11-15/14.txt")
@show day14(Rules, S, 10)
@show day14(Rules, S, 40)

using BenchmarkTools
@btime day14(Rules, S, 40)


text:
## ----------------------------------------------------------------------------
# RBook

function find_obj(label0)
    book = currbook()
    for obj in eachobj(book)
        (get_label(obj) == label0) && return obj
    end
    return nothing
end

filter_book(qp, qps...) = filter_match(eachobj(currbook()), qp, qps...)

## ----------------------------------------------------------------------------
# bookbib
findall_bookbib(qp, qps...) = _findall_bookbib(bookdir(), qp, qps...)
findfirst_bookbib(qp, qps...) = _findfirst_bookbib(bookdir(), qp, qps...)
filter_bookbib(qp, qps...) = _filter_bookbib(bookdir(), qp, qps...)

## ----------------------------------------------------------------------------
# crossrefs
function findall_crossrefs(qp, qps...)
    vec = references(crossrefs())
    findall_match(vec, qp, qps...)
end

function findfirst_crossrefs(qp, qps...)
    vec = references(crossrefs())
    ret = findfirst_match(vec, qp, qps...)
    isnothing(ret) ? ret : last(ret)
end

function filter_crossrefs(qp, qps...)
    vec = references(crossrefs())
    filter_match(vec, qp, qps...)
end


text:
__precompile__(true)
module ODBC_v06

using DataStreams, Missings, CategoricalArrays, WeakRefStrings, DataFrames
import Compat: Sys

export Data, DataFrame, odbcdf

if VERSION < v"0.7.0-DEV.2575"
    const Dates = Base.Dates
else
    import Dates
end

include("API.jl")

"just a block of memory; T is the element type, `len` is total # of **bytes** pointed to, and `elsize` is size of each element"
mutable struct Block{T}
    ptr::Ptr{T}    # pointer to a block of memory
    len::Int       # total # of bytes in block
    elsize::Int    # size between elements in bytes
end

"""
Block allocator:
    -Takes an element type, and number of elements to allocate in a linear block
    -Optionally specify an extra dimension of elements that make up each element (i.e. container types)
"""
function Block(::Type{T}, elements::Int, extradim::Integer=1) where {T}
    len = sizeof(T) * elements * extradim
    block = Block{T}(convert(Ptr{T}, Libc.malloc(len)), len, sizeof(T) * extradim)
    finalizer(block, x->Libc.free(x.ptr))
    return block
end

# used for getting messages back from ODBC driver manager; SQLDrivers, SQLError, etc.
Base.string(block::Block, len::Integer) = String(transcode(UInt8, unsafe_wrap(Array, block.ptr, len, false)))

struct ODBCError <: Exception
    msg::String
end

const BUFLEN = 1024

function ODBCError(handle::Ptr{Void}, handletype::Int16)
    i = Int16(1)
    state = ODBC_v06.Block(ODBC_v06.API.SQLWCHAR, 6)
    native = Ref{ODBC_v06.API.SQLINTEGER}()
    error_msg = ODBC_v06.Block(ODBC_v06.API.SQLWCHAR, BUFLEN)
    msg_length = Ref{ODBC_v06.API.SQLSMALLINT}()
    while ODBC_v06.API.SQLGetDiagRec(handletype, handle, i, state.ptr, native, error_msg.ptr, BUFLEN, msg_length) == ODBC_v06.API.SQL_SUCCESS
        st  = string(state, 5)
        msg = string(error_msg, msg_length[])
        println("[ODBC] $st: $msg")
        i += 1
    end
    return true
end

#Macros to to check if a function returned a success value or not
macro CHECK(handle, handletype, func)
    str = string(func)
    esc(quote
        ret = $func
        ret != ODBC_v06.API.SQL_SUCCESS && ret != ODBC_v06.API.SQL_SUCCESS_WITH_INFO && ODBCError($handle, $handletype) &&
            throw(ODBCError("$($str) failed; return code: $ret => $(ODBC_v06.API.RETURN_VALUES[ret])"))
        nothing
    end)
end

Base.@deprecate listdrivers ODBC_v06.drivers
Base.@deprecate listdsns ODBC_v06.dsns

"List ODBC drivers that have been installed and registered"
function drivers()
    descriptions = String[]
    attributes   = String[]
    driver_desc = Block(ODBC_v06.API.SQLWCHAR, BUFLEN)
    desc_length = Ref{ODBC_v06.API.SQLSMALLINT}()
    driver_attr = Block(ODBC_v06.API.SQLWCHAR, BUFLEN)
    attr_length = Ref{ODBC_v06.API.SQLSMALLINT}()
    dir = ODBC_v06.API.SQL_FETCH_FIRST
    while ODBC_v06.API.SQLDrivers(ENV, dir, driver_desc.ptr, BUFLEN, desc_length, driver_attr.ptr, BUFLEN, attr_length) == ODBC_v06.API.SQL_SUCCESS
        push!(descriptions, string(driver_desc, desc_length[]))
        push!(attributes,   string(driver_attr, attr_length[]))
        dir = ODBC_v06.API.SQL_FETCH_NEXT
    end
    return [descriptions attributes]
end

"List ODBC DSNs, both user and system, that have been previously defined"
function dsns()
    descriptions = String[]
    attributes   = String[]
    dsn_desc    = Block(ODBC_v06.API.SQLWCHAR, BUFLEN)
    desc_length = Ref{ODBC_v06.API.SQLSMALLINT}()
    dsn_attr    = Block(ODBC_v06.API.SQLWCHAR, BUFLEN)
    attr_length = Ref{ODBC_v06.API.SQLSMALLINT}()
    dir = ODBC_v06.API.SQL_FETCH_FIRST
    while ODBC_v06.API.SQLDataSources(ENV, dir, dsn_desc.ptr, BUFLEN, desc_length, dsn_attr.ptr, BUFLEN, attr_length) == ODBC_v06.API.SQL_SUCCESS
        push!(descriptions, string(dsn_desc, desc_length[]))
        push!(attributes,   string(dsn_attr, attr_length[]))
        dir = ODBC_v06.API.SQL_FETCH_NEXT
    end
    return [descriptions attributes]
end

"""
A DSN represents an established ODBC connection.
It is passed to most other ODBC methods as a first argument
"""
mutable struct DSN
    dsn::String
    dbc_ptr::Ptr{Void}
    stmt_ptr::Ptr{Void}
    stmt_ptr2::Ptr{Void}
end

Base.show(io::IO,conn::DSN) = print(io, "ODBC_v06.DSN($(conn.dsn))")

const dsn = DSN("", C_NULL, C_NULL, C_NULL)

"""
Construct a `DSN` type by connecting to a valid ODBC DSN or by specifying a valid connection string.
Takes optional 2nd and 3rd arguments for `username` and `password`, respectively.
1st argument `dsn` can be either the name of a pre-defined ODBC DSN or a valid connection string.
A great resource for building valid connection strings is [http://www.connectionstrings.com/](http://www.connectionstrings.com/).
"""
function DSN(connectionstring::AbstractString, username::AbstractString=String(""), password::AbstractString=String(""); prompt::Bool=true)
    dbc = ODBC_v06.ODBCAllocHandle(ODBC_v06.API.SQL_HANDLE_DBC, ODBC_v06.ENV)
    dsns = ODBC_v06.dsns()
    found = false
    for d in dsns[:,1]
        connectionstring == d && (found = true)
    end
    if found
        @CHECK dbc ODBC_v06.API.SQL_HANDLE_DBC ODBC_v06.API.SQLConnect(dbc, connectionstring, username, password)
    else
        connectionstring = ODBCDriverConnect!(dbc, connectionstring, prompt)
    end
    stmt = ODBCAllocHandle(ODBC_v06.API.SQL_HANDLE_STMT, dbc)
    stmt2 = ODBCAllocHandle(ODBC_v06.API.SQL_HANDLE_STMT, dbc)
    global dsn
    dsn.dsn = connectionstring
    dsn.dbc_ptr = dbc
    dsn.stmt_ptr = stmt
    dsn.stmt_ptr2 = stmt2
    return DSN(connectionstring, dbc, stmt, stmt2)
end

"disconnect a connected `DSN`"
function disconnect!(conn::DSN)
    ODBCFreeStmt!(conn.stmt_ptr)
    ODBCFreeStmt!(conn.stmt_ptr2)
    ODBC_v06.API.SQLDisconnect(conn.dbc_ptr)
    return nothing
end

mutable struct Statement
    dsn::DSN
    stmt::Ptr{Void}
    query::String
    task::Task
end

"An `ODBC_v06.Source` type executes a `query` string upon construction and prepares data for streaming to an appropriate `Data.Sink`"
mutable struct Source{T} <: Data.Source
    schema::Data.Schema
    dsn::DSN
    query::String
    columns::T
    status::Int
    rowsfetched::Ref{ODBC_v06.API.SQLLEN}
    rowoffset::Int
    boundcols::Vector{Any}
    indcols::Vector{Vector{ODBC_v06.API.SQLLEN}}
    sizes::Vector{ODBC_v06.API.SQLULEN}
    ctypes::Vector{ODBC_v06.API.SQLSMALLINT}
    jltypes::Vector{Type}
    supportsreset::Bool
end

Base.show(io::IO, source::Source) = print(io, "ODBC_v06.Source:\n\tDSN: $(source.dsn)\n\tstatus: $(source.status)\n\tschema: $(source.schema)")

include("Source.jl")
include("Sink.jl")
include("sqlreplmode.jl")

function __init__()
    global const ENV = ODBC_v06.ODBCAllocHandle(ODBC_v06.API.SQL_HANDLE_ENV, ODBC_v06.API.SQL_NULL_HANDLE)
    global const DECZERO = Float64(0)
    toggle_sql_repl()
end

# used to 'clear' a statement of bound columns, resultsets,
# and other bound parameters in preparation for a subsequent query
function ODBCFreeStmt!(stmt)
    ODBC_v06.API.SQLFreeStmt(stmt, ODBC_v06.API.SQL_CLOSE)
    ODBC_v06.API.SQLFreeStmt(stmt, ODBC_v06.API.SQL_UNBIND)
    ODBC_v06.API.SQLFreeStmt(stmt, ODBC_v06.API.SQL_RESET_PARAMS)
end

end #ODBC module


text:
# This example is reused in the documentation
# -> NLS_Solver.jl/docs/src/nonlinear_regression.md
#
using NLS_Solver
using Random
using Plots

# Be reproducible
rng = MersenneTwister(1234)

# Define model ================
#
function gaussian_peak(x::Real,
                       h::Real,
                       μ::Real,
                       σ::Real)
    @assert σ>0

    h*exp(-((x-μ)/σ)^2/2)
end

function exp_baseline(x::Real,
                      b::Real,
                      τ::Real)
    @assert τ>0

    b*exp(-x/τ)
end

function model(x::Real,θ::AbstractVector)
    @assert length(θ) == 5

    # one gaussian peaks + exp baseline
    #
    (h, μ, σ, b, τ) = (θ_i for θ_i in θ)

    exp_baseline(x,b,τ) + gaussian_peak(x,h,μ,σ)
end

function residue(X::AbstractVector,Y::AbstractVector,θ::AbstractVector)
    map(zip(X,Y)) do (X_i,Y_i)
        Y_i - model(X_i, θ)
    end
end

# Generate synthetic data ================
#
X=[1:0.1:30;]
θ_true=[2.0,15,3.0,3.0,8.0]

Y_true = map(X) do X_i
    model(X_i,θ_true)
end
Y_noisy = map(Y_true) do Y_i
    Y_i + 0.5 * randn()
end

plot(X,Y_noisy,label="noisy signal")
plot!(X,Y_true,linewidth=3,linestyle=:dot,label="true signal")

# Solve the problem ================
#
n_θ = length(θ_true)
n_sample = length(X)

# Bound constraints ----------------
#
#                         h,  μ,  σ,   b,  τ
θ_lowerbound = Float64[   1, 10,  1,   0,  1 ]
θ_upperbound = Float64[ Inf, 20, 10, Inf, 10 ]
bc = BoundConstraints(θ_lowerbound, θ_upperbound)

# Initial guess
#
θ_init = project!(ones(n_θ),bc)

# Solver ----------------
#
conf = LevenbergMarquardt_BC_Conf()

# Wrap residue ----------------
#
nls = create_NLS_problem_using_ForwardDiff(θ->residue(X,Y_noisy,θ),n_θ => n_sample)

# Solve it! ----------------
#
result = solve(nls,θ_init,bc,conf)

# Plot fitted model ================
#

# Initial model ----------------
#
# Y_init = map(X_i -> model(X_i,θ_init), X)
# plot!(X,Y_init,linewidth=3,label="initial model")

# Fitted model ----------------
#
θ_solution = solution(result)

Y_solution = map(X_i -> model(X_i,θ_solution), X)
plot!(X,Y_solution,linewidth=3,label="fitted model")


text:
module DiffEqSensitivity

using DiffEqBase, ForwardDiff, Tracker, FiniteDiff, Statistics
using DiffEqCallbacks, QuadGK, RecursiveArrayTools, LinearAlgebra
using DiffEqOperators
using Adapt
using LinearSolve
using Parameters: @unpack
using StochasticDiffEq
import DiffEqNoiseProcess
import RandomNumbers: Xorshifts
using Random
import ZygoteRules, Zygote, ReverseDiff
import ArrayInterfaceCore, ArrayInterfaceTracker
import Enzyme
import GPUArrays

using Cassette, DiffRules
using Core: CodeInfo, SlotNumber, SSAValue, ReturnNode, GotoIfNot
using EllipsisNotation

using Markdown

using Reexport
import ChainRulesCore: @thunk, NoTangent, @not_implemented
abstract type SensitivityFunction end
abstract type TransformedFunction end

include("hasbranching.jl")
include("sensitivity_algorithms.jl")
include("derivative_wrappers.jl")
include("sensitivity_interface.jl")
include("forward_sensitivity.jl")
include("adjoint_common.jl")
include("lss.jl")
include("nilss.jl")
include("nilsas.jl")
include("backsolve_adjoint.jl")
include("interpolating_adjoint.jl")
include("quadrature_adjoint.jl")
include("callback_tracking.jl")
include("concrete_solve.jl")
include("second_order.jl")
include("steadystate_adjoint.jl")
include("sde_tools.jl")

# AD Extensions
include("reversediff.jl")
include("tracker.jl")
include("zygote.jl")

export extract_local_sensitivities

export ODEForwardSensitivityFunction, ODEForwardSensitivityProblem, SensitivityFunction,
       ODEAdjointSensitivityProblem, ODEAdjointProblem, AdjointSensitivityIntegrand,
       SDEAdjointProblem, RODEAdjointProblem, SensitivityAlg,
       adjoint_sensitivities, adjoint_sensitivities_u0,
       ForwardLSSProblem, AdjointLSSProblem,
       NILSSProblem, NILSASProblem,
       shadow_forward, shadow_adjoint

export BacksolveAdjoint, QuadratureAdjoint, InterpolatingAdjoint,
       TrackerAdjoint, ZygoteAdjoint, ReverseDiffAdjoint,
       ForwardSensitivity, ForwardDiffSensitivity,
       ForwardDiffOverAdjoint,
       SteadyStateAdjoint,
       ForwardLSS, AdjointLSS, NILSS, NILSAS

export second_order_sensitivities, second_order_sensitivity_product

export TrackerVJP, ZygoteVJP, EnzymeVJP, ReverseDiffVJP

export StochasticTransformedFunction

end # module


text:
#Look up Table
include("lut/mt.jl")

"""
    tetIx(tIx, cubeindex)

Determines which case in the triangle table we are dealing with
"""
@inline function tetIx(tIx, cubeindex)
    @inbounds v1 = subTetsMask[tIx][1]
    @inbounds v2 = subTetsMask[tIx][2]
    ix = 0x01 & cubeindex | (0x40 & cubeindex) >> 0x03
    !iszero(v1 & cubeindex) && (ix |= 0x02)
    !iszero(v2 & cubeindex) && (ix |= 0x04)
    ix
end

"""
    vertId(e, x, y, z, nx, ny)

Determines a unique integer ID associated with the edge. This is used
as a key in the vertex dictionary. It needs to be both unambiguous (no
two edges get the same index) and unique (every edge gets the same ID
regardless of which of its neighboring voxels is asking for it) in order
for vertex sharing to be implemented properly.
"""
@inline function vertId(e, x, y, z, nx, ny)
    dx, dy, dz = voxCrnrPosInt[voxEdgeCrnrs[e][1]]
    voxEdgeDir[e]+7*(x-0x01+dx+nx*(y-0x01+dy+ny*(z-0x01+dz)))
end

"""
    vertPos(e, x, y, z, vals::NTuple{8,T}, iso::Real, eps::Real, ::Type{VertType}) where {T<:Real, VertType}

Assuming an edge crossing, determines the point in space at which it
occurs.
eps represents the "bump" factor to keep vertices away from voxel
corners (thereby preventing degeneracies).
"""
@inline function vertPos(e, x, y, z, scale, origin, vals::V, iso, eps, ::Type{VertType}) where {V, VertType}
    T = eltype(vals)

    ixs     = voxEdgeCrnrs[e]
    srcVal  = vals[ixs[1]]
    tgtVal  = vals[ixs[2]]
    a       = min(max((iso-srcVal)/(tgtVal-srcVal), eps), one(T)-eps)
    b       = one(T)-a
    c1 = voxCrnrPos(VertType)[ixs[1]]
    c2 = voxCrnrPos(VertType)[ixs[2]]

    ((VertType(x,y,z) + c1 .* b + c2.* a) .- 1) .* scale .+ origin
end

"""
    getVertId(e, x, y, z, nx, ny,
                vals, iso::Real,
                vts::Dict,
                vtsAry::Vector,
                eps::Real)

Gets the vertex ID, adding it to the vertex dictionary if not already
present.
"""
@inline function getVertId(e, x, y, z, nx, ny,
                           vals, iso::Real,
                           scale, origin,
                           vts::Dict,
                           vtsAry::Vector,
                           eps::Real,
                           reduceverts::Bool)

    VertType = eltype(vtsAry)
    if reduceverts
        vId = vertId(e, x, y, z, nx, ny)
        haskey(vts, vId) && return vts[vId]
    end

    # calculate vert position
    v = vertPos(e, x, y, z, scale, origin, vals, iso, eps, VertType)
    push!(vtsAry, v)

    # if deduplicting, push to dict
    if reduceverts
        vts[vId] = length(vtsAry)
    end

    return length(vtsAry)
end

"""
    voxEdgeId(subTetIx, tetEdgeIx)
Given a sub-tetrahedron case and a tetrahedron edge ID, determines the
corresponding voxel edge ID.
"""
@inline function voxEdgeId(subTetIx, tetEdgeIx)
    srcVoxCrnr = subTets[subTetIx][tetEdgeCrnrs[tetEdgeIx][1]]
    tgtVoxCrnr = subTets[subTetIx][tetEdgeCrnrs[tetEdgeIx][2]]
    return voxEdgeIx[srcVoxCrnr][tgtVoxCrnr]
end

"""
    procVox(vals, iso::Real, x, y, z, nx, ny,
                    vts::Dict, vtsAry::Vector, fcs::Vector,
                    eps::Real)

Processes a voxel, adding any new vertices and faces to the given
containers as necessary.
"""
function procVox(vals, iso::Real, x, y, z, nx, ny, scale, origin,
                 vts::Dict, vtsAry::Vector, fcs::Vector,
                 eps::Real, cubeindex, reduceverts)
    VertType = eltype(vtsAry)
    FaceType = eltype(fcs)
    # check each sub-tetrahedron in the voxel
    @inbounds for i = 1:6
        tIx = tetIx(i, cubeindex)
        (tIx == 0x00 || tIx == 0x0f) && continue

        e = tetTri[tIx]

        # add the face to the list
        push!(fcs, FaceType(
                    getVertId(voxEdgeId(i, e[1]), x, y, z, nx, ny, vals, iso, scale, origin, vts, vtsAry, eps, reduceverts),
                    getVertId(voxEdgeId(i, e[2]), x, y, z, nx, ny, vals, iso, scale, origin, vts, vtsAry, eps, reduceverts),
                    getVertId(voxEdgeId(i, e[3]), x, y, z, nx, ny, vals, iso, scale, origin, vts, vtsAry, eps, reduceverts)))

        # bail if there are no more faces
        iszero(e[4]) && continue
        push!(fcs, FaceType(
                    getVertId(voxEdgeId(i, e[4]), x, y, z, nx, ny, vals, iso, scale, origin, vts, vtsAry, eps, reduceverts),
                    getVertId(voxEdgeId(i, e[5]), x, y, z, nx, ny, vals, iso, scale, origin, vts, vtsAry, eps, reduceverts),
                    getVertId(voxEdgeId(i, e[6]), x, y, z, nx, ny, vals, iso, scale, origin, vts, vtsAry, eps, reduceverts)))
    end
end

function isosurface(sdf::AbstractArray{T, 3}, method::MarchingTetrahedra, ::Type{VertType}=SVector{3,Float64}, ::Type{FaceType}=SVector{3, Int};
                    origin=VertType(-1,-1,-1), widths=VertType(2,2,2)) where {T, VertType, FaceType}

    vts    = Dict{Int, Int}()
    fcs    = FaceType[]
    vtsAry = VertType[]
    sizehint!(vts, div(length(sdf),8))
    sizehint!(vtsAry, div(length(sdf),8))
    sizehint!(fcs, div(length(sdf),4))
    # process each voxel
    scale = widths ./ VertType(size(sdf) .- 1)
    nx::Int, ny::Int, nz::Int = size(sdf)

    @inbounds for k = 1:nz-1, j = 1:ny-1, i = 1:nx-1

        vals = (sdf[i,  j,  k  ],
                sdf[i,  j+1,k  ],
                sdf[i+1,j+1,k  ],
                sdf[i+1,j,  k  ],
                sdf[i,  j,  k+1],
                sdf[i,  j+1,k+1],
                sdf[i+1,j+1,k+1],
                sdf[i+1,j,  k+1])

        cubeindex = method.insidepositive ? _get_cubeindex_pos(vals, method.iso) : _get_cubeindex(vals, method.iso)

        _no_triangles(cubeindex) && continue

        procVox(vals, method.iso, i, j, k, nx, ny, scale, origin, vts, vtsAry, fcs, method.eps, cubeindex, method.reduceverts)
    end

    vtsAry,fcs
end

function isosurface(f::Function, method::MarchingTetrahedra,
                    ::Type{VertType}=SVector{3,Float64}, ::Type{FaceType}=SVector{3, Int};
                    origin=VertType(-1,-1,-1), widths=VertType(2,2,2),
                    samples::NTuple{3,T}=_DEFAULT_SAMPLES) where {T, VertType, FaceType}

    vts    = Dict{Int, Int}()
    fcs    = FaceType[]
    vtsAry = VertType[]
    sizehint!(vts, div(prod(samples),8))
    sizehint!(vtsAry, div(prod(samples),8))
    sizehint!(fcs, div(prod(samples),4))
    # process each voxel
    scale = widths ./ (VertType(samples...) .-1)
    nx::Int, ny::Int, nz::Int = samples
    zv = zero(eltype(VertType))
    vals = (zv,zv,zv,zv,zv,zv,zv,zv)

    @inbounds for k = 1:nz, j = 1:ny, i = 1:nx
        points = (VertType(i-1,j-1,k-1).* scale .+ origin,
                  VertType(i-1,j  ,k-1).* scale .+ origin,
                  VertType(i  ,j  ,k-1).* scale .+ origin,
                  VertType(i  ,j-1,k-1).* scale .+ origin,
                  VertType(i-1,j-1,k  ).* scale .+ origin,
                  VertType(i-1,j  ,k  ).* scale .+ origin,
                  VertType(i  ,j  ,k  ).* scale .+ origin,
                  VertType(i  ,j-1,k  ).* scale .+ origin)
        if i == 1
            vals = (f(points[1]),
                    f(points[2]),
                    f(points[3]),
                    f(points[4]),
                    f(points[5]),
                    f(points[6]),
                    f(points[7]),
                    f(points[8]))
        else
            vals = (vals[4],
                    vals[3],
                    f(points[3]),
                    f(points[4]),
                    vals[8],
                    vals[7],
                    f(points[7]),
                    f(points[8]))
        end

        cubeindex = method.insidepositive ? _get_cubeindex_pos(vals, method.iso) : _get_cubeindex(vals, method.iso)

        _no_triangles(cubeindex) && continue

        procVox(vals, method.iso, i, j, k, nx, ny, scale, origin, vts, vtsAry, fcs, method.eps, cubeindex, method.reduceverts)
    end

    vtsAry,fcs
end


text:
module Graphs

using DocStringExtensions
using UnPack
using LinearAlgebra
using ..Mesh

export AbstractGraph,
    Graph, build_face_graph, build_vertex_graph, laplacian_x_v!, fiedler_vector, split_graph

abstract type AbstractGraph end
struct Lanczos end
struct PowerIt end
struct Direct end
"""
    Graph{I,FT,IA1D} <: AbstractGraph

Connectivity graph for a mesh.

# Fields
$(DocStringExtensions.FIELDS)
"""
struct Graph{I<:Int,FT<:AbstractFloat,IA1D<:AbstractVector{I}} <: AbstractGraph
    "number of graph vertices"
    nverts::I
    "edge data for each vertex"
    edge_data::IA1D
    "edge data offsets for each vertex"
    edge_offset::IA1D
    "vertex numbers"
    vertices::IA1D
    "spectral radius estimate/bound of graph laplacian matrix"
    spectral_radius::FT
end

# constructor for an empty graph
Graph(::Type{I}, ::Type{FT}) where {I,FT} = Graph{I,FT,Vector{I}}(0, I[], I[], I[], FT(0))

"""
    build_vertex_graph(mesh::AbstractMesh{I,FT}) where {I,FT}

This function builds a graph, with adjacency based on mesh elements connected by a vertex.
"""
function build_vertex_graph(mesh::AbstractMesh{I,FT}) where {I,FT}
    @unpack nelems, unique_verts, uverts_conn, uverts_offset, elem_verts = mesh
    connect = [zeros(I, i) for i in zeros(I, nelems)]
    gsize = 0
    spectral_radius = FT(0)
    for i = 1:length(unique_verts)
        st, en = uverts_offset[i:i+1]
        for j = st:2:en-1
            el = uverts_conn[j]
            for k = st:2:en-1
                elc = uverts_conn[k]
                if elc ≠ el && findfirst(connect[el] .== elc) == nothing
                    push!(connect[el], elc)
                    gsize += 1
                end
            end
        end
    end
    # construct the graph object
    g_nverts = nelems
    g_edge_offset = Array{I}(undef, g_nverts + 1)
    g_edge_offset[1] = 1
    g_edge_data = Array{I}(undef, gsize)

    for el = 1:nelems
        sort!(connect[el])
        st = g_edge_offset[el]
        en = st + length(connect[el])
        g_edge_offset[el+1] = en
        g_edge_data[st:en-1] .= connect[el]
        spectral_radius = max(spectral_radius, FT(en - st))
    end
    spectral_radius *= 2
    return Graph(g_nverts, g_edge_data, g_edge_offset, Array(1:g_nverts), spectral_radius)
end
"""
    build_face_graph(mesh::AbstractMesh{I,FT}) where {I,FT}

This function builds a graph, with adjacency based on mesh elements connected by a face.
"""
function build_face_graph(mesh::AbstractMesh{I,FT}) where {I,FT}
    @unpack nelems, elem_faces, face_neighbors = mesh
    nlocalfaces = size(elem_faces, 2)

    g_nverts = nelems
    g_edge_offset = Array{I}(undef, g_nverts + 1)
    g_edge_offset[1] = 1
    g_edge_data = Array{I}(undef, nlocalfaces * nelems) # conservative initial estimate

    temp = zeros(I, nlocalfaces)
    spectral_radius = FT(0)

    for el = 1:nelems
        n, st = 0, g_edge_offset[el]
        for f = 1:nlocalfaces
            fc = elem_faces[el, f]
            if face_neighbors[fc, 1] == el
                neigh = face_neighbors[fc, 3]
            elseif face_neighbors[fc, 3] == el
                neigh = face_neighbors[fc, 1]
            end
            if neigh ≠ 0
                temp[n+1] = neigh
                n += 1
            end
        end
        en = st + n
        g_edge_offset[el+1] = en
        g_edge_data[st:en-1] .= sort(temp[1:n])
        spectral_radius = max(spectral_radius, FT(en - st))
    end
    spectral_radius *= 2
    return Graph(g_nverts, g_edge_data, g_edge_offset, Array(1:g_nverts), spectral_radius)
end
"""
    laplacian_x_v!(vout::FTA1D, vin::FTA1D, grph::Graph)

This function computed the product of graph Laplacian matrix with vector `vin`.
"""
function laplacian_x_v!(
    vout::FTA1D,
    vin::FTA1D,
    grph::Graph,
) where {FT<:AbstractFloat,FTA1D<:Vector{FT}}
    @unpack nverts, edge_data, edge_offset = grph
    @assert length(vout) == length(vin) == nverts

    for i = 1:nverts
        vout[i] = FT(0)
        for j = edge_offset[i]:edge_offset[i+1]-1
            vout[i] += vin[i] - vin[edge_data[j]]
        end
    end
    return nothing
end
"""
    fiedler_vector(graph::Graph{<:Any,FT}, ::Lanczos)

This function estimates the Fiedler vector using Lanczos method.

Reference:
Estimating the Largest Eigenvalue by the Power and Lanczos Algorithms with a Random Start
Kuczynski, Jacek; Wozniakowski, Henryk

https://academiccommons.columbia.edu/doi/10.7916/D8B56SVV
"""
function fiedler_vector(graph::Graph{<:Any,FT}, ::Lanczos) where {FT} # using Lanczos
    nv = graph.nverts
    sr = graph.spectral_radius # convervative estimate of spectral radius of graph laplacian
    rel_tol = FT(1e-2) #FT(1e-3)
    maxit = min(
        Int(ceil((2.575 / rel_tol)^0.5 * log(nv))), # Kuczynski et. al.
        nv,
        1024,
    )

    ev = Vector{FT}(undef, nv)
    start = similar(ev)
    q0 = similar(ev)
    q1 = similar(ev)
    qn = similar(ev)
    H = zeros(FT, maxit, maxit)

    [start[i] = i^(-1.0) for i = 1:nv] # initial guess
    # orthogonalize wrt constant vector
    mean = sum(start) / nv
    start .-= mean
    start ./= norm(start) # normalize

    q1 .= start
    #-----------------------------------------
    for i = 1:maxit
        laplacian_x_v!(qn, q1, graph)
        for j = 1:nv
            qn[j] = sr * q1[j] - qn[j]
        end
        # orthogonalizing wrt cnst vector
        mean = sum(qn) / nv
        qn .-= mean

        #fac = sqrt(sum(q1.*qn))
        fac = sum(q1 .* qn)
        H[i, i] = fac

        for j = 1:nv
            qn[j] -= H[i, i] * q1[j]
        end

        if i > 1
            for j = 1:nv
                qn[j] = qn[j] - H[i-1, i] * q0[j]
            end
        end

        fac = sqrt(sum(qn .* qn))

        if i < maxit
            H[i+1, i] = fac
            H[i, i+1] = fac
        end

        qn ./= fac
        q0 .= q1
        q1 .= qn
    end
    il = iu = maxit
    vl = vu = FT(0) # not used
    abstol = 1.0E-14

    eigval, eigvec = LinearAlgebra.LAPACK.syevr!('V', 'I', 'L', H, vl, vu, il, iu, abstol)
    lam = sr - eigval[1]
    # computing the largest eigen vector
    ev .= FT(0)
    q1 .= start

    for i = 1:maxit
        for j = 1:nv
            ev[j] += q1[j] * eigvec[i]
        end
        laplacian_x_v!(qn, q1, graph)
        for j = 1:nv
            qn[j] = sr * q1[j] - qn[j]
        end
        # orthogonalizing wrt constant vector
        mean = sum(qn) / nv
        qn .-= mean

        fac = sum(q1 .* qn)
        H[i, i] = fac
        #-----------------------
        for j = 1:nv
            qn[j] -= H[i, i] * q1[j]
        end

        if i > 1
            for j = 1:nv
                qn[j] = qn[j] - H[i-1, i] * q0[j]
            end
        end

        fac = sum(qn .* qn)^0.5

        if i < maxit
            H[i+1, i] = fac
            H[i, i+1] = fac
        end

        qn ./= fac
        q0 .= q1
        q1 .= qn
        #-----------------------
    end
    ev ./= norm(ev)
    #-----------------------------------------
    return ev
end
"""
    fiedler_vector(graph::Graph{<:Any,FT}, ::PowerIt)

This function estimates the Fiedler vector using power iteration.
"""
function fiedler_vector(graph::Graph{<:Any,FT}, ::PowerIt) where {FT} # using power_iteration
    maxit = 100
    toler = 1.0E-14
    nv = graph.nverts
    sr = graph.spectral_radius # convervative estimate of spectral radius of graph laplacian
    # initial guess for fiedler vector
    split = cld(nv, 2)
    ev = Vector{FT}(undef, nv)
    [ev[i] = -1.0 for i = 1:split]
    [ev[i] = 1.0 for i = split+1:nv]

    # orthogonlizing wrt constant (first) eigen vector
    shift = sum(ev) / FT(nv)
    ev .-= shift
    ev ./= norm(ev) # normalize

    ev_old = deepcopy(ev)
    ev_temp = similar(ev)

    for i = 1:maxit-1
        # compute (I*sr - lmat) * ev
        laplacian_x_v!(ev_temp, ev, graph)
        for j = 1:nv
            ev[j] = ev[j] * sr - ev_temp[j]
        end
        ev ./= norm(ev) # normalize

        # orthogonlizing wrt constant (first) eigen vector
        shift = sum(ev) / FT(nv)
        ev .-= shift
        ev ./= norm(ev) # normalize

        diff = FT(1) - abs(sum(ev .* ev_old))
        if diff ≤ toler
            break
        else
            ev_old .= ev
        end
    end
    return ev
end
"""
    fiedler_vector(graph::Graph{<:Any,FT}, ::Direct)

This function computes the Fiedler vector using Julia/LAPACK subroutines.
This is intended for testing purposes only. Lanczos method is the
recommended default.
"""
function fiedler_vector(graph::Graph{<:Any,FT}, ::Direct) where {FT}
    lmat = laplacian_matrix_full(graph)
    vl, vu = FT(1), FT(2) # dummies, not used
    il, iu = 2, 2 # second eigen vector (Fiedler vector)
    #abstol = FT(1e-12)
    abstol = FT(1e-14)
    _, z = LinearAlgebra.LAPACK.syevr!('V', 'I', 'L', lmat, vl, vu, il, iu, abstol)
    return z[:]
end

"""
    laplacian_matrix_full(graph::Graph{<:Any,FT})

Computed the full Graph laplacian matrix explicitly. This is
recommended for testing purposes only.
"""
function laplacian_matrix_full(graph::Graph{<:Any,FT}) where {FT}
    @unpack nverts, edge_data, edge_offset = graph
    lmat = zeros(nverts, nverts)
    for i = 1:nverts
        st, en = edge_offset[i:i+1]
        for j = st:en-1
            lmat[edge_data[j], i] = -1
        end
        lmat[i, i] = en - st
    end
    return lmat
end
"""
    split_graph(graph::Graph{I,FT}, n1)

This function partitions the graph such that the first partition
has `n1` vertices and second partition has `graph.nverts - n1`
vertices.
"""
function split_graph(graph::Graph{I,FT}, n1) where {I,FT}
    @unpack nverts, vertices, edge_offset, edge_data = graph
    if n1 == nverts
        return graph, Graph(I, FT)
    end

    lvertno = deepcopy(vertices)
    #order = sortperm(fiedler_vector_direct(graph))
    #order = sortperm(fiedler_vector_power_it(graph))
    order = sortperm(fiedler_vector(graph, Lanczos()))

    n2 = nverts - n1
    sr1 = sr2 = FT(0)
    sort_ord1 = sort(order[1:n1])
    sort_ord2 = sort(order[n1+1:nverts])
    vertices1 = vertices[sort_ord1]
    vertices2 = vertices[sort_ord2]

    lvertno[sort_ord1] .= 1:n1    # local vertex numbers
    lvertno[sort_ord2] .= -(1:n2) # for split graphs

    edge_offset1 = zeros(I, n1 + 1)
    edge_offset2 = zeros(I, n2 + 1)

    ed1 = [zeros(I, i) for i in zeros(I, n1)]
    ed2 = [zeros(I, i) for i in zeros(I, n2)]

    for i = 1:nverts
        v1 = lvertno[i]
        for j = edge_offset[i]:edge_offset[i+1]-1
            v2 = lvertno[edge_data[j]]
            if v1 > 0 && v2 > 0
                edge_offset1[v1+1] += 1
                push!(ed1[v1], v2)
            elseif v1 < 0 && v2 < 0
                edge_offset2[-v1+1] += 1
                push!(ed2[-v1], -v2)
            end
        end
    end
    edge_offset1[1], edge_offset2[1] = 1, 1

    [edge_offset1[i+1] += edge_offset1[i] for i = 1:n1]
    [edge_offset2[i+1] += edge_offset2[i] for i = 1:n2]

    edge_data1 = Array{I}(undef, edge_offset1[end] - 1)
    edge_data2 = Array{I}(undef, edge_offset2[end] - 1)

    for i = 1:n1
        edge_data1[edge_offset1[i]:edge_offset1[i+1]-1] .= ed1[i]
        sr1 = max(sr1, length(ed1[i]))
    end
    sr1 *= 2
    for i = 1:n2
        edge_data2[edge_offset2[i]:edge_offset2[i+1]-1] .= ed2[i]
        sr2 = max(sr2, length(ed2[i]))
    end
    sr2 *= 2

    return Graph(n1, edge_data1, edge_offset1, vertices1, sr1),
    Graph(n2, edge_data2, edge_offset2, vertices2, sr2)
end

end

max_stars_repo_name: CliMA/GCMMesh

max_issues_repo_name: CliMA/GCMMesh

max_forks_repo_name: CliMA/GCMMesh


text:
# interfaces for defining new compilers

# the definition of a new GPU compiler is typically split in two:
# - a generic compiler that lives in GPUCompiler.jl (e.g., emitting PTX, SPIR-V, etc)
# - a more specific version in a package that targets an environment (e.g. CUDA, ROCm, etc)
#
# the first level of customizability is found in the AbstractCompilerTarget hierarchy,
# with methods and interfaces that can only be implemented within GPUCompiler.jl.
#
# further customization should be put in a concrete instance of the AbstractCompilerParams
# type, and can be used to customize interfaces defined on CompilerJob.

## target

export AbstractCompilerTarget

# container for state handled by targets defined in GPUCompiler.jl

abstract type AbstractCompilerTarget end

source_code(::AbstractCompilerTarget) = "text"

llvm_triple(::AbstractCompilerTarget) = error("Not implemented")

function llvm_machine(target::AbstractCompilerTarget)
    triple = llvm_triple(target)

    t = Target(triple=triple)

    tm = TargetMachine(t, triple)
    asm_verbosity!(tm, true)

    return tm
end

llvm_datalayout(target::AbstractCompilerTarget) = DataLayout(llvm_machine(target))

## params

export AbstractCompilerParams

# container for state handled by external users of GPUCompiler.jl

abstract type AbstractCompilerParams end

## function specification

export FunctionSpec

# what we'll be compiling

struct FunctionSpec{F,TT}
    f::Base.Callable
    tt::DataType
    kernel::Bool
    name::Union{Nothing,String}
end

# put the function and argument types in typevars
# so that we can access it from generated functions
FunctionSpec(f, tt=Tuple{}, kernel=true, name=nothing) =
    FunctionSpec{typeof(f),tt}(f, tt, kernel, name)

function signature(spec::FunctionSpec)
    fn = something(spec.name, nameof(spec.f))
    args = join(spec.tt.parameters, ", ")
    return "$fn($(join(spec.tt.parameters, ", ")))"
end

function Base.show(io::IO, spec::FunctionSpec)
    spec.kernel ? print(io, "kernel ") : print(io, "function ")
    print(io, signature(spec))
end

## job

export CompilerJob

# a specific invocation of the compiler, bundling everything needed to generate code

Base.@kwdef struct CompilerJob{T,P}
    target::T
    source::FunctionSpec
    params::P

    CompilerJob(target::AbstractCompilerTarget, source::FunctionSpec, params::AbstractCompilerParams) =
        new{typeof(target), typeof(params)}(target, source, params)
end

Base.similar(job::CompilerJob, source::FunctionSpec) =
    CompilerJob(target=job.target, source=source, params=job.params)

function Base.show(io::IO, job::CompilerJob{T}) where {T}
    print(io, "CompilerJob of ", job.source, " for ", T)
end

## interfaces and fallback definitions

# the Julia module to look up target-specific runtime functions in (this includes both
# target-specific functions from the GPU runtime library, like `malloc`, but also
# replacements functions for operations like `Base.sin`)
runtime_module(::CompilerJob) = error("Not implemented")

# check if a function is an intrinsic that can assumed to be always available
isintrinsic(::CompilerJob, fn::String) = false

# does this target support throwing Julia exceptions with jl_throw?
# if not, calls to throw will be replaced with calls to the GPU runtime
can_throw(::CompilerJob) = false

# generate a string that represents the type of compilation, for selecting a compiled
# instance of the runtime library. this slug should encode everything that affects
# the generated code of this compiler job (with exception of the function source)
runtime_slug(::CompilerJob) = error("Not implemented")

# early processing of the newly generated LLVM IR module
process_module!(::CompilerJob, mod::LLVM.Module) = return

# early processing of the newly identified LLVM kernel function
process_kernel!(::CompilerJob, mod::LLVM.Module, kernel::LLVM.Function) = return kernel

# final processing of the IR module, right before validation and machine-code generation
finish_module!(::CompilerJob, mod::LLVM.Module) = return

add_lowering_passes!(::CompilerJob, pm::LLVM.PassManager) = return

add_optimization_passes!(::CompilerJob, pm::LLVM.PassManager) = return

link_libraries!(::CompilerJob, mod::LLVM.Module, undefined_fns::Vector{String}) = return


text:
#import MultivariateStats
import Combinatorics

"""
    sMSC(datacube, cycle_length)

subtract the median seasonal cycle from the datacube given the length of year `cycle_length`.

# Examples
```
julia> dc = hcat(rand(193) + 2* sin.(0:pi/24:8*pi), rand(193) + 2* sin.(0:pi/24:8*pi))
julia> sMSC_dc = sMSC(dc, 48)
```
"""
function sMSC(datacube, cycle_length)
    x_out = similar(datacube)
    removeMSC!(datacube, x_out, cycle_length, 1)
    return(x_out)
end

function removeMSC!(xin::AbstractArray{Tp,ndim},xout::AbstractArray{Tp,ndim},NpY::Integer,itimedim::Integer;imscstart::Int=1) where {Tp,ndim}
   # Determine length of reshape dimensions
   s=size(xin)
   l1=itimedim==1 ? 1 : prod(s[1:(itimedim-1)])
   l2=itimedim==ndim ? 1 : prod(s[(itimedim+1):end])
   ltime=s[itimedim]
   #Reshape the cube to squeeze unimportant variables
   xin2=reshape(xin,l1,ltime,l2);
   xout2=reshape(xout,l1,ltime,l2);
   msc=zeros(Float64,NpY); # This is the array where the temp msc is stored
   nmsc=zeros(Int,NpY); # This is for counting how many values were added
   #Start loop through all other variables
   for i2=1:l2,i1=1:l1
       imsc=imscstart
       fill!(msc,zero(Float64))
       fill!(nmsc,zero(Int))
       for itime=1:ltime
           curval=xin2[i1,itime,i2]
           if !isnan(curval)
               msc[imsc]  += curval
               nmsc[imsc] += 1
           end
           imsc =imsc==NpY ? 1 : imsc+1 # Increase msc time step counter
       end
       for i in 1:NpY msc[i] = nmsc[i] > 0 ? msc[i]/nmsc[i] : NaN end # Get MSC by dividing by number of points
       imsc=imscstart
       for i in 1:ltime
           xout2[i1,i,i2] = xin2[i1,i,i2]-msc[imsc]
           imsc =imsc==NpY ? 1 : imsc+1 # Increase msc time step counter
       end
   end
   #Copy data if necessary
   pointer(xout2) != pointer(xout) && copyto!(xout,xout2)
  return(xout)
end

#small helper
function func(i,x)
        return(x[i])
end

"""
    TDE(datacube::Array{tp, 4}, ΔT::Integer, DIM::Int = 3) where {tp}
    TDE(datacube::Array{tp, 3}, ΔT::Integer, DIM::Int = 3) where {tp}

returns an embedded datacube by concatenating lagged versions of the 2-, 3- or 4-dimensional datacube with `ΔT` time steps in the past up to dimension `DIM` (presetting: `DIM = 3`)

# Examples
```
julia> dc = randn(50,3)
julia> TDE(dc, 3, 2)
```
"""
function TDE(datacube::AbstractArray{tp, 4}, ΔT::Integer, DIM::Int = 3) where {tp}
    start = ((DIM-1)*ΔT+1)
    embedded_datacube = zeros(tp, (size(datacube, 1) - start +1, size(datacube, 2), size(datacube, 3), size(datacube, 4) * DIM))
    for dim = 1:DIM
    embedded_datacube[:,:,:,((dim-1)*size(datacube, 4)+1):(dim*size(datacube, 4))] =
        datacube[(start - ((dim-1) * ΔT)) : (size(datacube, 1) - ((dim-1) * ΔT)),:,:,:]
    end
    return(embedded_datacube)
end

function TDE(datacube::AbstractArray{tp, 3}, ΔT::Integer, DIM::Int = 3) where {tp}
    start = ((DIM-1)*ΔT+1)
    embedded_datacube = zeros(tp, (size(datacube, 1) - start +1, size(datacube, 2), size(datacube, 3) * DIM))
    for dim = 1:DIM
    embedded_datacube[:,:,((dim-1)*size(datacube, 3)+1):(dim*size(datacube, 3))] =
        datacube[(start - ((dim-1) * ΔT)) : (size(datacube, 1) - ((dim-1) * ΔT)),:,:]
    end
    return(embedded_datacube)
end

function TDE(datacube::AbstractArray{tp, 2}, ΔT::Integer, DIM::Int = 3) where {tp}
    start = ((DIM-1)*ΔT+1)
    embedded_datacube = zeros(tp, (size(datacube, 1) - start +1,  size(datacube, 2) * DIM))
    for dim = 1:DIM
    embedded_datacube[:,((dim-1)*size(datacube, 2)+1):(dim*size(datacube, 2))] =
        datacube[(start - ((dim-1) * ΔT)) : (size(datacube, 1) - ((dim-1) * ΔT)),:]
    end
    return(embedded_datacube)
end

function TDE(datacube::AbstractArray{tp, 1}, ΔT::Integer, DIM::Int = 3) where {tp}
    start = ((DIM-1)*ΔT+1)
    embedded_datacube = zeros(tp, (size(datacube, 1) - start +1,   DIM))
    for dim = 1:DIM
    embedded_datacube[:,((dim-1)+1):(dim)] =
        datacube[(start - ((dim-1) * ΔT)) : (size(datacube, 1) - ((dim-1) * ΔT))]
    end
    return(embedded_datacube)
end

"""
    mw_VAR(datacube::Array{tp,N}, windowsize::Int = 10) where {tp,N}

compute the variance in a moving window along the first dimension of the datacube (presetting: `windowsize = 10`).
Accepts N dimensional datacubes.

# Examples
```
julia> dc = randn(50,3,3,3)
julia> mw_VAR(dc, 15)
```
"""
function mw_VAR(datacube::AbstractArray{<:AbstractFloat,N}, windowsize::Int = 10) where {N}
  out = zeros(eltype(datacube), size(datacube))
  for k = 1:length(out) out[k] = NaN end
  datacube0mean = datacube .- mean(datacube, dims = 1)
  mw_VAR!(out, datacube0mean, windowsize)
  return(out)
end

"""
    mw_VAR!(out::Array{tp, N}, datacube0mean::Array{tp,N}, windowsize::Int = 10) where {tp,N}

mutating version for `mw_VAR()`. The mean of the input data `datacube0mean` has to be 0. Initialize out properly: `out = datacube0mean` leads to wrong results.

"""
function mw_VAR!(out::AbstractArray{tp, N}, datacube0mean::AbstractArray{tp,N}, windowsize::Int = 10) where {tp, N}
  T = size(datacube0mean, 1)
  for beg = 1:T:length(datacube0mean)
    inner_mw_VAR!(out, datacube0mean, windowsize, beg, T)
  end
   broadcast!(/, out, out, windowsize)
  return(out)
end

"""
    inner_mw_VAR!(out::Array{tp, N}, datacube0mean::Array{tp,N}, windowsize::Int = 10) where {tp,N}

internal function for mw_VAR!()

"""
function inner_mw_VAR!(out, datacube, windowsize, beg, T) # mean already subtracted
  out_beg = floor(Int, windowsize * 0.5) + beg
  for notavailable =  (beg):(out_beg-1)
    out[notavailable] = NaN
  end
  # init x
  x = sum(unsafe_wrap(Array, pointer(datacube, beg), windowsize).^2)
  out[out_beg] = x
  for i = 1:(T-windowsize)
    x = x - datacube[beg + i - 1]^2 + datacube[beg + i - 1 + windowsize]^2
    out[out_beg + i] = x
  end
  for notavailable = (out_beg + (T-windowsize +1) + 1):(beg+T-1)
    out[notavailable] = NaN
  end
  return(out)
end

"""
    mw_AVG(datacube::AbstractArray{tp,N}, windowsize::Int = 10) where {tp,N}

compute the average in a moving window along the first dimension of the datacube (presetting: `windowsize = 10`).
Accepts N dimensional datacubes.

# Examples
```
julia> dc = randn(50,3,3,3)
julia> mw_AVG(dc, 15)
```
"""
function mw_AVG(datacube::AbstractArray{<:AbstractFloat,N}, windowsize::Int = 10) where {N}
  out = zeros(eltype(datacube), size(datacube))
  mw_AVG!(out, datacube, windowsize)
  return(out)
end

"""
    mw_AVG!(out::Array{tp, N}, datacube::Array{tp,N}, windowsize::Int = 10) where {tp,N}

internal and mutating version for `mw_AVG()`.

"""
function mw_AVG!(out::AbstractArray{tp,N}, datacube::AbstractArray{tp,N}, windowsize::Int = 10) where {tp,N}
  for k = 1:length(out) out[k] = NaN end
  T = size(datacube, 1)
  for beg = 1:T:length(datacube)
    inner_mw_AVG!(out, datacube, windowsize, beg, T)
  end
   broadcast!(/, out, out, windowsize)
  return(out)
end

function inner_mw_AVG!(out::AbstractArray{tp, N}, datacube::AbstractArray{tp, N}, windowsize::Int, beg::Int, T::Int) where {tp,N}
  out_beg = floor(Int, windowsize * 0.5) + beg
  for notavailable =  (beg):(out_beg-1)
    out[notavailable] = NaN
  end
  # init x
  x = sum(unsafe_wrap(Array, pointer(datacube, beg), windowsize))
  out[out_beg] = x
  for i = 1:(T - windowsize)
    x = x - datacube[beg + i - 1] + datacube[beg + i - 1 + windowsize]
    out[out_beg + i] = x
  end
  for notavailable = (out_beg + (T-windowsize +1) + 1):(beg+T-1)
    out[notavailable] = NaN
  end
  return(out)
end

"""
    mw_COR(datacube::Array{tp, 4}, windowsize::Int = 10) where {tp}

compute the correlation in a moving window along the first dimension of the datacube (presetting: `windowsize = 10`).
Accepts 4-dimensional datacubes.

"""
function mw_COR(datacube::AbstractArray{tp, 4}, windowsize::Int = 10) where {tp}
  comb = collect(Combinatorics.combinations(1:size(datacube, 4), 2))
  out = zeros(Float64, size(datacube, 1), size(datacube, 2), size(datacube, 3), size(comb, 1))
  x = zeros(Float64, windowsize, 2)
  for icomb = 1:size(comb, 1)
    for lon = 1:size(datacube, 3)
      for lat = 1:size(datacube, 2)
        for t = 1:(size(datacube, 1)-windowsize)
          x = copyto!(x, datacube[t:(t+windowsize-1),lat,lon,comb[icomb]])
          copyto!(view(out, t + round.(Int, windowsize * 0.5),lat,lon,icomb)
                , cor(view(x,:,1), view(x,:,2)))
        end
      end
    end
  end
  return(comb, out)
end

function innerst_ewma!(z::AbstractArray{tp, N}, dat::AbstractArray{tp, N}, beg::Int, len::Int, lambda::Float64 = 0.15) where {tp, N}
    z[beg] = dat[beg]*lambda
    one_m_lambda = 1.0-lambda
    for i = (beg+1):(beg+len-1)
            z[i] = dat[i] * lambda + z[i-1] * one_m_lambda
    end
  return(z)
end

"""
    EWMA(dat,  λ)

Compute the exponential weighted moving average (EWMA) with the weighting parameter `λ` between 0 (full weighting) and 1 (no weighting) along the first dimension of `dat`.
Supports N-dimensional Arrays.

Lowry, C. A., & Woodall, W. H. (1992). A Multivariate Exponentially Weighted Moving Average Control Chart. Technometrics, 34, 46–53.

# Examples
```
julia> dc = rand(100,3,2)
julia> ewma_dc = EWMA(dc, 0.1)
```
"""
function EWMA(dat::AbstractArray{tp, N}, lambda::Float64 = 0.15) where {tp, N}
  Z = similar(dat)
  EWMA!(Z, dat, lambda)
  return(Z)
end

"""
    EWMA!(Z, dat,  λ)

use a preallocated output Z. `Z = similar(dat)` or `dat = dat` for overwriting itself.

# Examples
```
julia> dc = rand(100,3,2)
julia> EWMA!(dc, dc, 0.1)
```
"""
function EWMA!(Z::AbstractArray{tp, N}, dat::AbstractArray{tp, N}, lambda::Float64 = 0.15) where {tp, N}
  T = size(dat, 1)
  LENGTH = length(dat)
  for istart = 1:T:LENGTH
    innerst_ewma!(Z, dat, istart, T, lambda)
  end
  return(Z)
end

"""
    init_MedianCycle(dat::Array{tp}, cycle_length::Int = 46)
    init_MedianCycle(temporal_length::Int[, cycle_length::Int = 46])

initialises an init_MC object to be used as input for `get_MedianCycle!()`. Input is either some sample data or the temporal lenght of the expected input vector
and the length of the annual cycle (presetting: `cycle_length = 46`)
"""
function init_MedianCycle(dat::AbstractArray{tp}, cycle_length::Int = 46) where {tp}
  complete_years = Int(floor(size(dat, 1)/cycle_length))
  cycle_dat = zeros(tp, cycle_length, complete_years)
  cycle_medians = zeros(tp, cycle_length)
  init_MC = (cycle_length, complete_years, cycle_medians, cycle_dat)
  return(init_MC)
end

function init_MedianCycle(temporal_length::Int, cycle_length::Int = 46)
  complete_years = Int(floor(temporal_length/cycle_length))
  cycle_dat = zeros(tp, cycle_length, complete_years)
  cycle_medians = zeros(tp, cycle_length)
  init_MC = (cycle_length, complete_years, cycle_medians, cycle_dat)
  return(init_MC)
end

"""
    get_MedianCycle!(init_MC, dat::Array{tp,1})

Memory efficient version of `get_MedianCycle()`, returning the median cycle in `init_MC[3]`. The `init_MC` object should be created with `init_MedianCycle`.
Can deal with some NaN values.

# Examples
```
julia> using MultivariateAnomalies
julia> dat = rand(193) + 2* sin(0:pi/24:8*pi)
julia> dat[100] = NaN
julia> init_MC = init_MedianCycle(dat, 48)
julia> get_MedianCycle!(init_MC, dat)
julia> init_MC[3]
```
"""
function get_MedianCycle!(init_MC::Tuple{Int64,Int64,Array{tp,1},Array{tp,2}}, dat::Array{tp,1}) where {tp}
  copyto!(init_MC[4], view(dat, 1:(init_MC[2] * init_MC[1])))
  for i = 1:init_MC[1]
    if(all(.!isnan.(view(init_MC[4], i, :))))
      copyto!(view(init_MC[3], i), median(view(init_MC[4], i, :)))
    elseif(all(isnan.(view(init_MC[4], i, :))))
        copyto!(view(init_MC[3], i), NaN)
    else
      copyto!(view(init_MC[3], i), median(view(init_MC[4], i, squeeze(!isnan.(view(init_MC[4], i,:)), 1))))
    end
  end
  return(init_MC[3])
end

"""
    get_MedianCycle(dat::Array{tp,1}, cycle_length::Int = 46)

returns the median annual cycle of a one dimensional data array, given the length of the annual cycle (presetting: `cycle_length = 46`).
Can deal with some NaN values.

# Examples
```
julia> using MultivariateAnomalies
julia> dat = randn(90) + x = sind.(0:8:719)
julia> cycles = get_MedianCycle(dat, 48)
```
"""
function get_MedianCycle(dat::AbstractArray{tp,1}, cycle_length::Int = 46) where {tp}
  init_MC = init_MedianCycle(dat, cycle_length)
  get_MedianCycle!(init_MC, dat)
  return(init_MC[3])
end

"""
    get_MedianCycles(datacube, cycle_length::Int = 46)

returns the median annual cycle of a datacube, given the length of the annual cycle (presetting: `cycle_length = 46`).
The datacube can be 2, 3, 4-dimensional, time is stored along the first dimension.

# Examples
```
julia> using MultivariateAnomalies
julia> dc = hcat(rand(193) + 2* sin(0:pi/24:8*pi), rand(193) + 2* sin(0:pi/24:8*pi))
julia> cycles = get_MedianCycles(dc, 48)
```
"""
function get_MedianCycles(datacube::AbstractArray{tp,4}, cycle_length::Int = 46) where {tp}
dat = zeros(Float64, size(datacube, 1));
init_MC = init_MedianCycle(dat, cycle_length);
med_cycles_out = zeros(Float64, cycle_length, size(datacube, 2), size(datacube, 3), size(datacube, 4));
# loop
  for var = 1:size(datacube, 4)
    for lon = 1:size(datacube, 3)
      for lat = 1:size(datacube, 2)
        copyto!(dat, view(datacube, :,lat,lon,var))
        copyto!(view(med_cycles_out, :, lat, lon, var), get_MedianCycle!(init_MC, dat))
      end
    end
  end
  return(med_cycles_out)
end

function get_MedianCycles(datacube::AbstractArray{tp,3}, cycle_length::Int = 46) where {tp}
dat = zeros(tp, size(datacube, 1));
init_MC = init_MedianCycle(dat, cycle_length);
med_cycles_out = zeros(tp, cycle_length, size(datacube, 2), size(datacube, 3));
# loop
  for var = 1:size(datacube, 3)
      for lat = 1:size(datacube, 2)
        copyto!(dat, view(datacube, :, lat,var))
        copyto!(view(med_cycles_out, :, lat, var), get_MedianCycle!(init_MC, dat))
      end
  end
  return(med_cycles_out)
end

function get_MedianCycles(datacube::AbstractArray{tp,2}, cycle_length::Int = 46) where {tp}
dat = zeros(tp, size(datacube, 1));
init_MC = init_MedianCycle(dat, cycle_length);
med_cycles_out = zeros(tp, cycle_length, size(datacube, 2));
# loop
  for var = 1:size(datacube, 2)
        copyto!(dat, view(datacube, :, var))
        copyto!(view(med_cycles_out, :, var), get_MedianCycle!(init_MC, dat))
  end
  return(med_cycles_out)
end

# get moving window loop with MWobj
mutable struct MWobj
    ObsPerYear::Int
    windowsize::Int
    edgecut::Int
    startidx::Int
    numoutvars::Int
    minwindow::Int
    maxwindow::Int
    count::Array{Int, 1}
    iterator_windowcenter::StepRange{Int,Int}
    mw_indat
    mw_outdat
    mw_idx::Array{Int, 1}
    xin
    xout
end

function init_MovingWindow(xin::AbstractArray{tp, 2}; ObsPerYear::Int = 46, windowsize::Int = 11, edgecut::Int = 0, startidx::Int = 1, numoutvars::Int = 0) where {tp}
    if 2*edgecut >= windowsize error("2 * edgecut has to be smaller windowsize, but is $edgecut, windowsize = $windowsize") end
    if !isodd(windowsize) error("windowsize has to be odd, but is $windowsize") end
    if round(Int, size(xin, 1) / ObsPerYear) * ObsPerYear != size(xin, 1) error("ObsPerYear multiplied by some integer is not matching size(xin, 1)") end
    if numoutvars == 0 numoutvars = size(xin, 2) end
    if numoutvars > 1
        mwobj = MWobj(
        ObsPerYear
        ,windowsize
        ,edgecut
        ,startidx
        ,numoutvars
        ,-floor(Int, windowsize / 2.0) # minwindow
        ,floor(Int, windowsize / 2.0) # maxwindow
        , [0] # count
        ,startidx:(windowsize-2*edgecut):ObsPerYear #iterator_windowcenter
        ,zeros(eltype(xin), windowsize * floor(Int, size(xin, 1) / ObsPerYear), size(xin, 2)) #mw_indat
        ,zeros(eltype(xin), windowsize * floor(Int, size(xin, 1) / ObsPerYear), numoutvars) #mw_outdat
        ,zeros(Int, windowsize * floor(Int, size(xin, 1) / ObsPerYear)) # mw_idx
        ,xin
        ,zeros(eltype(xin), size(xin, 1), numoutvars) #xout
        )
    else
        mwobj = MWobj(
        ObsPerYear
        ,windowsize
        ,edgecut
        ,startidx
        ,numoutvars
        ,-floor(Int, windowsize / 2.0) # minwindow
        ,floor(Int, windowsize / 2.0) # maxwindow
        , [0] # count
        ,startidx:(windowsize-2*edgecut):ObsPerYear #iterator_windowcenter
        ,zeros(eltype(xin), windowsize * floor(Int, size(xin, 1) / ObsPerYear), size(xin, 2)) #mw_indat
        ,zeros(eltype(xin), windowsize * floor(Int, size(xin, 1) / ObsPerYear)) #mw_outdat
        ,zeros(Int, windowsize * floor(Int, size(xin, 1) / ObsPerYear)) # mw_idx
        ,xin
        ,zeros(eltype(xin), size(xin, 1), numoutvars) #xout
        )
    end
    return mwobj
end

function getMWData!(mwobj::MWobj, windowcenter::Int)
    xin = mwobj.xin
    mwdata = mwobj.mw_indat
    mwidx = mwobj.mw_idx
    ObsPerYear = mwobj.ObsPerYear
    windowsize = mwobj.windowsize
    minwindow =  mwobj.minwindow
    maxwindow = mwobj.maxwindow
    count = mwobj.count

    count[1] = 0
    for i = windowcenter:ObsPerYear:(size(xin, 1)+windowsize)
      for j = minwindow:maxwindow
        tempidx = i+j
        if tempidx > 0 && tempidx <= size(xin, 1)
          count[1] += 1
          mwidx[count[1]] = tempidx
          for varidx = 1:size(xin, 2)
            mwdata[count[1], varidx] = xin[tempidx, varidx]
          end
        end
      end
    end

  return view(mwobj.mw_indat,1:mwobj.count[1],:)
end

function pushMWResultsBack!(mwobj::MWobj, windowcenter::Int)
    mwdata = mwobj.mw_outdat
    mwidx = mwobj.mw_idx
    xout = mwobj.xout
    ObsPerYear=mwobj.ObsPerYear
    windowsize = mwobj.windowsize
    edgecut = mwobj.edgecut
    minwindow = mwobj.minwindow
    maxwindow = mwobj.maxwindow
    count = mwobj.count

  # assumes Time * Vars
    count[1] = 0
    for i = windowcenter:ObsPerYear:(size(xout, 1)+windowsize)
      for j = minwindow:maxwindow
        tempidx = i+j
        if tempidx > 0 && tempidx <= size(xout, 1)
          count[1] += 1
          if mwidx[count[1]] != tempidx error("not matching mwidx[count] = $(mwidx[count]), check windowsize") end
          if j >= (minwindow + edgecut) && j <= (maxwindow - edgecut)
              for varidx = 1:size(xout, 2)
                xout[tempidx, varidx] = mwdata[count[1], varidx]
              end
          end
        end
      end
    end

    return view(mwobj.mw_outdat, 1:mwobj.count[1], :)
end

"""
    mapMovingWindow(function2mw, x, args...; ObsPerYear::Int = 46, windowsize::Int = 9, edgecut::Int = 0, startidx::Int = 1, numoutvars::Int = 0)
    mapMovingWindow(function2mw, x; ...)
apply a function (function2mw) in a moving window encompassing all years and running along the time, e.g. apply the function to all summers, then summers + 1 timestep ...
results are written to the center of the respective windowsize. Input axes are time or time-variables. The number of output variables can be different from the number of input variables, specified in numoutvars.
e.g. transforming the variables in normalised ranks between zero and one to get rid of heteroscedasticity would look like:
    using MultivariateAnomalies
    x = randn(10*46, 3)
    mapMovingWindow(get_quantile_scores, x, numoutvars = size(x, 2))
"""
function mapMovingWindow(function2mw, x, args...; ObsPerYear::Int = 46, windowsize::Int = 9, edgecut::Int = 0, startidx::Int = 1, numoutvars::Int = 0)
  mwobj = init_MovingWindow(x, ObsPerYear = ObsPerYear, windowsize = windowsize, edgecut = edgecut, numoutvars = numoutvars, startidx = startidx)
  if numoutvars > 1
      for windowcenter = mwobj.iterator_windowcenter
        getMWData!(mwobj, windowcenter)
        # do something with mwobj.mw_indat and write the results to mwobj.mw_outdat
        xout = view(mwobj.mw_outdat, 1:mwobj.count[1], :)
        xin = view(mwobj.mw_indat, 1:mwobj.count[1], :)
        xout[:] = function2mw(xin, args...)
        pushMWResultsBack!(mwobj, windowcenter)
    end
  else
      for windowcenter = mwobj.iterator_windowcenter
        getMWData!(mwobj, windowcenter)
        # do something with mwobj.mw_indat and write the results to mwobj.mw_outdat
        xout = view(mwobj.mw_outdat, 1:mwobj.count[1])
        xin = view(mwobj.mw_indat, 1:mwobj.count[1], :)
        xout[:] = function2mw(xin, args...)
        pushMWResultsBack!(mwobj, windowcenter)
    end
  end
  return mwobj.xout
end

function mapMovingWindow(function2mw, x; ObsPerYear::Int = 46, windowsize::Int = 9, edgecut::Int = 0, startidx::Int = 1, numoutvars::Int = 0)
  mwobj = init_MovingWindow(x, ObsPerYear = ObsPerYear, windowsize = windowsize, edgecut = edgecut, numoutvars = numoutvars, startidx = startidx)
  if numoutvars > 1
      for windowcenter = mwobj.iterator_windowcenter
        getMWData!(mwobj, windowcenter)
        # do something with mwobj.mw_indat and write the results to mwobj.mw_outdat
        xout = view(mwobj.mw_outdat, 1:mwobj.count[1], :)
        xin = view(mwobj.mw_indat, 1:mwobj.count[1], :)
        xout[:] = function2mw(xin)
        pushMWResultsBack!(mwobj, windowcenter)
    end
  else
      for windowcenter = mwobj.iterator_windowcenter
        getMWData!(mwobj, windowcenter)
        # do something with mwobj.mw_indat and write the results to mwobj.mw_outdat
        xout = view(mwobj.mw_outdat, 1:mwobj.count[1])
        xin = view(mwobj.mw_indat, 1:mwobj.count[1], :)
        xout[:] = function2mw(xin)
        pushMWResultsBack!(mwobj, windowcenter)
    end
  end
  return mwobj.xout
end

###################################
#end