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

type BatchProcessingError <: Exception
    data
    ex
end

"""
    pgenerate([::WorkerPool], f, c...) -> iterator

Apply `f` to each element of `c` in parallel using available workers and tasks.

For multiple collection arguments, apply f elementwise.

Results are returned in order as they become available.

Note that `f` must be made available to all worker processes; see
[Code Availability and Loading Packages](:ref:`Code Availability
and Loading Packages <man-parallel-computing-code-availability>`)
for details.
"""
function pgenerate(p::WorkerPool, f, c)
    if length(p) == 0
        return AsyncGenerator(f, c; ntasks=()->nworkers(p))
    end
    batches = batchsplit(c, min_batch_count = length(p) * 3)
    return Iterators.flatten(AsyncGenerator(remote(p, b -> asyncmap(f, b)), batches))
end
pgenerate(p::WorkerPool, f, c1, c...) = pgenerate(p, a->f(a...), zip(c1, c...))
pgenerate(f, c) = pgenerate(default_worker_pool(), f, c)
pgenerate(f, c1, c...) = pgenerate(a->f(a...), zip(c1, c...))

"""
    pmap([::AbstractWorkerPool], f, c...; distributed=true, batch_size=1, on_error=nothing, retry_n=0, retry_max_delay=DEFAULT_RETRY_MAX_DELAY, retry_on=DEFAULT_RETRY_ON) -> collection

Transform collection `c` by applying `f` to each element using available
workers and tasks.

For multiple collection arguments, apply f elementwise.

Note that `f` must be made available to all worker processes; see
[Code Availability and Loading Packages](:ref:`Code Availability
and Loading Packages <man-parallel-computing-code-availability>`)
for details.

If a worker pool is not specified, all available workers, i.e., the default worker pool
is used.

By default, `pmap` distributes the computation over all specified workers. To use only the
local process and distribute over tasks, specify `distributed=false`. This is equivalent to `asyncmap`.

`pmap` can also use a mix of processes and tasks via the `batch_size` argument. For batch sizes
greater than 1, the collection is split into multiple batches, which are distributed across
workers. Each such batch is processed in parallel via tasks in each worker. The specified
`batch_size` is an upper limit, the actual size of batches may be smaller and is calculated
depending on the number of workers available and length of the collection.

Any error stops pmap from processing the remainder of the collection. To override this behavior
you can specify an error handling function via argument `on_error` which takes in a single argument, i.e.,
the exception. The function can stop the processing by rethrowing the error, or, to continue, return any value
which is then returned inline with the results to the caller.

Failed computation can also be retried via `retry_on`, `retry_n`, `retry_max_delay`, which are passed through
to `retry` as arguments `retry_on`, `n` and `max_delay` respectively. If batching is specified, and an entire batch fails,
all items in the batch are retried.

The following are equivalent:

* `pmap(f, c; distributed=false)` and `asyncmap(f,c)`
* `pmap(f, c; retry_n=1)` and `asyncmap(retry(remote(f)),c)`
* `pmap(f, c; retry_n=1, on_error=e->e)` and `asyncmap(x->try retry(remote(f))(x) catch e; e end, c)`
"""
function pmap(p::AbstractWorkerPool, f, c; distributed=true, batch_size=1, on_error=nothing,
                                           retry_n=0,
                                           retry_max_delay=DEFAULT_RETRY_MAX_DELAY,
                                           retry_on=DEFAULT_RETRY_ON,
                                           # deprecated keyword args:
                                           err_retry=nothing, err_stop=nothing, pids=nothing)
    #15409
    if err_retry !== nothing
        depwarn("err_retry is deprecated, use pmap(retry(f), c...).", :pmap)
        if err_retry == true
            f = retry(f)
        end
    end
    if pids !== nothing
        depwarn("pids is deprecated, use pmap(::WorkerPool, f, c...).", :pmap)
        p = WorkerPool(pids)
    end
    if err_stop !== nothing
        depwarn("err_stop is deprecated, use pmap(f, c...; on_error = error_handling_func).", :pmap)
        if err_stop === false
            on_error = e->e
        end
    end

    f_orig = f
    # Don't do remote calls if there are no workers.
    if (length(p) == 0) || (length(p) == 1 && fetch(p.channel) == myid())
        distributed = false
    end

    # Don't do batching if not doing remote calls.
    if !distributed
        batch_size = 1
    end

    # If not batching, do simple remote call.
    if batch_size == 1
        if distributed
            f = remote(p, f)
        end

        if retry_n > 0
            f = wrap_retry(f, retry_on, retry_n, retry_max_delay)
        end
        if on_error !== nothing
            f = wrap_on_error(f, on_error)
        end

        return collect(AsyncGenerator(f, c; ntasks=()->nworkers(p)))
    else
        batches = batchsplit(c, min_batch_count = length(p) * 3,
                                max_batch_size = batch_size)

        # During batch processing, We need to ensure that if on_error is set, it is called
        # for each element in error, and that we return as many elements as the original list.
        # retry, if set, has to be called element wise and we will do a best-effort
        # to ensure that we do not call mapped function on the same element more than retry_n.
        # This guarantee is not possible in case of worker death / network errors, wherein
        # we will retry the entire batch on a new worker.
        if (on_error !== nothing) || (retry_n > 0)
            f = wrap_on_error(f, (x,e)->BatchProcessingError(x,e); capture_data=true)
        end
        f = wrap_batch(f, p, on_error)
        results = collect(Iterators.flatten(AsyncGenerator(f, batches; ntasks=()->nworkers(p))))
        if (on_error !== nothing) || (retry_n > 0)
            process_batch_errors!(p, f_orig, results, on_error, retry_on, retry_n, retry_max_delay)
        end

        return results
    end
end

pmap(p::AbstractWorkerPool, f, c1, c...; kwargs...) = pmap(p, a->f(a...), zip(c1, c...); kwargs...)
pmap(f, c; kwargs...) = pmap(default_worker_pool(), f, c; kwargs...)
pmap(f, c1, c...; kwargs...) = pmap(a->f(a...), zip(c1, c...); kwargs...)

function wrap_on_error(f, on_error; capture_data=false)
    return x -> begin
        try
            f(x)
        catch e
            if capture_data
                on_error(x, e)
            else
                on_error(e)
            end
        end
    end
end

wrap_retry(f, retry_on, n, max_delay) = retry(f, retry_on; n=n, max_delay=max_delay)

function wrap_batch(f, p, on_error)
    f = asyncmap_batch(f)
    return batch -> begin
        try
            remotecall_fetch(f, p, batch)
        catch e
            if on_error !== nothing
                return Any[BatchProcessingError(batch[i], e) for i in 1:length(batch)]
            else
                rethrow(e)
            end
        end
    end
end

asyncmap_batch(f) = batch -> asyncmap(f, batch)

function process_batch_errors!(p, f, results, on_error, retry_on, retry_n, retry_max_delay)
    # Handle all the ones in error in another pmap, with batch size set to 1
    if (on_error !== nothing) || (retry_n > 0)
        reprocess = []
        for (idx, v) in enumerate(results)
            if isa(v, BatchProcessingError)
                push!(reprocess, (idx,v))
            end
        end

        if length(reprocess) > 0
            errors = [x[2] for x in reprocess]
            exceptions = [x.ex for x in errors]
            if (retry_n > 0) && all([retry_on(ex) for ex in exceptions])
                retry_n = retry_n - 1
                error_processed = pmap(p, f, [x.data for x in errors];
                                                    on_error=on_error,
                                                    retry_on=retry_on,
                                                    retry_n=retry_n,
                                                    retry_max_delay=retry_max_delay)
            elseif on_error !== nothing
                error_processed = map(on_error, exceptions)
            else
                throw(CompositeException(exceptions))
            end

            for (idx, v) in enumerate(error_processed)
                results[reprocess[idx][1]] = v
            end
        end
    end
    nothing
end

"""
    head_and_tail(c, n) -> head, tail

Returns `head`: the first `n` elements of `c`;
and `tail`: an iterator over the remaining elements.

```jldoctest
julia> a = 1:10
1:10

julia> b, c = Base.head_and_tail(a, 3)
([1,2,3],Base.Iterators.Rest{UnitRange{Int64},Int64}(1:10,4))

julia> collect(c)
7-element Array{Any,1}:
  4
  5
  6
  7
  8
  9
 10
```
"""
function head_and_tail(c, n)
    head = Vector{eltype(c)}(n)
    s = start(c)
    i = 0
    while i < n && !done(c, s)
        i += 1
        head[i], s = next(c, s)
    end
    return resize!(head, i), Iterators.rest(c, s)
end

"""
    batchsplit(c; min_batch_count=1, max_batch_size=100) -> iterator

Split a collection into at least `min_batch_count` batches.

Equivalent to `partition(c, max_batch_size)` when `length(c) >> max_batch_size`.
"""
function batchsplit(c; min_batch_count=1, max_batch_size=100)
    if min_batch_count < 1
        throw(ArgumentError("min_batch_count must be ≥ 1, got $min_batch_count"))
    end

    if max_batch_size < 1
        throw(ArgumentError("max_batch_size must be ≥ 1, got $max_batch_size"))
    end

    # Split collection into batches, then peek at the first few batches
    batches = Iterators.partition(c, max_batch_size)
    head, tail = head_and_tail(batches, min_batch_count)

    # If there are not enough batches, use a smaller batch size
    if length(head) < min_batch_count
        batch_size = max(1, div(sum(length, head), min_batch_count))
        return Iterators.partition(collect(Iterators.flatten(head)), batch_size)
    end

    return Iterators.flatten((head, tail))
end

max_stars_repo_name: habemus-papadum/julia

max_issues_repo_name: habemus-papadum/julia

max_forks_repo_name: habemus-papadum/julia


text:
using StochasticDiffEq, Zygote
using DiffEqSensitivity, Test, ForwardDiff

abstol = 1e-12
reltol = 1e-12
savingtimes = 0.5

function test_SDE_callbacks()
  function dt!(du, u, p, t)
    x, y = u
    α, β, δ, γ = p
    du[1] = dx = α * x - β * x * y
    du[2] = dy = -δ * y + γ * x * y
  end

  function dW!(du, u, p, t)
    du[1] = 0.1u[1]
    du[2] = 0.1u[2]
  end

  u0 = [1.0, 1.0]
  tspan = (0.0, 10.0)
  p = [2.2, 1.0, 2.0, 0.4]
  prob_sde = SDEProblem(dt!, dW!, u0, tspan, p)

  condition(u, t, integrator) = integrator.t > 9.0 #some condition
  function affect!(integrator)
    #println("Callback")  #some callback
  end
  cb = DiscreteCallback(condition, affect!, save_positions=(false, false))

  function predict_sde(p)
    return Array(solve(prob_sde, EM(), p=p, saveat=savingtimes, sensealg=ForwardDiffSensitivity(), dt=0.001, callback=cb))
  end

  loss_sde(p) = sum(abs2, x - 1 for x in predict_sde(p))

  loss_sde(p)
  @time dp = gradient(p) do p
    loss_sde(p)
  end

  @test !iszero(dp[1])
end

@testset "SDEs" begin
  println("SDEs")
  test_SDE_callbacks()
end

ext: jl


text:
module SelfEnergies

export SelfEnergy

using ..BoundaryConditions
using ..Domains
using SparseArrays

struct SelfEnergy{N,M,TF}
    Σ0::NTuple{M,SparseMatrixCSC{ComplexF64,Int}}
    Σ1::NTuple{M,SparseMatrixCSC{ComplexF64,Int}}
    Σ2::SparseMatrixCSC{ComplexF64,Int}
    f::TF
end

include("1D/SelfEnergies.jl")
# include("2D/SelfEnergies.jl")
# include("3D/SelfEnergies.jl")

function Base.show(io::IO,Σ::SelfEnergy{N,M}) where {N,M}
    print(io,"SelfEnergy in $(N)D with $(M) parts")
end

end


text:
module linkedlist

    using DataStructures

    #=
    a few utility method for unit testing ease
    =#

    export are_equal

    function are_equal(l₁::LinkedList{T}, l₂::LinkedList{T})::Bool where T
        length(l₁) != length(l₂) && return false

        return collect(l₁) == collect(l₂)
    end

    export to_mutable_list

    function to_mutable_list(values::Vector{T}) where T
        list = MutableLinkedList()

        for value in values
            append!(list, value)
        end

        return list
    end

    #=
        Remove Dups.

        Write code to remove duplicates from an unsorted linked list.
    =#

    export dedupe

    #=
        dedupe 1

        if we get an immutable linked list, then
        we're forced to create a new one with the
        distinct elements from the passed in list.

        O(N) run-time
        O(N) memory
    =#
    function dedupe(list::LinkedList{T}) where T
        length(list) == 0 && return list

        seen = Set{T}()
        result = MutableLinkedList()

        for value in list

            value ∈ seen && continue

            append!(result, value)
            seen = union(seen, value)
        end

        return result
    end

    #=
        dedupe 2

        if we're passed a mutable linked list then
        we iterate over all the elements and mark the
        index of the duplicates in the first pass. Then
        we can iterate the indices of the duplicates and
        remove them seperately.

        O(N + D) run-time where D is the number of duplicates,
        O(1) memory
    =#
    function dedupe(list::MutableLinkedList{T}) where T

        length(list) == 0 && return list

        seen = Set{T}()
        to_delete = Vector{Integer}()

        for (index, value) in enumerate(list)
            value ∈ seen && append!(to_delete, index)
            seen = union(seen, value)
        end

        for (index, value) in enumerate(to_delete)
            position = value - (index - 1)
            delete!(list, position)
        end

        return list
    end

    #=
        Kth to Last.

        Implement an algorithm to find the kth to last
        element of a singly linked list.
    =#

    export kth_last

    #=
    kth last

    simple approach that uses two forward iterators, one
    as the main iterator and one as a runner that only starts
    iterating if the main iterator index is greater than k. Once
    the main iterator reaches the end of the list, the runner will
    be pointed to the kth last element, we just need to peel off
    the value.
    =#

    # O(N) where N = length of the list, O(1)/constant space.
    function kth_last(lst::LinkedList{T}, k::Integer) where T
        result, index = nothing, 0

        let iter = iterate(lst)
            let runner = iterate(lst)
                while iter !== nothing
                    (_, main_state) = iter
                    (_, run_state) = runner

                    if k < index
                        runner = iterate(lst, run_state)
                    end

                    iter = iterate(lst, main_state)
                    index += 1
                end

                result, _ = Iterators.peel(runner)
            end
        end

        return result
    end

    #=
        Delete Middle.

        Implement an algorithm to delete a node in the middle (i.e.,
        any node the first and last node, not necessarily the exact
        middle) of a singly linked list.

        EXAMPLE

        Input: the node c from the linked list
            a → b → c → d → e → f
        Result:
            a → b → d → e → f
    =#

    export delete_middle

    # O(n) where N → length of list
    function delete_middle(lst::LinkedList, k::Integer)
        T = typeof(lst).parameters[1]
        let iter = iterate(lst)
            index, new = 1, nil(T)
            while iter != nothing
                (head, tail) = iter
                index += 1
                if index <= k
                    new = cat(new, list(head))
                    iter = iterate(lst, tail)
                else
                    new = cat(new, tail)
                    break
                end
            end
            return new
        end
    end

    #=
        Partition.

        Write code to partition a linked list around a value x,
        such that all nodes less than x come before all nodes
        greater than or equal to x. If x is contained within the list,
        the values of x can appear anywhere in the "right partition";
        it does not need to appear between the left and right partitions.

        EXAMPLE

        Input: 3 → 5 → 8 → 5 → 10 → 2 → 1 [partition = 5]
        Output: 3 → 1 → 2 → 10 → 5 → 5 → 8
    =#

    export partition

    # O(n) where N → length of the list
    function partition(lst::LinkedList, value::Integer)::LinkedList
        let iter = iterate(lst)
            head, tail = nil(), nil()
            while iter != nothing
                v, s = iter

                if v < value
                    head = cons(v, head)
                else
                    tail = cons(v, tail)
                end

                iter = iterate(lst, s)
            end
            return cat(head, tail)
        end
    end

end # end module

max_stars_repo_name: bmoretz/Daily-Coding-Problem

max_issues_repo_name: bmoretz/Daily-Coding-Problem

max_forks_repo_name: bmoretz/Daily-Coding-Problem


text:
#=
  From https://math.stackexchange.com/questions/3633307/bayes-rule-broken
  """
  Bayes' Rule broken?!?!

  This question has been driving me CRAZY for 4 days now. The question comes from
  the textbook 'One Thousand Exercises in Probability', specifically Exercise 3 in
  section 1.4. The solution does not make sense! The question goes as follows:

     'A man possesses five coins, two double-headed, two normal and one double-tailed.
     The man shuts his eyes, picks a coin at random, and tosses the coin. He opens
     his eyes, sees a head: what is the probability the lower face is also a head?.

  The book gives an answer that is 2/3
  ....
  """

  Distributions of variable coin (num:0)
  1.00000 =>    6585  (0.658500)
  2.00000 =>    3415  (0.341500)

  Distributions of variable flip_is_head (num:0)
  1.00000 =>    6585  (0.658500)
  0.00000 =>    3415  (0.341500)

   Cf ~/blog/five_coins.blog
      ~/webppl/family_out_problem.wppl
=#

using Turing, StatsPlots, DataFrames
include("jl_utils.jl")

@model function five_coins()
    doubleHead = 1
    normalCoin = 2
    doubleTail = 3

    head = 1
    tail = 2

    coin ~ Categorical([2/5,2/5,1/5]) # [doubleHead,normalCoin,doubleTail]

    #  Toss the coin
    toss ~
        (coin == doubleHead) ? Dirac(head) :
        (coin == doubleTail) ? Dirac(tail) :
        Categorical([1/2,1/2]) # [head,tail]

    #  Flip the coin (i.e. turn it around)
    flip ~ Dirac((coin == doubleHead) ? head :
                 (coin == doubleTail) ? tail :
                 (toss == head) ? tail : head
                 )

    #  We observe that the toss give a head.
    true ~ Dirac(toss == head);

    #  What is the probability that the other side (flip side)
    #  is head?
    flip_is_head ~ Dirac(flip == head)
end

model = five_coins()

num_chains = 4

# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 10_000)
chns = sample(model, PG(5), 10_000)
# chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), 10_000)

display(chns)
# display(plot(chns))

show_var_dist_pct(chns,:coin)
show_var_dist_pct(chns,:flip_is_head)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
function scheduling(processing_times, due_dates, ::Val{:raw})
    m = model(; kind = :scheduling)
    n = length(processing_times) # number of jobs
    max_time = sum(processing_times)

    d = domain(0:max_time)
    foreach(_ -> variable!(m, d), 1:n) # C (1:n)
    foreach(_ -> variable!(m, d), 1:n) # S (n+1:2n)

    minus_eq_param(param) = x -> abs(x[1] - x[2] - param)
    less_than_param(param) = x -> abs(x[1] - param)
    sequential_tasks(x) = x[1] ≤ x[2] || x[3] ≤ x[4]

    for i in 1:n
        constraint!(m, minus_eq_param(processing_times[i]), [i, i + n])
        constraint!(m, less_than_param(due_dates[i]), [i])
    end
    for i in 1:n, j in 1:n
        i == j && continue
        constraint!(m, sequential_tasks, [i, j + n, j, i + n])
    end

    return m
end

function scheduling(processing_times, due_dates, ::Val{:JuMP})
    model = JuMP.Model(CBLS.Optimizer)
    n = length(processing_times) # number of jobs
    max_time = sum(processing_times)

    @variable(model, 0 ≤ C[1:n] ≤ max_time, Int) # completion
    @variable(model, 0 ≤ S[1:n] ≤ max_time, Int) # start

    for i in 1:n
        @constraint(model, [C[i], S[i]] in MinusEqualParam(processing_times[i]))
        @constraint(model, [C[i]] in LessThanParam(due_dates[i]))
    end
    for i in 1:n, j in 1:n
        i == j && continue
        @constraint(model, [C[i], S[j], C[j], S[i]] in SequentialTasks())
    end

    return model, C, S
end

"""
    scheduling(processing_time, due_date; modeler=:JuMP)

Create a model for the n-queens problem with `n` queens. The `modeler` argument accepts :JuMP (default), which refer to the JuMP model.

!!! warning

    The model seems to have a flaw. Needs to be investigated.

"""
scheduling(processing_times, due_dates; modeler=:JuMP) = scheduling(processing_times, due_dates, Val(modeler))


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

mutable struct ContainerServiceAgentPoolProfile <: SwaggerModel
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    count::Any # spec type: Union{ Nothing, Int32 } # spec name: count
    vmSize::Any # spec type: Union{ Nothing, String } # spec name: vmSize
    dnsPrefix::Any # spec type: Union{ Nothing, String } # spec name: dnsPrefix
    fqdn::Any # spec type: Union{ Nothing, String } # spec name: fqdn

    function ContainerServiceAgentPoolProfile(;name=nothing, count=1, vmSize=nothing, dnsPrefix=nothing, fqdn=nothing)
        o = new()
        validate_property(ContainerServiceAgentPoolProfile, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(ContainerServiceAgentPoolProfile, Symbol("count"), count)
        setfield!(o, Symbol("count"), count)
        validate_property(ContainerServiceAgentPoolProfile, Symbol("vmSize"), vmSize)
        setfield!(o, Symbol("vmSize"), vmSize)
        validate_property(ContainerServiceAgentPoolProfile, Symbol("dnsPrefix"), dnsPrefix)
        setfield!(o, Symbol("dnsPrefix"), dnsPrefix)
        validate_property(ContainerServiceAgentPoolProfile, Symbol("fqdn"), fqdn)
        setfield!(o, Symbol("fqdn"), fqdn)
        o
    end
end # type ContainerServiceAgentPoolProfile

const _property_map_ContainerServiceAgentPoolProfile = Dict{Symbol,Symbol}(Symbol("name")=>Symbol("name"), Symbol("count")=>Symbol("count"), Symbol("vmSize")=>Symbol("vmSize"), Symbol("dnsPrefix")=>Symbol("dnsPrefix"), Symbol("fqdn")=>Symbol("fqdn"))
const _property_types_ContainerServiceAgentPoolProfile = Dict{Symbol,String}(Symbol("name")=>"String", Symbol("count")=>"Int32", Symbol("vmSize")=>"String", Symbol("dnsPrefix")=>"String", Symbol("fqdn")=>"String")
Base.propertynames(::Type{ ContainerServiceAgentPoolProfile }) = collect(keys(_property_map_ContainerServiceAgentPoolProfile))
Swagger.property_type(::Type{ ContainerServiceAgentPoolProfile }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ContainerServiceAgentPoolProfile[name]))}
Swagger.field_name(::Type{ ContainerServiceAgentPoolProfile }, property_name::Symbol) =  _property_map_ContainerServiceAgentPoolProfile[property_name]

const _allowed_ContainerServiceAgentPoolProfile_vmSize = ["Standard_A0", "Standard_A1", "Standard_A2", "Standard_A3", "Standard_A4", "Standard_A5", "Standard_A6", "Standard_A7", "Standard_A8", "Standard_A9", "Standard_A10", "Standard_A11", "Standard_D1", "Standard_D2", "Standard_D3", "Standard_D4", "Standard_D11", "Standard_D12", "Standard_D13", "Standard_D14", "Standard_D1_v2", "Standard_D2_v2", "Standard_D3_v2", "Standard_D4_v2", "Standard_D5_v2", "Standard_D11_v2", "Standard_D12_v2", "Standard_D13_v2", "Standard_D14_v2", "Standard_G1", "Standard_G2", "Standard_G3", "Standard_G4", "Standard_G5", "Standard_DS1", "Standard_DS2", "Standard_DS3", "Standard_DS4", "Standard_DS11", "Standard_DS12", "Standard_DS13", "Standard_DS14", "Standard_GS1", "Standard_GS2", "Standard_GS3", "Standard_GS4", "Standard_GS5"]

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

function validate_property(::Type{ ContainerServiceAgentPoolProfile }, name::Symbol, val)
    if name === Symbol("count")
        Swagger.validate_param(name, "ContainerServiceAgentPoolProfile", :maximum, val, 100, false)
        Swagger.validate_param(name, "ContainerServiceAgentPoolProfile", :minimum, val, 1, false)
    end
    if name === Symbol("vmSize")
        Swagger.validate_param(name, "ContainerServiceAgentPoolProfile", :enum, val, _allowed_ContainerServiceAgentPoolProfile_vmSize)
    end
end


text:
@everywhere push!(LOAD_PATH,"solver")
@everywhere push!(LOAD_PATH,"mdp")
@everywhere using ParallelValueIteration
@everywhere using VerticalCAS
using HDF5

### OPTIONS ###
saveFile = "./VertCAS_qvals_parallel_v7_tauMax40.h5"
# Run in parallel on all available processors
nProcs = nprocs()-1
mdp = VerticalCAS_MDP()
approx_solver = ParallelValueIterationSolver(interp, verbose=true, max_iterations=1,n_procs=nProcs)

# Solve
println("Solving...")
qvals = solve(approx_solver, mdp)

# Write Qvalues to a file
println("Writing Qvalues")
h5open(saveFile, "w") do file
    write(file, "q", qvals)
end
end

max_stars_repo_name: bethelhall/CoCoCAS

max_issues_repo_name: bethelhall/CoCoCAS

max_forks_repo_name: bethelhall/CoCoCAS


text:
module TestExamples2DUnstructuredQuad

using Test
using Trixi

include("test_trixi.jl")

# pathof(Trixi) returns /path/to/Trixi/src/Trixi.jl, dirname gives the parent directory
EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "2d")

@testset "Unstructured Curve Mesh for Euler" begin

  @testset "elixir_euler_unstructured_quad_periodic.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_unstructured_quad_periodic.jl"),
      l2   = [0.00010978828464875207, 0.00013010359527356914, 0.00013010359527326057, 0.0002987656724828824],
      linf = [0.00638626102818618, 0.009804042508242183, 0.009804042508253286, 0.02183139311614468])
  end

  @testset "elixir_euler_unstructured_quad_free_stream.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_unstructured_quad_free_stream.jl"),
      l2   = [3.3596013848623337e-14, 2.4397818684917506e-13, 1.4659046683587473e-13, 4.673119220881569e-13],
      linf = [1.024846874031482e-11, 6.85111412046524e-11, 4.405278919428213e-11, 1.418225537008766e-10],
      tspan = (0.0, 0.1))
  end

  @testset "elixir_euler_unstructured_quad_wall_bc.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_unstructured_quad_wall_bc.jl"),
      l2   = [0.0401951408663404, 0.042562446022296446, 0.03734280979760256, 0.10058806400957201],
      linf = [0.24502910304890335, 0.298141188019722, 0.29446571031375074, 0.5937151600027155],
      tspan = (0.0, 0.25))
  end

  @testset "elixir_euler_unstructured_quad_basic.jl" begin
    @test_trixi_include(default_example_unstructured(),
      l2   = [0.0007258658867098887, 0.000676268065087451, 0.0006316238024054346, 0.0014729738442086392],
      linf = [0.004476908674416524, 0.0052614635050272085, 0.004926298866533951, 0.018058026023565432],
      tspan = (0.0, 1.0))
  end

  @testset "elixir_euler_unstructured_quad_basic_restart.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_unstructured_quad_basic_restart.jl"),
      l2   = [0.0007258658867098887, 0.000676268065087451, 0.0006316238024054346, 0.0014729738442086392],
      linf = [0.004476908674416524, 0.0052614635050272085, 0.004926298866533951, 0.018058026023565432])
  end

  @testset "elixir_euler_unstructured_quad_ec.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_unstructured_quad_ec.jl"),
      l2   = [0.06594600495903137, 0.10803914821786433, 0.10805946357846291, 0.1738171782368222],
      linf = [0.31880214280781305, 0.3468488554333352, 0.34592958184413264, 0.784555926860546],
      tspan = (0.0, 1.0))
  end
end

end # module


text:
# ------------------------------------------
# Adjoint julia script
# ------------------------------------------
using Kinetic, BenchmarkTools

begin
    u = collect(-5.0:0.1:5.0)
    nu = length(u)
    weights = ones(nu) .* 0.5

    fw = zeros(3)
    fh = zeros(nu)
    fb = zeros(nu)

    inK = 2
    γ = 5.0 / 3.0
    primL = [1., 0., 1.]
    wL = prim_conserve(primL, γ)
    hL = maxwellian(u, primL) |> Array;
    bL = hL .* 2 ./ (2.)
    shL = zeros(nu)
    sbL = zeros(nu)
    lenL = 0.1

    primR = [0.5, 0., 1.]
    wR = prim_conserve(primR, γ)
    hR = maxwellian(u, primR) |> Array;
    bR = hR .* 2 ./ (2.)
    shR = zeros(nu)
    sbR = zeros(nu)
    lenR = 0.1

    muref = 0.001
    omega = 0.72
    prandtl = 1.0
    dt = 1e-4
end

#--- kfvs ---#
#@btime
ccall(
    (:__kinetic_MOD_flux_kfvs_2f1v, "../src/fortran/kitmod.so"),
    Nothing,
    (
        Ref{Float64},
        Ref{Float64},
        Ref{Float64},
        Ref{Float64},
        Ref{Float64},
        Ref{Float64},
        Ref{Float64},
        Ref{Float64},
        Ref{Float64},
        Ref{Int},
        Ref{Float64},
        Ref{Float64},
        Ref{Float64},
        Ref{Float64},
        Ref{Float64},
    ),
    fw,
    fh,
    fb,
    hL,
    bL,
    hR,
    bR,
    u,
    weights,
    nu,
    dt,
    shL,
    sbL,
    shR,
    sbR,
)

#@btime
flux_kfvs!(fw, fh, fb, hL, bL, hR, bR, u, weights, dt, shL, sbL, shR, sbR)

#--- ugks ---#
#@btime
ccall(
    (:__kinetic_MOD_flux_ugks_2f1v, "../src/fortran/kitmod.so"),
    Nothing,
    (
        Ref{Float64},
        Ref{Float64},
        Ref{Float64},
        Ref{Float64},
        Ref{Float64},
        Ref{Float64},
        Ref{Float64},
        Ref{Float64},
        Ref{Float64},
        Ref{Float64},
        Ref{Float64},
        Ref{Int},
        Ref{Float64},
        Ref{Float64},
        Ref{Float64},
        Ref{Float64},
        Ref{Float64},
        Ref{Float64},
        Ref{Float64},
        Ref{Float64},
        Ref{Float64},
        Ref{Float64},
        Ref{Float64},
        Ref{Float64},
    ),
    fw,
    fh,
    fb,
    wL,
    hL,
    bL,
    wR,
    hR,
    bR,
    u,
    weights,
    nu,
    inK,
    γ,
    muref,
    omega,
    prandtl,
    dt,
    lenL,
    lenR,
    shL,
    sbL,
    shR,
    sbR,
)

#@btime
flux_ugks!(fw, fh, fb, wL, hL, bL, wR, hR, bR, u, weights, inK, γ, muref, omega, prandtl, dt, lenL, lenR, shL, sbL, shR, sbR)

"""
Results on my NUC with i7-8559U
@101 velocity points

Kinetic.jl
- kfvs ~ 6.747 μs (13 allocations: 11.38 KiB)
- ugks ~ 13.344 μs (123 allocations: 20.94 KiB) on my NUC

Kifort
- kfvs ~ 5.421 μs (37 allocations: 800 bytes)
- ugks ~ 11.413 μs (55 allocations: 1.09 KiB) on my NUC

"""

max_stars_repo_name: vavrines/Kifort

max_issues_repo_name: vavrines/Kifort

max_forks_repo_name: vavrines/Kifort


text:
function hist2spec(h::Histogram{<:Any, N}) where N
    ptype = [:BarPlot, :Heatmap, :Volume][N]
    f(edges) = edges[1:end-1] .+ diff(edges)./2
    kwargs = N == 1 ? (; width = step(h.edges[1])) : NamedTuple()
    return mapping(map(f, h.edges)..., Float64.(h.weights)) * visual(ptype; kwargs...)
end

to_weights(v) = weights(v)
to_weights(v::AbstractWeights) = v

function trim(range, min, max, closed)
    if closed == :left
        i1 = searchsortedlast(range, min)
        i2 = searchsortedfirst(range, nextfloat(max))
    else
        i1 = searchsortedlast(range, prevfloat(min))
        i2 = searchsortedfirst(range, max)
    end
    return range[i1:i2]
end

function compute_edges(data, extrema, bins::Tuple{Vararg{Integer}}, closed)
    ranges = map(extrema, bins) do (min, max), n
        histrange(min, max, n, closed)
    end
    # trim axis
    map(ranges, data) do range, d
        trim(range, Base.extrema(d)..., closed)
    end
end
compute_edges(data, extrema, bins::Tuple{Vararg{AbstractArray}}, closed) = bins

function _histogram(data...; bins = sturges(length(data[1])), wts = automatic,
    normalization = :none, extrema = map(extrema, data), closed = :left)

    bins_tuple = bins isa Tuple ? bins : map(_ -> bins, data)
    edges = compute_edges(data, extrema, bins_tuple, closed)
    weights = wts === automatic ? () : (to_weights(wts),)
    h = fit(Histogram, data, weights..., edges)
    hn = normalize(h, mode = normalization)
    return hist2spec(hn)
end

"""
    histogram(data...; bins = automatic, wts = automatic, normalization = :none)

Plot a histogram of `values`. `bins` can be an `Int` to create that
number of equal-width bins over the range of `values`.
Alternatively, it can be a sorted iterable of bin edges. The histogram
can be normalized by setting `normalization`. Possible values are:
*  `:pdf`: Normalize by sum of weights and bin sizes. Resulting histogram
   has norm 1 and represents a PDF.
* `:density`: Normalize by bin sizes only. Resulting histogram represents
   count density of input and does not have norm 1.
* `:probability`: Normalize by sum of weights only. Resulting histogram
   represents the fraction of probability mass for each bin and does not have
   norm 1.
*  `:none`: Do not normalize.
Weighted data is supported via the keyword `wts`.
"""
const histogram = Analysis(_histogram)

function global_options(::typeof(histogram), d::AlgebraicList)
    combine(extr1, extr2) = map(extr1, extr2) do (min1, max1), (min2, max2)
        min(min1, min2), max(max1, max2)
    end
    extr = mapfoldl(combine, d) do spec
        map(extrema, positional(spec.mapping.value))
    end
    return (extrema = extr,)
end


text:
using Pkg
Pkg.add("Images")
Pkg.add("ImageMagick")
Pkg.add("LinearAlgebra")
Pkg.add("Printf")
Pkg.add("FileIO")
Pkg.add("NLopt")
Pkg.add("SpecialFunctions")
Pkg.add("Gadfly")


text:
#       _/   _/   _/      _/ Kinh Nguyen
#      _/  _/    _/_/    _/  nguyenkinh@ytecongcong.com
#     _/ _/     _/ _/   _/
# TEXTS RUNNNING IN THE MIDDLE, NO, THIS IS NOT DEFAULT
#   _/ _/     _/   _/ _/
#  _/   _/   _/    _/_/
"""
Function to count for a cross-table for Epi package.
> Internal
"""
function crossCount(i, j, data) # Look, data dependent
	count(x->x==j, data[findall(x->x==i, data[:, 1]), 2])
end

"""
    removeMissing(A; cols=missing)

Removing missing values in an array `A`. By default remove all rows with missing values in **any of the columns**. Otherwise removing only missing data only in those columns specified by `cols`. `A` can be `R x 1` as well.

Return the modified `A`.
"""
function removeMissing(A; cols=missing)
	if ismissing(cols) cols=1:size(A, 2) end
	rowMissing = collect([findall(ismissing, A[:, i]) for i in cols])
	rowMissing = unlist(rowMissing, unify=true)
	A[setdiff(1:size(A, 1), rowMissing), :]
end

"""
    unlist(A; unify=false)

R's style 'unlist''. `A` can be a 1-dimensional array of arrays (the motivation of this function), or a 2-dimensional array (not tested exhaustively).

`unify=true` wraps the result in the `unique()` function before return.
"""
function unlist(A; unify=false)
	foo = Vector()
	for i in eachindex(A) append!(foo, A[i]) end
	unify ? unique(foo) : foo
end

"""
    removeNaN(A; cols=missing)

Removing non number values (`NaN`) in an array `A`. By default remove all rows with missing values in **any of the columns**. Otherwise removing only missing data only in those columns specified by `cols`. `A` can be `R x 1` as well

Return the modified `A`.
"""
function removeNaN(A; cols=missing)
	if ismissing(cols) cols=1:size(A, 2) end
	rowMissing = collect([findall(isnan, A[:, i]) for i in cols])
	rowMissing = unlist(rowMissing, unify=true)
	A[setdiff(1:size(A, 1), rowMissing), :]
end


text:
using LinearAlgebra

struct ValueOne; end
ValueOne()

# Compute X <- a X + b I.
function matfun_axpby!(X,a,b,Y::UniformScaling)
    m,n=size(X)
    if ~(a isa ValueOne)
        rmul!(X,a)
    end
    @inbounds for i=1:n
        X[i,i]+=(b isa ValueOne) ? 1 : b
    end
end

@inline function exp_ps_m7_opt_rho3_59(A)
    T=promote_type(eltype(A),Float64)
    A_copy=similar(A,T); A_copy .= A;
    return exp_ps_m7_opt_rho3_59!(A_copy)
end

@inline function exp_ps_m7_opt_rho3_59!(A)
    T=promote_type(eltype(A),Float64) # Make it work for many 'bigger' types (matrices and scalars)
    # max_memslots=10
    n=size(A,1)
    # The first slots are precomputed nodes [:A]
    memslots2 = similar(A,T)
    memslots3 = similar(A,T)
    memslots4 = similar(A,T)
    memslots5 = similar(A,T)
    memslots6 = similar(A,T)
    memslots7 = similar(A,T)
    memslots8 = similar(A,T)
    memslots9 = similar(A,T)
    memslots10 = similar(A,T)
    # Assign precomputed nodes memslots
    memslots1=A # overwrite A
    # Uniform scaling is exploited.
    # No matrix I explicitly allocated.
    # Computation order: Bb2 Ba2 B2 Ba3 Bb3 B3 Ba4 Bb4 B4 Ba5 Bb5 B5 Ba6 Bb6 B6 Ba7 Bb7 B7 Ba8 Bb8 B8 T2k10
    # Computing Bb2 = x*I+x*A
    coeff1=0.40288943046645087
    coeff2=0.3243718504289327
    memslots2 .= coeff2.*memslots1
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Ba2 = x*I+x*A
    coeff1=0.40288943046645087
    coeff2=0.3243718504289327
    memslots3 .= coeff2.*memslots1
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B2 with operation: mult
    mul!(memslots4,memslots3,memslots2)
    # Deallocating Ba2 in slot 3
    # Deallocating Bb2 in slot 2
    # Computing Ba3 = x*I+x*A+x*B2
    coeff1=0.07870552092616893
    coeff2=0.9039871669928049
    coeff3=0.3209056390860863
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb3 = x*I+x*A+x*B2
    coeff1=0.8599464743752632
    coeff2=0.29073009539473676
    coeff3=0.0877838768781113
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B3 with operation: mult
    mul!(memslots5,memslots2,memslots3)
    # Deallocating Ba3 in slot 2
    # Deallocating Bb3 in slot 3
    # Computing Ba4 = x*I+x*A+x*B2+x*B3
    coeff1=-0.07403200916806128
    coeff2=0.22059858577526156
    coeff3=0.6031002164926261
    coeff4=0.5359493824375996
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb4 = x*I+x*A+x*B2+x*B3
    coeff1=0.7557993691382382
    coeff2=0.16186458120340602
    coeff3=0.024014443649892782
    coeff4=0.010967880705307256
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B4 with operation: mult
    mul!(memslots6,memslots2,memslots3)
    # Deallocating Ba4 in slot 2
    # Deallocating Bb4 in slot 3
    # Computing Ba5 = x*I+x*A+x*B2+x*B3+x*B4
    coeff1=-0.00013699934593016893
    coeff2=0.14529476087609983
    coeff3=0.3261057160948171
    coeff4=0.710963649841826
    coeff5=0.34521134017746485
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb5 = x*I+x*A+x*B2+x*B3+x*B4
    coeff1=0.5028563878759467
    coeff2=0.5066127059710326
    coeff3=0.19548297448489585
    coeff4=0.04481784891151489
    coeff5=0.023615787302546792
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B5 with operation: mult
    mul!(memslots7,memslots2,memslots3)
    # Deallocating Ba5 in slot 2
    # Deallocating Bb5 in slot 3
    # Computing Ba6 = x*I+x*A+x*B2+x*B3+x*B4+x*B5
    coeff1=-9.869136468294225e-10
    coeff2=-7.512320521145497e-5
    coeff3=-0.0044353505165274365
    coeff4=-0.020309346530828727
    coeff5=0.195679851580968
    coeff6=0.9850416366804396
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb6 = x*I+x*A+x*B2+x*B3+x*B4+x*B5
    coeff1=9.211440361111637e-7
    coeff2=-1.3293424291101191e-8
    coeff3=-3.064863890686237e-9
    coeff4=7.122641395116041e-10
    coeff5=7.64822032891094e-12
    coeff6=1.3796783760471413e-13
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B6 with operation: mult
    mul!(memslots8,memslots2,memslots3)
    # Deallocating Ba6 in slot 2
    # Deallocating Bb6 in slot 3
    # Computing Ba7 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6
    coeff1=0.00041041842903290026
    coeff2=0.0029599104645556024
    coeff3=0.03421823700975471
    coeff4=-0.07655965971121473
    coeff5=-0.45442773831109645
    coeff6=0.9367257958392154
    coeff7=-4.0465312386473035e-8
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb7 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6
    coeff1=0.00660875329127505
    coeff2=9.739631514054463e-5
    coeff3=5.134838062997168e-5
    coeff4=1.9444505682154752e-5
    coeff5=9.194990346464199e-6
    coeff6=-8.932667565096493e-8
    coeff7=1.0000000001479474
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B7 with operation: mult
    mul!(memslots9,memslots2,memslots3)
    # Deallocating Ba7 in slot 2
    # Deallocating Bb7 in slot 3
    # Computing Ba8 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6+x*B7
    coeff1=0.7868458282010143
    coeff2=-0.044300593493478396
    coeff3=0.09816190377505705
    coeff4=0.45590771888971265
    coeff5=0.6380495620724085
    coeff6=0.7701122107742214
    coeff7=3.0616077790644636e-8
    coeff8=-0.040362632774053595
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8 .+ coeff8.*memslots9
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb8 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6+x*B7
    coeff1=0.1230391215271295
    coeff2=0.038029956154608906
    coeff3=0.010801521452913285
    coeff4=0.019497085762494887
    coeff5=0.0036174776941494966
    coeff6=-0.0021685857628971537
    coeff7=-1.2133490010511278e-5
    coeff8=1.002439604236121
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8 .+ coeff8.*memslots9
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B8 with operation: mult
    mul!(memslots10,memslots2,memslots3)
    # Deallocating Ba8 in slot 2
    # Deallocating Bb8 in slot 3
    # Computing T2k10 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6+x*B7+x*B8
    coeff1=0.7726601458443281
    coeff2=0.27826200797600853
    coeff3=0.3803045282099576
    coeff4=0.11581791434832714
    coeff5=0.3407182508161634
    coeff6=0.07858151588080117
    coeff7=1.6667149787880624e-7
    coeff8=0.00042398516888089884
    coeff9=1.009492506685841
    # Smart lincomb recycle A
    memslots1 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8 .+ coeff8.*memslots9 .+ coeff9.*memslots10
    mul!(memslots1,true,I*coeff1,true,true)
    # Deallocating B2 in slot 4
    # Deallocating B3 in slot 5
    # Deallocating B4 in slot 6
    # Deallocating B5 in slot 7
    # Deallocating B6 in slot 8
    # Deallocating B7 in slot 9
    # Deallocating B8 in slot 10
    return memslots1 # Returning T2k10
end

max_stars_repo_name: matrixfunctions/GraphMatFunData

max_issues_repo_name: matrixfunctions/GraphMatFunData

max_forks_repo_name: matrixfunctions/GraphMatFunData


text:
using FEniCS

mesh = UnitSquareMesh(8, 8)
element = FiniteElement("Lagrange", ufl_cell(mesh))

family_element = family(element)
cell_element = cell(element)
degree_element = degree(element)
variant_element = variant(element)
reconstructed_element = reconstruct(element)
sobolev_space_element = sobolev_space(element)
facet_element = FacetNormal(mesh)

V = FunctionSpace(mesh, "P", 1)
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(-6.0)

A_1 = inner(grad(u),grad(v))*dx
A_2 = dot(grad(u),grad(v))*dx
A_3 = cross(grad(u),grad(v))
A_4 = outer(grad(u),grad(v))
ng = nabla_grad(u)

B_1 = assemble(A_1)
B_2 = assemble(A_2)

measure_x = dx
measure_s = ds
measure_S = dS
measure_P = dP

true


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

module FakePTYs

if Sys.iswindows()
    pushfirst!(LOAD_PATH, Sys.STDLIB)
    using Sockets
    popfirst!(LOAD_PATH)
end

function open_fake_pty()
    @static if Sys.iswindows()
        # Fake being cygwin
        pid = string(getpid(), base=16, pad=16)
        pipename = """\\\\?\\pipe\\cygwin-$pid-pty10-abcdefg"""
        server = listen(pipename)
        pts = connect(pipename)
        @assert ccall(:jl_ispty, Cint, (Ptr{Cvoid},), pts.handle) == 1
        ptm = accept(server)
        close(server)
        # extract just the file descriptor
        fds = Libc.dup(Base._fd(pts))
        close(pts)
        pts = fds
        # convert pts handle to a TTY
        #fds = pts.handle
        #pts.status = Base.StatusClosed
        #pts.handle = C_NULL
        #pts = Base.TTY(fds, Base.StatusOpen)
    else
        O_RDWR = Base.Filesystem.JL_O_RDWR
        O_NOCTTY = Base.Filesystem.JL_O_NOCTTY

        fdm = ccall(:posix_openpt, Cint, (Cint,), O_RDWR | O_NOCTTY)
        fdm == -1 && error("Failed to open ptm")
        rc = ccall(:grantpt, Cint, (Cint,), fdm)
        rc != 0 && error("grantpt failed")
        rc = ccall(:unlockpt, Cint, (Cint,), fdm)
        rc != 0 && error("unlockpt")

        fds = ccall(:open, Cint, (Ptr{UInt8}, Cint),
            ccall(:ptsname, Ptr{UInt8}, (Cint,), fdm), O_RDWR | O_NOCTTY)

            pts = RawFD(fds)
        # pts = fdio(fds, true)
        # pts = Base.Filesystem.File(RawFD(fds))
        # pts = Base.TTY(RawFD(fds); readable = false)
        ptm = Base.TTY(RawFD(fdm))
    end
    return pts, ptm
end

function with_fake_pty(f)
    pts, ptm = open_fake_pty()
    try
        f(pts, ptm)
    finally
        close(ptm)
    end
    nothing
end

end

max_stars_repo_name: vanillajonathan/julia

max_issues_repo_name: vanillajonathan/julia

max_forks_repo_name: vanillajonathan/julia


text:
using MLFS
using PrettyPrint
using Setfield
using MLFS.HM
using MLFS.TypedIO
using MLStyle
using JSON

it = Symbol(:i, sizeof(Int) * 8)
g = empty(GlobalTC)
l = empty(LocalTC)
l = @set l.typeEnv = l.typeEnv[
    [
    :Int => T(Nom(it)),
    :Str => T(Nom(:str)),
    :Bool => T(Nom(:bool)),
    :Type => T(Nom(:Type))
    ]
]

stmts = @surf_toplevel begin
    implicitStr :: implicit[Str]
    implicitStr = "114514"
    implicitInt :: implicit[implicit[Str] -> Int]
    implicitInt = @julia(x -> parse(Int, x))
    (+) :: ((Int, Int) -> Int)
    (+) = @julia(((x, y),) -> x + y)
    fun :: (implicit[Int] -> Int -> Int)
    fun = x -> y -> x + y
    var :: Int
    var = fun(1)
end;
begin
    results = inferModule(g, l, stmts)[1];

    prune = g.tcstate.prune

    function postInfer(root::IR.Decl)
        @match root begin
            IR.Perform() => IR.gTrans(postInfer, root)
            IR.Assign(a, t, e) =>
                IR.Assign(a, prune(t), postInfer(e))
        end
    end

    function postInfer(root::IR.Expr)
        @match root begin
            IR.Expr(ln, ty, expr) =>
                let expr = postInferE(expr, ln),
                    ty = ty === nothing ? nothing : prune(ty)

                    IR.Expr(ln, ty, expr)
                end
        end
    end

    function postInferE(root::IR.ExprImpl, ln::LineNumberNode)
        @match root begin
            IR.EIm(expr, t, insts) =>
                let expr = postInfer(expr),
                    t = prune(t)

                    IR.EApp(expr, instanceResolve(g, insts, t, ln))
                end
            IR.ETypeVal(t) => IR.ETypeVal(prune(t))
            _ => IR.gTrans(postInfer, root)
        end
    end

    pprintln(results)
    pprintln(g.globalImplicits)
    solvedIRs = [postInfer(result) for result in results]
end;
ir = IR.ELet(solvedIRs,IR.Expr(LineNumberNode(1),nothing, IR.EInt(0, 8)))

open("testim_gen_dump.json", "w") do f
    write(f, json(toVec(IR.ExprImpl, ir)))
end

a = open("testim_gen_dump.json") do f
    fromVec(IR.ExprImpl, JSON.Parser.parse(read(f, String)))
end
println(a)

open("testim_gen.jl", "w") do f; println(f, "using FunctionWrappers; Fn = FunctionWrappers.FunctionWrapper"); println(f,Expr(:let, Expr(:block),
    irToJulia(ir, LineNumberNode(1))))
end


text:
immutable Point
  x::Int
  y::Int
end

t = Point(10, 10)
println(t)

function norm(p::Point)
  return sqrt(p.x^2+p.y^2)
end

println(norm(Point(10, 20)))


text:
struct ZeroOrderHoldLinearization{Dx,Du,T,TA<:SMatrix{Dx,Dx},TB<:SMatrix{Dx,Du},Tc<:SVector{Dx}}
    dt::T
    A::TA
    B::TB
    c::Tc
end
struct FirstOrderHoldLinearization{Dx,Du,T,TA<:SMatrix{Dx,Dx},TB0<:SMatrix{Dx,Du},TBf<:SMatrix{Dx,Du},Tc<:SVector{Dx}}
    dt::T
    A::TA
    B0::TB0
    Bf::TBf
    c::Tc
end

for (Tx, Tu) in ((:(StaticVector{Dx}), :(StaticVector{Du})),    # needed to avoid method ambiguity
                 (:(SVector{Dx}),      :(SVector{Du})))
    @eval function linearize(f::LinearDynamics{Dx,Du}, x::$Tx, u::$Tu;
                             keep_state_dims=:, keep_control_dims=:) where {Dx,Du}
        keep_state_dims === (:) && keep_control_dims === (:) && return f
        X, U = keep_state_dims, keep_control_dims
        Ã, B̃ = f.A[X,X], f.B[X,U]
        f̃, x̃, ũ = f(x,u)[X], x[X], u[U]
        LinearDynamics(Ã, B̃, f̃ - Ã*x̃ - B̃*ũ)
    end
end
function linearize(f::ZeroOrderHoldLinearization{Dx,Du}, x::StaticVector{Dx}, SC::StepControl{Du};
                   keep_state_dims=:, keep_control_dims=:) where {Dx,Du}
    @assert f.dt == SC.t
    keep_state_dims === (:) && keep_control_dims === (:) && return f
    X, U = keep_state_dims, keep_control_dims
    Ã, B̃ = f.A[X,X], f.B[X,U]
    f̃, x̃, ũ = f(x,SC)[X], x[X], SC.u[U]
    ZeroOrderHoldLinearization(f.dt, Ã, B̃, f̃ - Ã*x̃ - B̃*ũ)
end
function linearize(f::FirstOrderHoldLinearization{Dx,Du}, x::StaticVector{Dx}, RC::RampControl{Du};
                   keep_state_dims=:, keep_control_dims=:) where {Dx,Du}
    @assert f.dt == RC.t
    keep_state_dims === (:) && keep_control_dims === (:) && return f
    X, U = keep_state_dims, keep_control_dims
    Ã, B̃0, B̃f = f.A[X,X], f.B0[X,U], f.Bf[X,U]
    f̃, x̃, ũ0, ũf = f(x,RC)[X], x[X], RC.u0[U], RC.uf[U]
    FirstOrderHoldLinearization(f.dt, Ã, B̃0, B̃f, f̃ - Ã*x̃ - B̃0*ũ0 - B̃f*ũf)
end

for (Tx, TSC) in ((:(StaticVector{Dx}), :(StepControl{Du})),    # needed to avoid method ambiguity
                  (:(SVector{Dx}),      :(StepControl{Du,<:Any,<:SVector{Du}})))
    @eval function linearize(f::LinearDynamics{Dx,Du}, x::$Tx, SC::$TSC;
                             keep_state_dims=:, keep_control_dims=:) where {Dx,Du}
        X, U = keep_state_dims, keep_control_dims
        B, u = f.B[:,U], SC.u[U]
        eᴬᵗ, ∫eᴬᵗB = integrate_expAt_B(f.A, B, SC.t)
        _  , ∫eᴬᵗc = integrate_expAt_B(f.A, f.B*SC.u - B*u + f.c, SC.t)
        linearize(ZeroOrderHoldLinearization(SC.t, eᴬᵗ, ∫eᴬᵗB, ∫eᴬᵗc), x, StepControl(SC.t, u),
                  keep_state_dims=keep_state_dims)
    end
end
for (Tx, TRC) in ((:(StaticVector{Dx}), :(RampControl{Du})),    # needed to avoid method ambiguity
                  (:(SVector{Dx}),      :(RampControl{Du,<:Any,<:SVector{Du},<:SVector{Du}})))
    @eval function linearize(f::LinearDynamics{Dx,Du}, x::$Tx, RC::$TRC;
                             keep_state_dims=:, keep_control_dims=:) where {Dx,Du}
        if keep_control_dims === (:)
            u0, uf = RC.u0, RC.uf
            eᴬᵗ, ∫eᴬᵗB, ∫eᴬᵗBtdt⁻¹ = integrate_expAt_Bt_dtinv(f.A, f.B, RC.t)
            _  , ∫eᴬᵗc = integrate_expAt_B(f.A, f.c, RC.t)
            A  = eᴬᵗ
            B0 = ∫eᴬᵗBtdt⁻¹
            Bf = ∫eᴬᵗB - ∫eᴬᵗBtdt⁻¹
            c  = ∫eᴬᵗc
        else
            X, U = keep_state_dims, keep_control_dims
            B, u0, uf = f.B[:,U], RC.u0[U], RC.uf[U]
            b0 = f.B*RC.u0 - B*u0
            bf = f.B*RC.uf - B*uf
            B̂, û0, ûf = [B (bf - b0)], [u0; 0], [uf; 1]
            eᴬᵗ, ∫eᴬᵗB̂, ∫eᴬᵗB̂tdt⁻¹ = integrate_expAt_Bt_dtinv(f.A, B̂, RC.t)
            _  , ∫eᴬᵗc = integrate_expAt_B(f.A, f.c + b0, RC.t)
            B̂0 = ∫eᴬᵗB̂tdt⁻¹
            B̂f = ∫eᴬᵗB̂ - ∫eᴬᵗB̂tdt⁻¹
            A  = eᴬᵗ
            B0 = B̂0[:,SUnitRange(1,length(U))]
            Bf = B̂f[:,SUnitRange(1,length(U))]
            c  = ∫eᴬᵗc + B̂f[:,end]
        end
        linearize(FirstOrderHoldLinearization(RC.t, A, B0, Bf, c), x, RampControl(RC.t, u0, uf),
                  keep_state_dims=keep_state_dims)
    end
end

function linearize(f::DifferentialDynamics, x::State, u::Control; keep_state_dims=:, keep_control_dims=:)
    linearize(f, SVector(x...), SVector(u...),
              keep_state_dims=keep_state_dims, keep_control_dims=keep_control_dims)
end
function linearize(f::DifferentialDynamics, x::SVector, u::SVector;
                   keep_state_dims=:, keep_control_dims=:)
    A = ForwardDiff.jacobian(y -> f(y, u), x)    # TODO: incorporate state/control reduction before differentiating
    B = ForwardDiff.jacobian(w -> f(x, w), u)
    linearize(LinearDynamics(A, B, f(x,u) - A*x - B*u), x, u,
              keep_state_dims=keep_state_dims, keep_control_dims=keep_control_dims)
end

function linearize(f::DifferentialDynamics, x::StaticVector, SC::StepControl;
                   keep_state_dims=:, keep_control_dims=:)
    linearize(f, SVector(x), StepControl(SC.t, SVector(SC.u)),
              keep_state_dims=keep_state_dims, keep_control_dims=keep_control_dims)
end
function linearize(f::DifferentialDynamics, x::SVector, SC::StepControl{<:Any,<:Any,<:SVector};
                   keep_state_dims=:, keep_control_dims=:)
    A = ForwardDiff.jacobian(y -> propagate(f, y, SC), x)
    B = ForwardDiff.jacobian(w -> propagate(f, x, StepControl(SC.t, w)), SC.u)
    linearize(ZeroOrderHoldLinearization(SC.t, A, B, propagate(f, x, SC) - A*x - B*SC.u), x, SC,
              keep_state_dims=keep_state_dims, keep_control_dims=keep_control_dims)
end

function linearize(f::DifferentialDynamics, x::StaticVector, RC::RampControl;
                   keep_state_dims=:, keep_control_dims=:)
    linearize(f, SVector(x), RampControl(RC.t, SVector(RC.u0), SVector(RC.uf)),
              keep_state_dims=keep_state_dims, keep_control_dims=keep_control_dims)
end
function linearize(f::DifferentialDynamics, x::SVector, RC::RampControl{<:Any,<:Any,<:SVector,<:SVector};
                   keep_state_dims=:, keep_control_dims=:)
    A  = ForwardDiff.jacobian(y -> propagate(f, y, RC), x)
    B0 = ForwardDiff.jacobian(w -> propagate(f, x, RampControl(RC.t, w, RC.uf)), RC.u0)
    Bf = ForwardDiff.jacobian(w -> propagate(f, x, RampControl(RC.t, RC.u0, w)), RC.uf)
    linearize(FirstOrderHoldLinearization(RC.t, A, B0, Bf, propagate(f, x, RC) - A*x - B0*RC.u0 - Bf*RC.uf), x, RC,
              keep_state_dims=keep_state_dims, keep_control_dims=keep_control_dims)
end

function (f::ZeroOrderHoldLinearization{Dx,Du})(x::StaticVector{Dx}, SC::StepControl{Du}) where {Dx,Du}
    @assert f.dt == SC.t
    f.A*x + f.B*SC.u + f.c
end
function (f::FirstOrderHoldLinearization{Dx,Du})(x::StaticVector{Dx}, RC::RampControl{Du}) where {Dx,Du}
    @assert f.dt == RC.t
    f.A*x + f.B0*RC.u0 + f.Bf*RC.uf + f.c
end
(f::FirstOrderHoldLinearization{Dx,Du})(x::StaticVector{Dx}, SC::StepControl{Du}) where {Dx,Du} = f(RampControl(SC))


text:
abstract type  AbstractMaxEntropyModel   end

# Large version UBCM
struct UBCM{T} <: AbstractMaxEntropyModel where {T<:Real}
    x0::Vector{T}         # initial solution estimate
    x::Vector{T}          # solution if converged
    F::Vector{T}          # vector for values of the system of equations
    f!::Function          # sytem that needs to be solved
end

show(io::IO, model::UBCM{T}) where {T} = print(io, "UBCM ($(T)) model ($(length(model.x0)) nodes)")

function UBCM(G::T) where {T<:AbstractGraph}
    # assert validity of the method
    !is_connected(G) && throw(DomainError("Graph is not connected"))
    has_self_loops(G) && throw(DomainError("Graph has self loops"))
    # get degree sequence
    k = degree(G)
    # estimate for x0
    x0 = k ./ maximum(k)
    return UBCM(x0, similar(x0), similar(x0), (F::Vector, x::Vector) -> UBCM!(F, x, k))
end

function UBCM!(F::Vector, x::Vector, k::Vector)
    Threads.@threads for i in eachindex(x)
        @inbounds F[i] = -f(x[i]*x[i]) - k[i]
        @simd for j in eachindex(x)
            @inbounds F[i] += f(x[i]*x[j])
        end
    end
    return F
end

function solve!(model::UBCM)
    res =  nlsolve(model.f!, model.x0, autodiff=:forward)
    !res.f_converged && @warn "Solver did not converge"
    if res.f_converged
        any(res.zero .< 0) && @warn "Negative values in solution"
        @inbounds for i in eachindex(model.x)
            model.x[i] = res.zero[i]
        end
    end
    return res
end

# Compact version UBCM
mutable struct UBCMCompact{T, I} <: AbstractMaxEntropyModel where {T<:Real, I<:Integer}
    x0::Vector{T}         # initial solution estimate
    xs::Vector{T}         # solution if converged
    x::Union{IndirectArray, Nothing}      # expanded vector (using IndirectArrays)
    F::Vector{T}          # vector for values of the system of equations
    f!::Function          # sytem that needs to be solved
    #freqs::Dict{Int, Int} # dict holding the mapping (not required, but nice to have)
    #nodemap::Dict{Int, Int} # dict holding the mapping
    revmap::Vector{I}        # dict holding the reverse mapping
end

show(io::IO, model::UBCMCompact{T,I}) where {T,I} = print(io, "Compact UBCM ($(T),$(I)) model ($(length(model.revmap)) nodes reduced to $(length(model.x0)))")

function UBCMCompact(G::T) where {T<:AbstractGraph}
    # assert validity of the method
    !is_connected(G) && throw(DomainError("Graph is not connected"))
    has_self_loops(G) && throw(DomainError("Graph has self loops"))
    # get unique degrees and frequencies
    d = degree(G)
    freqs = countmap(d)
    κ = collect(keys(freqs))
    # get reverse mapping
    #revmap = Dict(i => findall(x->isequal(x, κ[i]), d) for i in eachindex(κ))
    #nodemap = Dict(keys => val for (key, val) in revmap)
    indtype = length(d) < typemax(UInt8) ? UInt8 : length(d) < typemax(UInt16) ? UInt16 : length(d) < typemax(UInt32) ? UInt32 : UInt64
    revmap = Vector{indtype}(undef, length(d))
	for i in eachindex(κ)
		revmap[findall(x->isequal(x, κ[i]), d)] .= indtype(i) # link values of large vector to compact one
	end

    # estimate for x0
    x0 = sqrt.(κ) ./ sum(κ)

    return UBCMCompact(x0, similar(x0), nothing, similar(x0), (F::Vector, x::Vector) -> UBCMCompact!(F, x, κ, freqs), revmap)
end

function UBCMCompact!(F::Vector, x::Vector, κ::Vector, freqs::Dict)
    Threads.@threads for i in eachindex(x)
        @inbounds F[i] = -f(x[i]*x[i]) - κ[i]
        @simd for j in eachindex(x)
            @inbounds F[i] += freqs[κ[j]]*f(x[i]*x[j])
        end
    end
    return F
end

function solve!(model::UBCMCompact)
    res =  nlsolve(model.f!, model.x0, autodiff=:forward)
    !res.f_converged && @warn "Solver did not converge"
    if res.f_converged
        # expansion of compact solution to large one happens here
        any(res.zero .< 0) && @warn "Negative values in solution"
        @inbounds for i in eachindex(model.x0)
            model.xs[i] = res.zero[i]
        end
        # generation of large solution happens here (using IndirectArrays)
        model.x = IndirectArray(model.revmap, model.xs)
    end
    return res
end

# supporting functions for computations of UBCM
f(x) = x / (1 + x) # helper function
a(i::Int,j::Int,x::T) where {T<:AbstractVector} = @inbounds f(x[i]*x[j]) # single element of adjacency matrix
"""Obtain full adjacency matrix from a model"""
function adjacency_matrix(model::UBCM)
    A = zeros(length(model.x), length(model.x))
    for i in eachindex(model.x)
        for j in eachindex(model.x)
            if i≠j
                A[i,j] = f(model.x[i]*model.x[j])
            end
        end
    end
    return A
end

"""helper function for adjacency matrix"""
function indmap(model::UBCMCompact,T::DataType,i,j)
    return model.revmap[i] <= model.revmap[j] ? T(model.revmap[i]+model.revmap[j]*(model.revmap[j]-1)÷2) : T(model.revmap[j]+model.revmap[i]*(model.revmap[i]-1)÷2)
end

"""Obtain full adjacency matrix from a compact model (as an IndirectArray)"""
function adjacency_matrix(model::UBCMCompact)
    # generate vector of all possible unique values in the adjacency matrix
    values = @inbounds [f(model.xs[i]*model.xs[j]) for j = 1:length(model.xs) for i = 1:j]
    push!(values,zero(eltype(values))) # include the zero for later
    # infer storage index type
    T = length(values) < typemax(UInt8) ? UInt8 : length(values) < typemax(UInt16) ? UInt16 : UInt32
    # generate indices in matrix from xs
    indices = [indmap(model,T,i,j) for i in 1:length(model.revmap), j in 1:length(model.revmap)]
    # set diagonal to zero
    indices[diagind(indices)] .= T(length(values))
    # return indirectarray
    return IndirectArray(indices, values)
end


text:
using ITensors

# TODO: this uses all of the tests to make
# precompile statements, but takes a long time
# (e.g. 700 seconds).
# Try again with later versions of PackageCompiler
#
# include(joinpath(joinpath(dirname(dirname(@__DIR__)),
#                           test"),
#                  "runtests.jl"))

N = 6
sweeps = Sweeps(3)
maxdim!(sweeps, 10)
cutoff!(sweeps, 1E-13)

ampo = OpSum()
for j in 1:(N - 1)
  ampo .+= "Sz", j, "Sz", j + 1
  ampo .+= 0.5, "S+", j, "S-", j + 1
  ampo .+= 0.5, "S-", j, "S+", j + 1
end

sites = siteinds("S=1", N)
H = MPO(ampo, sites)
psi0 = randomMPS(sites; linkdims=2)
dmrg(H, psi0, sweeps)

sites_qn = siteinds("S=1", N; conserve_qns=true)
if !hasqns(sites_qn[1])
  throw(ErrorException("Index does not have QNs in part of precompile script"))
end
H_qn = MPO(ampo, sites_qn)
psi0_qn = randomMPS(sites_qn, [isodd(n) ? "Up" : "Dn" for n in 1:N]; linkdims=2)
dmrg(H_qn, psi0_qn, sweeps)


text:
using LombScargle
using Measurements
using Base.Test

ntimes = 1001
# Observation times
t = linspace(0.01, 10pi, ntimes)
# Randomize times
t += step(t)*rand(ntimes)
t = collect(t)
# The signal
s = sinpi(t) + cospi(2t) + rand(ntimes)
# Frequency grid
nfreqs = 10000
freqs = linspace(0.01, 3, nfreqs)
# Randomize frequency grid
freqs += step(freqs)*rand(nfreqs)
freqs = collect(freqs)
# Use "freqpower" just to call that function and increase code coverage.
# "autofrequency" function is tested below.
@test_approx_eq_eps freqpower(lombscargle(t, s, frequencies=freqs, fit_mean=false))[2] freqpower(lombscargle(t, s, frequencies=freqs, fit_mean=true))[2] 5e-3

# Simple signal, without any randomness
t = collect(linspace(0.01, 10pi, ntimes))
s = sin(t)
pgram1 = lombscargle(t, s, fit_mean=false)
pgram2 = lombscargle(t, s, fit_mean=true)
@test_approx_eq_eps power(pgram1) power(pgram2) 2e-7
pgram3 = lombscargle(t, s, center_data=false, fit_mean=false)
pgram4 = lombscargle(t, s, center_data=false, fit_mean=true)
@test_approx_eq_eps power(pgram3) power(pgram4) 3e-7

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

# Test signal with uncertainties
err = collect(linspace(0.5, 1.5, ntimes))
@test_approx_eq power(lombscargle(t, s, err, frequencies=0.1:0.1:1, fit_mean=true))  [0.06659683848818687,0.09361438921056377,0.006815919926284516,0.0016347568319229223,0.0005385706045724484,0.00021180745624003642,0.00010539881897690343,7.01610752020905e-5,5.519295593372065e-5,4.339157565349008e-5]
@test_approx_eq power(lombscargle(t, s, err, frequencies=0.1:0.1:1, fit_mean=false)) [0.0692080444168825,0.09360343748833044,0.006634919855866448,0.0015362074096358814,0.000485825178683968,0.000181798596773626,8.543735242380012e-5,5.380000205539795e-5,4.010727072660524e-5,2.97840883747593e-5]
@test_approx_eq power(lombscargle(t, s, err, frequencies=0.1:0.1:1, fit_mean=false, center_data=false)) [0.06920814049261209,0.09360344864985352,0.006634919960009565,0.0015362072871144769,0.0004858250831632676,0.00018179850370583626,8.543727416919218e-5,5.379994730581837e-5,4.0107232867763e-5,2.9784059487535237e-5]
@test power(lombscargle(t, s, err)) ==
    power(lombscargle(t, measurement(s, err)))

# Test fast method
t = linspace(0.01, 10pi, ntimes)
s = sin(t)
pgram5 = lombscargle(t, s, maximum_frequency=30, fast=true)
pgram6 = lombscargle(t, s, maximum_frequency=30, fast=false)
@test_approx_eq_eps power(pgram5) power(pgram6) 2e-6
@test_approx_eq power(lombscargle(t, s, frequencies=0.2:0.2:1, fast=true, fit_mean=true))  [0.029886871262324963,0.0005453105325981758,1.9499330722168046e-5,2.0859593514888897e-6,1.0129019249708592e-5]
@test_approx_eq power(lombscargle(t, s, frequencies=0.2:0.2:1, fast=true, fit_mean=false)) [0.029886867760422008,0.0005453104197620392,1.9499329579010375e-5,2.085948496002562e-6,1.0128073536975395e-5]
@test_approx_eq power(lombscargle(t, s, frequencies=0.2:0.2:1, fast=true, center_data=false, fit_mean=false)) [0.029886868328967718,0.0005453105220405559,1.949931928224576e-5,2.0859802347505357e-6,1.0127777365273726e-5]
@test_approx_eq power(lombscargle(t, s, err, frequencies=0.2:0.2:1, fast=true, fit_mean=true))  [0.09230959166317655,0.0015654929813132702,0.00019405185108843607,6.0898671943944786e-5,6.0604505038256276e-5]
@test_approx_eq power(lombscargle(t, s, err, frequencies=0.2:0.2:1, fast=true, fit_mean=false)) [0.09219168665786258,0.0015654342453078724,0.00019403694017215876,6.088944186950046e-5,6.05771360378885e-5]
@test_approx_eq power(lombscargle(t, s, err, frequencies=0.2:0.2:1, fast=true, center_data=false, fit_mean=false)) [0.09360344864985332,0.0015354489715019735,0.0001784388515190763,4.744247354697125e-5,3.240223498703448e-5]
@test power(lombscargle(t, s, err)) ==
    power(lombscargle(t, measurement(s, err)))

##################################################
### Testing utilities

# Test findmaxfreq and findmaxpower
@test_approx_eq findmaxfreq(pgram1)        [31.997145470342]
@test_approx_eq findmaxfreq(pgram1, 0.965) [0.15602150741832602,31.685102455505348,31.997145470342,63.52622641842902,63.838269433265665]
t = linspace(0.01, 10pi, 1001)
s = sinpi(2t) + cospi(4t)
p = lombscargle(t, s, maximum_frequency=4)
@test_approx_eq findmaxfreq(p, [0.9, 1.1]) 1.0029954048320957
@test_approx_eq findmaxfreq(p, [1.9, 2.1]) 2.002806697267899
@test_approx_eq findmaxpower(pgram1) 0.9695017551608017

# Test autofrequency function
@test_approx_eq LombScargle.autofrequency(t)                       0.003184112396292367:0.006368224792584734:79.6824127172165
@test_approx_eq LombScargle.autofrequency(t, minimum_frequency=0)                   0.0:0.006368224792584734:79.6792286048202
@test_approx_eq LombScargle.autofrequency(t, maximum_frequency=10) 0.003184112396292367:0.006368224792584734:9.99492881196174

# Test probabilities and FAP
t = collect(linspace(0.01, 10pi, 101))
s = sin(t)
for norm in ("standard", "Scargle", "HorneBaliunas", "Cumming")
    P = lombscargle(t, s, normalization = norm)
    for z_0 in (0.1, 0.5, 0.9)
        @test_approx_eq prob(P, probinv(P, z_0)) z_0
        @test_approx_eq fap(P,  fapinv(P, z_0))  z_0
    end
end
P = lombscargle(t, s, normalization = "log")
@test_throws ErrorException prob(P, 0.5)
@test_throws ErrorException probinv(P, 0.5)

# Test model function
@test_approx_eq s LombScargle.model(t, s, 1/2pi, center_data=false, fit_mean=false)
s = sinpi(t) + pi*cospi(t) + e
@test_approx_eq s LombScargle.model(t, measurement(s, ones(s)), 0.5)

# Test extirpolation function
x = linspace(0, 10)
y = sin(x)
@test_approx_eq LombScargle.extirpolate(x, y) [0.39537718210649553,3.979484140636793,4.833090108345013,0.506805556164743,-3.828112427525919,-4.748341359084166,-1.3022050566901917,3.3367666084342256,5.070478111668922,1.291245296032218,-0.8264466821981216,0.0,0.0]
@test_approx_eq LombScargle.extirpolate(x, y, 11) LombScargle.extirpolate(x, y)[1:11]

# Test trig_sum
C, S = LombScargle.trig_sum(x, y, 1, 10)
@test_approx_eq S [0.0,0.3753570125888358,0.08163980192703546,-0.10139634351774979,-0.4334223744905633,-2.7843373311769777,0.32405810159838055,0.05729663600471602,-0.13191736591325876,-0.5295781583202946]
@test_approx_eq C [8.708141477890015,-0.5402668064176129,-0.37460815054027985,-0.3793457539084364,-0.5972351546196192,14.612204307982497,-0.5020253140297526,-0.37724493022381034,-0.394096831764578,-0.6828241623474718]


text:
abstract type AbstractFreeListAllocator{T, I} end
abstract type AbstractAllocator{T, I} <: AbstractFreeListAllocator{T, I} end
abstract type AbstractFreeListSOAllocator{T, I, S} <: AbstractFreeListAllocator{T, I} end
abstract type AbstractSOAllocator{T, I, S} <: AbstractFreeListSOAllocator{T, I, S} end

const N = 16

# Allocator
mutable struct Allocator{T, I} <: AbstractAllocator{T, I}
    n::I
    resizable::Bool
    store::Vector{T}
    first::I
    last::I
    Allocator{T, I}(::Nothing, init=0) where {T, I} =
        new{T, I}(N, true, Vector{T}(undef, N), init + 1, init)
    Allocator{T, I}(n::I, init=0) where {T, I} =
        new{T, I}(n, false, Vector{T}(undef, n), init + 1, init)
end

@inline Base.@propagate_inbounds Base.getindex(alloc::Allocator{T, I}, i::I) where {T, I} =
    @inbounds alloc.store[i]
@inline Base.@propagate_inbounds Base.setindex!(alloc::Allocator{T, I}, t::T, i::I) where {T, I} =
    @inbounds alloc.store[i] = t

@noinline function resizeend!(alloc::Allocator{T, I}, n::I) where {T, I}
    if n > alloc.n
        Base._growend!(alloc.store, n - alloc.n)
    elseif n < alloc.n
        Base._deleteend!(alloc.store, alloc.n - n)
    end
    alloc.n = n
end
@noinline function resizebegin!(alloc::Allocator{T, I}, n::I) where {T, I}
    if n > alloc.n
        Base._growbeg!(alloc.store, n - alloc.n)
    elseif n < alloc.n
        Base._deletebeg!(alloc.store, alloc.n - n)
    end
    alloc.first += n - alloc.n
    alloc.last += n - alloc.n
    alloc.n = n
end

@inline Base.@propagate_inbounds function allocateend!(alloc::Allocator{T, I}, t::T) where {T, I}
    if alloc.resizable && alloc.last >= alloc.n
        resizeend!(alloc, 2 * alloc.n)
    end
    alloc.last += 1
    @inbounds alloc.store[alloc.last] = t
    alloc.last
end
@inline Base.@propagate_inbounds function allocatebegin!(alloc::Allocator{T, I}, t::T) where {T, I}
    if alloc.resizable && alloc.first <= 1
        resizebegin!(alloc, 2 * alloc.n)
    end
    alloc.first -= 1
    @inbounds alloc.store[alloc.first] = t
    alloc.first
end
@inline function deallocateend!(alloc::Allocator{T, I}) where {T, I}
    alloc.last -= 1
    if alloc.resizable && max(N, 4 * alloc.last) < alloc.n
        resizeend!(alloc, alloc.n ÷ 2)
    end
end
@inline function deallocatebegin!(alloc::Allocator{T, I}) where {T, I}
    alloc.first += 1
    if alloc.resizable && max(N, 4 * (alloc.n - alloc.first)) < alloc.n
        resizebegin!(alloc, alloc.n ÷ 2)
    end
end
function deallocate!(alloc::Allocator{T, I}, i::I) where {T, I}
    # nope
end

@inline Base.isempty(alloc::Allocator{T, I}) where {T, I} =
    alloc.first == alloc.last + 1
@inline function emptyend!(alloc::Allocator{T, I}) where {T, I}
    alloc.first = 1
    alloc.last = 0
end
@inline function emptybegin!(alloc::Allocator{T, I}) where {T, I}
    alloc.first = alloc.n + 1
    alloc.last = alloc.n
end

# FreeListAllocator
struct FreeListAllocator{T, I} <: AbstractFreeListAllocator{T, I}
    used::Allocator{T, I}
    free::Allocator{I, I}
    FreeListAllocator{T, I}(::Nothing, init=0) where {T, I} =
        new{T, I}(Allocator{T, I}(nothing, init), Allocator{I, I}(nothing, 0))
    FreeListAllocator{T, I}(n::I, init=0) where {T, I} =
        new{T, I}(Allocator{T, I}(n, init), Allocator{I, I}(n, 0))
end

@inline Base.@propagate_inbounds Base.getindex(alloc::FreeListAllocator{T, I}, i::I) where {T, I} =
    @inbounds alloc.used[i]
@inline Base.@propagate_inbounds Base.setindex!(alloc::FreeListAllocator{T, I}, t::T, i::I) where {T, I} =
    @inbounds alloc.used[i] = t

@inline Base.@propagate_inbounds function allocateend!(alloc::FreeListAllocator{T, I}, t::T) where {T, I}
    if isempty(alloc.free)
        allocateend!(alloc.used, t)
    else
        # deallocateend!(alloc.free)
        @inbounds i = alloc.free[alloc.free.last]
        alloc.free.last -= 1
        @inbounds alloc.used[i] = t
        i
    end
end
@inline Base.@propagate_inbounds function allocatebegin!(alloc::FreeListAllocator{T, I}, t::T) where {T, I}
    if isempty(alloc.free)
        allocatebegin!(alloc.used, t)
    else
        # deallocateend!(alloc.free)
        @inbounds i = alloc.free[alloc.free.last]
        alloc.free.last -= 1
        @inbounds alloc.used[i] = t
        i
    end
end
#= @inline function deallocateend!(alloc::FreeListAllocator{T, I}) where {T, I}
    deallocateend!(alloc.used)
end
@inline function deallocatebegin!(alloc::FreeListAllocator{T, I}) where {T, I}
    deallocatebegin!(alloc.used)
end =#
@inline function deallocate!(alloc::FreeListAllocator{T, I}, i::I) where {T, I}
    allocateend!(alloc.free, i)
end

@inline function Base.isempty(alloc::FreeListAllocator{T, I}) where {T, I}
    alloc.free.last + 1 - alloc.free.first == alloc.used.last + 1 - alloc.used.first
end

function emptyend!(alloc::FreeListAllocator{T, I}) where {T, I}
    emptyend!(alloc.used)
    emptyend!(alloc.free)
end
function emptybegin!(alloc::FreeListAllocator{T, I}) where {T, I}
    emptybegin!(alloc.used)
    emptyend!(alloc.free)
end

# SOAllocator
mutable struct SOAllocator{T, I, S} <: AbstractSOAllocator{T, I, S}
    n::I
    resizable::Bool
    store::S
    first::I
    last::I
    function SOAllocator{T, I, S}(store::S, ::Nothing, init=0) where {T, I, S}
        @assert length(store) == N
        new{T, I, S}(N, true, store, init + 1, init)
    end
    function SOAllocator{T, I, S}(store::S, n::I, init=0) where {T, I, S}
        @assert length(store) == n
        new{T, I, S}(n, false, store, init + 1, init)
    end
end

@inline Base.@propagate_inbounds Base.getindex(alloc::SOAllocator{T, I, S}, i::I) where {T, I, S} =
    @inbounds alloc.store[i]
@inline Base.@propagate_inbounds Base.setindex!(alloc::SOAllocator{T, I, S}, t::T, i::I) where {T, I, S} =
    @inbounds alloc.store[i] = t
@inline Base.@propagate_inbounds Base.getindex(alloc::SOAllocator{T, I, S}, i::I, j::I) where {T <: Tuple, I, S <: TupleVector{T}} =
    @inbounds alloc.store.vectors[j][i]
@inline Base.@propagate_inbounds Base.setindex!(alloc::SOAllocator{T, I, S}, t::U, i::I, j::I) where {T <: Tuple, I, S <: TupleVector{T}, U} =
    @inbounds alloc.store.vectors[j][i] = t

@noinline function resizeend!(alloc::SOAllocator{T, I, S}, n::I) where {T, I, S}
    resizeend!(alloc.store, n)
    alloc.n = n
end
@noinline function resizebegin!(alloc::SOAllocator{T, I, S}, n::I) where {T, I, S}
    resizebegin!(alloc.store, n)
    alloc.first += n - alloc.n
    alloc.last += n - alloc.n
    alloc.n = n
end

@inline Base.@propagate_inbounds function allocateend!(alloc::SOAllocator{T, I, S}, t::T) where {T, I, S}
    if alloc.resizable && alloc.last >= alloc.n
        resizeend!(alloc, 2 * alloc.n)
    end
    alloc.last += 1
    @inbounds alloc.store[alloc.last] = t
    alloc.last
end
@inline Base.@propagate_inbounds function allocatebegin!(alloc::SOAllocator{T, I, S}, t::T) where {T, I, S}
    if alloc.resizable && alloc.first <= 1
        resizebegin!(alloc, 2 * alloc.n)
    end
    alloc.first -= 1
    @inbounds alloc.store[alloc.first] = t
    alloc.first
end
@inline function deallocateend!(alloc::SOAllocator{T, I, S}) where {T, I, S}
    alloc.last -= 1
    if alloc.resizable && max(N, 4 * alloc.last) < alloc.n
        resizeend!(alloc, alloc.n ÷ 2)
    end
end
@inline function deallocatebegin!(alloc::SOAllocator{T, I, S}) where {T, I, S}
    alloc.first += 1
    if alloc.resizable && max(N, 4 * (alloc.n - alloc.first)) < alloc.n
        resizebegin!(alloc, alloc.n ÷ 2)
    end
end
@inline function deallocate!(alloc::SOAllocator{T, I, S}, i::I) where {T, I, S}
    # nope
end

@inline Base.isempty(alloc::SOAllocator{T, I, S}) where {T, I, S} =
    alloc.first == alloc.last + 1
@inline function emptyend!(alloc::SOAllocator{T, I, S}) where {T, I, S}
    alloc.first = 1
    alloc.last = 0
end
@inline function emptybegin!(alloc::SOAllocator{T, I, S}) where {T, I, S}
    alloc.first = alloc.n + 1
    alloc.last = alloc.n
end

# FreeListSOAllocator
struct FreeListSOAllocator{T, I, S} <: AbstractFreeListSOAllocator{T, I, S}
    used::SOAllocator{T, I, S}
    free::Allocator{I, I}
    function FreeListSOAllocator{T, I, S}(used::S, ::Nothing, init=0) where {T, I, S}
        new{T, I, S}(SOAllocator{T, I, S}(used, nothing), Allocator{I, I}(nothing, 0))
    end
    function FreeListSOAllocator{T, I, S}(used::S, n::I, init=0) where {T, I, S}
        new{T, I, S}(SOAllocator{T, I, S}(used, n), Allocator{I, I}(n, 0))
    end
end

@inline Base.@propagate_inbounds Base.getindex(alloc::FreeListSOAllocator{T, I, S}, i::I) where {T, I, S} =
    @inbounds alloc.used[i]
@inline Base.@propagate_inbounds Base.setindex!(alloc::FreeListSOAllocator{T, I, S}, t::T, i::I) where {T, I, S} =
    @inbounds alloc.used[i] = t
@inline Base.@propagate_inbounds Base.getindex(alloc::FreeListSOAllocator{T, I, S}, i::I, j::I) where {T <: Tuple, I, S} =
    @inbounds alloc.used.store.vectors[j][i]
@inline Base.@propagate_inbounds Base.setindex!(alloc::FreeListSOAllocator{T, I, S}, t::U, i::I, j::I) where {T <: Tuple, I, S, U} =
    @inbounds alloc.used.store.vectors[j][i] = t

@inline Base.@propagate_inbounds function allocateend!(alloc::FreeListSOAllocator{T, I, S}, t::T) where {T, I, S}
    if isempty(alloc.free)
        allocateend!(alloc.used, t)
    else
        # deallocateend!(alloc.free)
        @inbounds i = alloc.free[alloc.free.last]
        alloc.free.last -= 1
        @inbounds alloc.used[i] = t
        i
    end
end
@inline Base.@propagate_inbounds function allocatebegin!(alloc::FreeListSOAllocator{T, I, S}, t::T) where {T, I, S}
    if isempty(alloc.free)
        allocatebegin!(alloc.used, t)
    else
        # deallocateend!(alloc.free)
        @inbounds i = alloc.free[alloc.free.last]
        alloc.free.last -= 1
        @inbounds alloc.used[i] = t
        i
    end
end
#= @inline function deallocateend!(alloc::FreeListSOAllocator{T, I, S}) where {T, I, S}
    deallocateend!(alloc.used)
end
@inline function deallocatebegin!(alloc::FreeListSOAllocator{T, I, S}) where {T, I, S}
    deallocatebegin!(alloc.used)
end =#
@inline function deallocate!(alloc::FreeListSOAllocator{T, I, S}, i::I) where {T, I, S}
    allocateend!(alloc.free, i)
end

@inline function Base.isempty(alloc::FreeListSOAllocator{T, I, S}) where {T, I, S}
    alloc.free.last + 1 - alloc.free.first == alloc.used.last + 1 - alloc.used.first
end

@inline function emptyend!(alloc::FreeListSOAllocator{T, I, S}) where {T, I, S}
    emptyend!(alloc.used)
    emptyend!(alloc.free)
end
@inline function emptybegin!(alloc::FreeListSOAllocator{T, I, S}) where {T, I, S}
    emptybegin!(alloc.used)
    emptyend!(alloc.free)
end


text:
using Distributions
using Random, Test

function test_location_scale_normal(μ::Float64,σ::Float64,μD::Float64,σD::Float64)
    ρ = Normal(μD,σD)
    d = LocationScale(μ,σ,ρ)
    dref = Normal(μ+σ*μD,σ*σD)

    @test minimum(d) == minimum(dref)
    @test maximum(d) == maximum(dref)
    @test extrema(d) == (minimum(d), maximum(d))

    #### Conversions

    @test convert(LocationScale{Float16}, μ, σ, ρ) == LocationScale{Float16,typeof(ρ)}(Float16(μ), Float16(σ), ρ)
    @test convert(LocationScale{Float16}, d) == LocationScale{Float16,typeof(ρ)}(Float16(μ), Float16(σ), ρ)

    @test location(d) == μ
    @test scale(d) == σ
    @test params(d) == (μ,σ,ρ)
    @test partype(d) == Float64

    #### Statistics

    @test mean(d) ≈ mean(dref)
    @test median(d) ≈ median(dref)
    @test mode(d) ≈ mode(dref)
    @test modes(d) ≈ modes(dref)

    @test var(d) ≈ var(dref)
    @test std(d) ≈ std(dref)
    @test skewness(d) ≈ skewness(dref)
    @test kurtosis(d) ≈ kurtosis(dref)

    @test isplatykurtic(d) == isplatykurtic(dref)
    @test isleptokurtic(d) == isleptokurtic(dref)
    @test ismesokurtic(d) == ismesokurtic(dref)

    @test entropy(d) ≈ entropy(dref)
    @test mgf(d,-0.1) ≈ mgf(dref,-0.1)

    #### Evaluation & Sampling

    insupport(d,0.4) == insupport(dref,0.4)
    @test pdf(d,0.1) ≈ pdf(dref,0.1)
    @test pdf.(d,2:4) ≈ pdf.(dref,2:4)
    @test logpdf(d,0.4) ≈ logpdf(dref,0.4)
    @test loglikelihood(d,[0.1,0.2,0.3]) ≈ loglikelihood(dref,[0.1,0.2,0.3])
    @test cdf(d,μ-0.4) ≈ cdf(dref,μ-0.4)
    @test logcdf(d,μ-0.4) ≈ logcdf(dref,μ-0.4)
    @test ccdf(d,μ-0.4) ≈ ccdf(dref,μ-0.4) atol=1e-100
    @test logccdf(d,μ-0.4) ≈ logccdf(dref,μ-0.4) atol=1e-16
    @test quantile(d,0.1) ≈ quantile(dref,0.1)
    @test quantile(d,0.5) ≈ quantile(dref,0.5)
    @test quantile(d,0.9) ≈ quantile(dref,0.9)
    @test cquantile(d,0.1) ≈ cquantile(dref,0.1)
    @test cquantile(d,0.5) ≈ cquantile(dref,0.5)
    @test cquantile(d,0.9) ≈ cquantile(dref,0.9)
    @test invlogcdf(d,log(0.2)) ≈ invlogcdf(dref,log(0.2))
    @test invlogcdf(d,log(0.5)) ≈ invlogcdf(dref,log(0.5))
    @test invlogcdf(d,log(0.8)) ≈ invlogcdf(dref,log(0.8))
    @test invlogccdf(d,log(0.2)) ≈ invlogccdf(dref,log(0.2))
    @test invlogccdf(d,log(0.5)) ≈ invlogccdf(dref,log(0.5))
    @test invlogccdf(d,log(0.8)) ≈ invlogccdf(dref,log(0.8))

    r = Array{Float64}(undef, 100000)
    rand!(d,r)
    @test mean(r) ≈ mean(dref) atol=0.01
    @test std(r) ≈ std(dref) atol=0.01
    @test cf(d, -0.1) ≈ cf(dref,-0.1)
    @test gradlogpdf(d, 0.1) ≈ gradlogpdf(dref, 0.1)

end

println("\ttesting LocationScale")
test_location_scale_normal(0.3,0.2,0.1,0.2)
test_location_scale_normal(-0.3,0.1,-0.1,0.3)
test_location_scale_normal(1.3,0.4,-0.1,0.5)


text:
function plot_transform_frame(transforms::Array, vecmesh::Array, frame::Int64)

    transforms = transform_steps([0 -1; 1 0]; n_steps=240)
    vecsₜ = [transforms[frame] * i for i in vecs]

    x = vec(first.(vecsₜ))
    y = vec(last.(vecsₜ))

    fig = Figure()
    ax = Axis(fig[1,1])

    hlines!(ax, 0, linewidth=3)
    vlines!(ax, 0, linewidth=3)

    scatter!(x, y, color = cols)

    xlims!(-10, 10)
    ylims!(-10, 10)

    fig
end


text:
module PokerDeck

export make_deck, ALL_CARDS_AS_STRING, draw!

using Random
using ..PokerCard

"""
Representing a deck. The first time we create, we seed the static
deck with the list of unique card integers.
"""
ALL_CARDS_AS_STRING = String[]

function make_deck()
    Deck = Integer[]
    for rank in PokerCard.STR_RANKS
        for suit in ['s', 'h', 'd', 'c']
            push!(ALL_CARDS_AS_STRING, rank * suit)
            push!(Deck, PokerCard.newCard(rank * suit))
        end
    end
    Random.shuffle!(Deck)
    return Deck
end

function draw!(number_of_cards::Integer, deck::Array)
    cards_drawn = Integer[]
    for i in 1:number_of_cards
        push!(cards_drawn, popfirst!(deck) )
    end
    return cards_drawn
end

end # end Deck

max_stars_repo_name: HigginTown/Quads

max_issues_repo_name: HigginTown/Quads

max_forks_repo_name: HigginTown/Quads


text:
function curvatureforces!(node::Int64, params::Parameters)
  # Static version
  # From top to bit (node == 1): :down, from bit to top: :up
  for kk in 1:2
    if !((kk == 1 && node == 1) || (kk == 2 && node == params.noofnodes))
      k = -3 + 2 * kk
      i2 = node + kk - 1
      i1 = i2 - 1
      l = params.nodes[i2].z - params.nodes[i1].z
      tn = params.elements[i1].ea * (params.del[3, 1, i2] - params.del[3, 1, i1]) / l
      ts = params.elements[i1].gj * (params.del[6, 1, i2] - params.del[6, 1, i1]) / l
      avthx = (params.del[4, 1, i2] + params.del[4, 1, i1]) / 2.0
      avthy = (params.del[5, 1, i2] + params.del[5, 1, i1]) / 2.0
      params.f[1, 1, node] += tn * avthy * k
      params.f[2, 1, node] -= tn * avthx * k
      params.f[4, 1, node] += ts * avthy * k
      params.f[5, 1, node] -= ts * avthx * k
    end
  end
end


text:
using Polynomials
using SemiLagrangian: getpoly
function test_cmplx(n, decref::CartesianIndex{2}=CartesianIndex(0,0))
#    fctref = Polynomial(big.(rationalize.(rand(n^2-1), tol=0.001)+ im*rationalize.(rand(n^2-1), tol=0.001)))
    tabref = [big(i+decref.I[1]+im*(j+decref.I[2])) + 0.1*(rand()+im*rand()-0.5-0.5*im) for i=1:n, j=1:n]
    moy = sum(tabref)/length(tabref)
    tabref = tabref .- moy

    p = getpoly(tabref)
    @show p
    for i=1:n, j=1:n
        res = p(tabref[i,j]) - i -im*j
#        @show n,decref,i,j,res
        @test isapprox(p( tabref[i,j] ), i +im*j, atol=1e-27)
    end
end

@testset "Lagrange complexe" begin
    for n=3:2:11
        println("")
        test_cmplx(n)
        test_cmplx(n, CartesianIndex(35,47))
        test_cmplx(n, CartesianIndex(-3557,4702))
    end
end


text:
module ScalarTime

import Dates

include("ExpressionSets.jl")
using .ExpressionSets

# import ExpressionSets: morse_code
export morse_code

export supdog

#export binary, senary, decimal, alphanumeric, modgin_int, modgin_fract

io = IOBuffer()

dt = Dates.format(Dates.now(), "H:MM e, Y-m-d")

greeting = "  'Sup, Dog?      --      $dt"

supdog() = print(io,greeting)

end    # module


text:
debug_level = lowercase(get(ENV, "MODERNGL_DEBUGGING", "false")) == "true"
open(joinpath(@__DIR__, "deps.jl"), "w") do io
    println(io, "const enable_opengl_debugging = $debug_level")
end


text:
#!/usr/bin/julia

# Copyright (C) 2021, 2022 Matheus Fernandes Bigolin <mfrdrbigolin@disroot.org>
# SPDX-License-Identifier: MIT

"Day Thirteen, Transparent Origami."

# This certainly was a enjoying day, especially part II.

include("Utils.jl")

using .Utils

"Fold the `dots` according to the `folds`’ instructions."
function fold_dots(dots :: Set{Tuple{Int64, Int64}}, folds :: Array{Tuple{Char, Int64}}) :: Set{Tuple{Int64, Int64}}
    folded_dots = Set{Tuple{Int, Int}}(copy(dots))

    for (axis, pos) in folds
        for (x, y) in folded_dots
            if axis == 'x' && x > pos
                push!(folded_dots, (2*pos - x, y))
                delete!(folded_dots, (x,y))
            elseif axis == 'y' && y > pos
                push!(folded_dots, (x, 2*pos - y))
                delete!(folded_dots, (x,y))
            end
        end
    end

    return folded_dots
end

"Calculate the number of `dots` after one `fold`."
function solve1(dots :: Set{Tuple{Int64, Int64}}, fold :: Tuple{Char, Int64}) :: Int64
    folded_dots = fold_dots(dots, [fold])

    return length(folded_dots)
end

"Apply the `folds` to the `dots` and print the secret code."
function solve2(dots :: Set{Tuple{Int64, Int64}}, folds :: Array{Tuple{Char, Int64}})
    folded_dots = fold_dots(dots, folds)

    max_x, max_y = maximum(first.(folded_dots)), maximum(last.(folded_dots))

    for i in 0:max_y
        for j in 0:max_x
            if (j,i) in folded_dots
                print('█')
            else
                print('⋅')
            end
        end

        println()
    end
end

const DOTS_REX = r"(\d+),(\d+)"
const FOLDS_REX = r"fold along (\w)=(\d+)"

function main()
    usage_and_exit(length(ARGS) != 1)

    input_data = readlines(first(ARGS))

    dots_data = filter(line -> !startswith(line, "fold along") && line != "", input_data)
    folds_data = filter(line -> startswith(line, "fold along"), input_data)

    dots_untyped = [match(DOTS_REX, dot).captures for dot in dots_data]
    dots = Set([Tuple(parse.(Int, dot)) for dot in dots_untyped])

    folds_untyped = [match(FOLDS_REX, fold).captures for fold in folds_data]
    folds = [(collect(axis)[1], parse(Int, pos)) for (axis, pos) in folds_untyped]

    println(solve1(dots, folds[1]))
    solve2(dots, folds)
end

if abspath(PROGRAM_FILE) == @__FILE__
    main()
end


text:
include("header.jl")

using SpecialFunctions
using AutoGrad: lgamma, loggamma # TODO: delete after everyone has SpecialFunctions 0.8

@testset "specialfunctions" begin
    o = (:delta=>0.0001,:rtol=>0.02,:atol=>0.01)
    ϵ = 0.1
    val_0_2(x)=rand() * (2-2ϵ) + ϵ
    val_gt_0(x)=abs(x) + ϵ
    val_gt_m1(x)=abs(x) - 1 + ϵ
    val_lt_1(x)=-abs(x) + 1 - ϵ
    abs_lt_1(x)=rand() * (2-2ϵ) - (1-ϵ)
    val_gamma(x)=(x < ϵ && abs(x-round(x)) < ϵ ? x+0.5 : x) # avoid <=0 integers

    # `airy(k,x)` is deprecated, use `airyai(x)`, `airyaiprime(x)`, `airybi(x)` or `airybiprime(x)` instead.
    @test randcheck(airyai; o...) # @primitive airyai(x),dy (dy.*(airyaiprime.(x)))
    @test randcheck(airyaiprime; o...) # @primitive airyaiprime(x),dy (dy.*(x .* airyai.(x)))
    # airyaiprimex
    # airyaix
    @test randcheck(airybi,val_lt_1; o...) # @primitive airybi(x),dy (dy.*(airybiprime.(x)))
    @test randcheck(airybiprime,val_lt_1; o...) # @primitive airybiprime(x),dy (dy.*(x .* airybi.(x)))
    # airybiprimex
    # airybix
    # `airyprime(z::Number)` is deprecated, use `airyaiprime(z)` instead.
    # `airyx(k,x)` is deprecated, use `airyaix(x)`, `airyaiprimex(x)`, `airybix(x)` or `airybiprimex(x)` instead.
    # besselh
    # besselhx
    # besseli
    # besselix
    # besselj
    @test randcheck(besselj0; o...) # @primitive besselj0(x),dy (dy.*(-(besselj1.(x))))
    @test randcheck(besselj1; o...) # @primitive besselj1(x),dy (dy.*((besselj0.(x) - besselj.(2,x)) / 2))
    # besseljx
    # besselk
    # besselkx
    # bessely
    @test randcheck(bessely0,val_gt_0; o...) # @primitive bessely0(x),dy (dy.*(-(bessely1.(x))))
    @test randcheck(bessely1,val_gt_0; o...) # @primitive bessely1(x),dy (dy.*((bessely0.(x) - bessely.(2,x)) / 2))
    # besselyx
    # beta
    # cosint
    @test randcheck(dawson; o...) # @primitive dawson(x),dy,y (dy.*((-2y) .* x + 1))
    @test randcheck(digamma,val_gamma; o...) # @primitive digamma(x),dy,y (dy.*(trigamma.(x))) ## avoid <=0 ints
    @test randcheck(erf; o...) # @primitive erf(x),dy,y (dy.*(exp.(-(abs2.(x))) * convert(eltype(x), 2 / √π)))
    @test randcheck(erfc; o...) # @primitive erfc(x),dy,y (dy.*(-(exp.(-(abs2.(x)))) * convert(eltype(x), 2 / √π)))
    @test randcheck(erfcinv,val_0_2; o...) # @primitive erfcinv(x),dy,y (dy.*(-(exp.(abs2.(y))) * convert(eltype(x), √π / 2)))
    @test randcheck(erfcx,val_gt_m1; o...) # @primitive erfcx(x),dy,y (dy.*((2y) .* x - convert(eltype(x), 2 / √π)))
    @test randcheck(erfi,abs_lt_1; o...) # @primitive erfi(x),dy,y (dy.*(exp.(abs2.(x)) * convert(eltype(x), 2 / √π)))
    @test randcheck(erfinv,abs_lt_1; o...) # @primitive erfinv(x),dy,y (dy.*(exp.(abs2.(y)) * convert(eltype(x), √π / 2)))
    # eta
    @test randcheck(gamma,val_gamma; o...) # @primitive gamma(x),dy,y (dy.*(y .* digamma.(x))) ## avoid <=0 ints
    # hankelh1
    # hankelh1x
    # hankelh2
    # hankelh2x
    @test randcheck(invdigamma, val_lt_1; o...) # @primitive invdigamma(x),dy,y (dy.*(1 ./ trigamma.(y)))
    # lbeta
    # `lfact` is deprecated, use `lfactorial` instead.
    # lfactorial
    @test randcheck(loggamma,val_gt_0; o...) # @primitive loggamma(x),dy,y (dy.*(digamma.(x))) ## avoid <=0
    # @test randcheck(polygamma) # @primitive polygamma(x1,x2),dy,y  nothing  unbroadcast(x2,dy.*polygamma(x1+1,x2))
    # sinint
    @test randcheck(trigamma,val_gamma; o...) # @primitive trigamma(x),dy,y (dy.*(polygamma.(2,x)))
    # zeta
end


text:
using EchelleCCFs
using RvSpectMLBase
using Test

@testset "Convenience functions" begin  # TODO repeat for multiple mask shapes
    Δv = 1000
    λ1 = 5025.0
    λ2 = 5075.0
    d1 = 0.5
    d2 = 0.25
    σ_v = 5000.0
    λlo = 5000
    λhi = 5100
    npix = 1000
    @test_nowarn BasicLineList([λ1,λ2], [d1,d2])
    ll = BasicLineList([λ1,λ2], [d1,d2])
    @test_nowarn BasicCCFPlan(line_list=ll, mask_shape=TopHatCCFMask(Δv))
    ccfpl = BasicCCFPlan(line_list=ll, mask_shape=TopHatCCFMask(Δv))
    @test_nowarn calc_ccf_v_grid(ccfpl)
    v_grid = calc_ccf_v_grid(ccfpl)
    @test length(v_grid) == calc_length_ccf_v_grid(ccfpl)
    λs = exp.(range(log(λlo),stop=log(λhi),length=npix))
    flux = ones(size(λs))
    flux .*= 1 .- d1.*exp.(-0.5.*((λs.-λ1)./λ1.*(EchelleCCFs.speed_of_light_mps./σ_v)).^2)
    flux .*= 1 .- d2.*exp.(-0.5.*((λs.-λ2)./λ2.*(EchelleCCFs.speed_of_light_mps./σ_v)).^2)
    @testset "calc_ccf_chunk" begin
        chunk = ChunkOfSpectrum(λs,flux,ones(size(flux)))
        @test_nowarn calc_ccf_chunk(chunk, ccfpl)
        @test_nowarn calc_ccf_and_var_chunk(chunk, ccfpl)
    end
    @testset "calc_ccf_chunklist" begin
        chunk = ChunkOfSpectrum(λs,flux,ones(size(flux)))
        chunklist = ChunkList([chunk,chunk,chunk],[1,2,3])
        @test_nowarn calc_ccf_chunklist(chunklist, fill(ccfpl,3) )
        @test_nowarn calc_ccf_and_var_chunklist(chunklist, fill(ccfpl,3) )
    end
    @testset "calc_order_ccfs_chunklist" begin
        chunk = ChunkOfSpectrum(λs,flux,ones(size(flux)))
        chunklist = ChunkList([chunk,chunk,chunk],[1,2,3])
        @test_nowarn calc_order_ccfs_chunklist(chunklist, fill(ccfpl,3)  )
        @test_nowarn calc_order_ccf_and_vars_chunklist(chunklist, fill(ccfpl,3)  )
    end

    @testset "Compare calc_ccfs_chunklist and calc_order_ccfs_chunklist" begin
        chunk = ChunkOfSpectrum(λs,flux,ones(size(flux)))
        chunklist = ChunkList([chunk,chunk,chunk],[1,2,3])
        ccf = calc_ccf_chunklist(chunklist, fill(ccfpl,3) )
        order_ccf = calc_order_ccfs_chunklist(chunklist, fill(ccfpl,3)  )
        num_lines = 6
        @test all(ccf .≈ sum(order_ccf,dims=2).*num_lines)
    end

    #= Not working
    @testset "calc_ccf_chunklist_timeseries" begin
        chunk = ChunkOfSpectrum(λs,flux,ones(size(flux)))
        chunklist = ChunkList([chunk,chunk,chunk],[1,2,3])
        #clt = ChunkListTimeseries(rand(3), [chunklist,chunklist,chunklist] )
        #@test_nowarn calc_ccf_chunklist_timeseries(clt, fill(ccfpl,3)  )
    end
    =#

    #=  Not working
    @testset "calc_ccf_chunklist_timeseries" begin
        using RvSpectMLBase.TheoreticalInstrument
        using DataFrames
        num_pixels = 1000
        num_orders = 2
        num_obs = 3
        lambda2d = collect(reshape(range(λlo,stop=λhi,length=num_pixels*num_orders),num_pixels,num_orders))
        inst = TheoreticalInstrument2D(lambda2d)
        times = sort(rand(num_obs))
        rvs = range(0.0, stop=10.0, length=length(times))
        line_list = DataFrame(:lambda=>[λ1,λ2], :weight=>[d1,d2])
        @test_nowarn TheoreticalInstrument.generate_spectra_timeseries(times,line_list,inst, rvs)
        spectra = TheoreticalInstrument.generate_spectra_timeseries(times,line_list,inst, rvs)
        #chunk_list = DataFrame(:lambda_lo=>[λlo, (λlo+λhi)/2], :lambda_hi=>[(λlo+λhi)/2, λhi])
        chunk_list = DataFrame(:lambda_lo=>0.9.*[λ1, λ2], :lambda_hi=>1.1.*[λ1, λ2])
        @test_nowarn make_chunk_list_timeseries(spectra,chunk_list)
        clt = make_chunk_list_timeseries(spectra,chunk_list)
        @test_nowarn calc_ccf_chunklist_timeseries(clt, fill(ccfpl,num_obs)  )
    end
    =#

end

max_stars_repo_name: RvSpectML/EchelleCCFs

max_issues_repo_name: RvSpectML/EchelleCCFs

max_forks_repo_name: RvSpectML/EchelleCCFs