text:
using JuMP
using GLPK

# Create a model
sudoku = Model(with_optimizer(GLPK.Optimizer))

# Create our variables
@variable(sudoku, x[i=1:9, j=1:9, k=1:9], Bin);

for i = 1:9, j = 1:9  # Each row and each column
    # Sum across all the possible digits
    # One and only one of the digits can be in this cell,
    # so the sum must be equal to one
    @constraint(sudoku, sum(x[i,j,k] for k in 1:9) == 1)
end

for ind = 1:9  # Each row, OR each column
    for k = 1:9  # Each digit
        # Sum across columns (j) - row constraint
        @constraint(sudoku, sum(x[ind,j,k] for j in 1:9) == 1)
        # Sum across rows (i) - column constraint
            @constraint(sudoku, sum(x[i,ind,k] for i in 1:9) == 1)
    end
end

for i = 1:3:7, j = 1:3:7, k = 1:9
    # i is the top left row, j is the top left column
    # We'll sum from i to i+2, e.g. i=4, r=4, 5, 6
    @constraint(sudoku, sum(x[r,c,k] for  r in i:i+2, c in j:j+2) == 1)
end

# The given digits
init_sol = [ 5 3 0 0 7 0 0 0 0;
             6 0 0 1 9 5 0 0 0;
             0 9 8 0 0 0 0 6 0;
             8 0 0 0 6 0 0 0 3;
             4 0 0 8 0 3 0 0 1;
             7 0 0 0 2 0 0 0 6;
             0 6 0 0 0 0 2 8 0;
             0 0 0 4 1 9 0 0 5;
             0 0 0 0 8 0 0 7 9]
for i = 1:9, j = 1:9
    # If the space isn't empty
    if init_sol[i,j] != 0
        # Then the corresponding variable for that digit
        # and location must be 1
        @constraint(sudoku, x[i,j,init_sol[i,j]] == 1)
    end
end

# solve problem
optimize!(sudoku)

# Extract the values of x
x_val = value.(x)
# Create a matrix to store the solution
sol = zeros(Int,9,9)  # 9x9 matrix of integers
for i in 1:9, j in 1:9, k in 1:9
    # Integer programs are solved as a series of linear programs
    # so the values might not be precisely 0 and 1. We can just
    # round them to the nearest integer to make it easier
    if round(Int,x_val[i,j,k]) == 1
        sol[i,j] = k
    end
end
# Display the solution
sol


text:
### Edge and Gradient Related Image Operations ###

# Phase (angle of steepest gradient ascent), calculated from X and Y gradient images
"""
    phase(grad_x, grad_y) -> p

Calculate the rotation angle of the gradient given by `grad_x` and
`grad_y`. Equivalent to `atan(-grad_y, grad_x)`, except that when both `grad_x` and
`grad_y` are effectively zero, the corresponding angle is set to zero.
"""
function phase(grad_x::T, grad_y::T, tol=sqrt(eps(T))) where T<:Number
    atan(-grad_y, grad_x) * ((abs(grad_x) > tol) | (abs(grad_y) > tol))
end
phase(grad_x::Number,   grad_y::Number)   = phase(promote(grad_x, grad_y)...)
phase(grad_x::NumberLike, grad_y::NumberLike) = phase(gray(grad_x), gray(grad_y))

phase(grad_x::AbstractRGB, grad_y::AbstractRGB) = phase(vecsum(grad_x), vecsum(grad_y))

magnitude_phase(grad_x::NumberLike, grad_y::NumberLike) =
    hypot(grad_x, grad_y), phase(grad_x, grad_y)

function magnitude_phase(grad_x::AbstractRGB, grad_y::AbstractRGB)
    gx, gy = vecsum(grad_x), vecsum(grad_y)
    magnitude_phase(gx, gy)
end

vecsum(c::AbstractRGB) = float(red(c)) + float(green(c)) + float(blue(c))

## TODO? orientation seems nearly redundant with phase, deprecate?

"""
    orientation(grad_x, grad_y) -> orient

Calculate the orientation angle of the strongest edge from gradient images
given by `grad_x` and `grad_y`.  Equivalent to `atan(grad_x, grad_y)`.  When
both `grad_x` and `grad_y` are effectively zero, the corresponding angle is set to
zero.
"""
function orientation(grad_x::T, grad_y::T, tol=sqrt(eps(T))) where T<:Number
    atan(grad_x, grad_y) * ((abs(grad_x) > tol) | (abs(grad_y) > tol))
end
orientation(grad_x::Number,   grad_y::Number)   = orientation(promote(grad_x, grad_y)...)
orientation(grad_x::NumberLike, grad_y::NumberLike) = orientation(gray(grad_x), gray(grad_y))

orientation(grad_x::AbstractRGB, grad_y::AbstractRGB) = orientation(vecsum(grad_x), vecsum(grad_y))

# Magnitude of gradient, calculated from X and Y image gradients
"""
```
m = magnitude(grad_x, grad_y)
```

Calculates the magnitude of the gradient images given by `grad_x` and `grad_y`.
Equivalent to `sqrt(grad_x.^2 + grad_y.^2)`.

Returns a magnitude image the same size as `grad_x` and `grad_y`.
"""
magnitude(grad_x::AbstractArray, grad_y::AbstractArray) = hypot.(grad_x, grad_y)

Base.hypot(x::AbstractRGB, y::AbstractRGB) = hypot(vecsum(x), vecsum(y))

phase(grad_x::AbstractArray, grad_y::AbstractArray) = phase.(grad_x, grad_y)

# Orientation of the strongest edge at a point, calculated from X and Y gradient images
# Note that this is perpendicular to the phase at that point, except where
# both gradients are close to zero.

orientation(grad_x::AbstractArray{T}, grad_y::AbstractArray{T}) where {T} = orientation.(grad_x, grad_y)

# Return both the magnitude and phase in one call
"""
    magnitude_phase(grad_x, grad_y) -> m, p

Convenience function for calculating the magnitude and phase of the gradient
images given in `grad_x` and `grad_y`.  Returns a tuple containing the magnitude
and phase images.  See `magnitude` and `phase` for details.
"""
function magnitude_phase(grad_x::AbstractArray{T}, grad_y::AbstractArray{T}) where T
    m = similar(grad_x, eltype(T))
    p = similar(m)
    for I in eachindex(grad_x, grad_y)
        m[I], p[I] = magnitude_phase(grad_x[I], grad_y[I])
    end
    m, p
end

# Return the magnitude and phase of the gradients in an image
function magnitude_phase(img::AbstractArray, method::Function=KernelFactors.ando3, border::AbstractString="replicate")
    grad_x, grad_y = imgradients(img, method, border)
    return magnitude_phase(grad_x, grad_y)
end

# Return the x-y gradients and magnitude and phase of gradients in an image
"""
```
grad_y, grad_x, mag, orient = imedge(img, kernelfun=KernelFactors.ando3, border="replicate")
```

Edge-detection filtering. `kernelfun` is a valid kernel function for
[`imgradients`](@ref), defaulting to [`KernelFactors.ando3`](@ref).
`border` is any of the boundary conditions specified in `padarray`.

Returns a tuple `(grad_y, grad_x, mag, orient)`, which are the horizontal
gradient, vertical gradient, and the magnitude and orientation of the strongest
edge, respectively.
"""
function imedge(img::AbstractArray, kernelfun=KernelFactors.ando3, border::AbstractString="replicate")
    grad_y, grad_x = imgradients(img, kernelfun, border)
    mag = magnitude(grad_y, grad_x)
    orient = orientation(grad_y, grad_x)
    return (grad_y, grad_x, mag, orient)
end

# Thin edges
"""
```
thinned = thin_edges(img, gradientangle, [border])
thinned, subpix = thin_edges_subpix(img, gradientangle, [border])
thinned, subpix = thin_edges_nonmaxsup(img, gradientangle, [border]; [radius::Float64=1.35], [theta=pi/180])
thinned, subpix = thin_edges_nonmaxsup_subpix(img, gradientangle, [border]; [radius::Float64=1.35], [theta=pi/180])
```

Edge thinning for 2D edge images.  Currently the only algorithm available is
non-maximal suppression, which takes an edge image and its gradient angle, and
checks each edge point for local maximality in the direction of the gradient.
The returned image is non-zero only at maximal edge locations.

`border` is any of the boundary conditions specified in `padarray`.

In addition to the maximal edge image, the `_subpix` versions of these functions
also return an estimate of the subpixel location of each local maxima, as a 2D
array or image of `Graphics.Point` objects.  Additionally, each local maxima is
adjusted to the estimated value at the subpixel location.

Currently, the `_nonmaxsup` functions are identical to the first two function
calls, except that they also accept additional keyword arguments.  `radius`
indicates the step size to use when searching in the direction of the gradient;
values between 1.2 and 1.5 are suggested (default 1.35).  `theta` indicates the
step size to use when discretizing angles in the `gradientangle` image, in
radians (default: 1 degree in radians = pi/180).

Example:

```
g = rgb2gray(rgb_image)
gx, gy = imgradients(g)
mag, grad_angle = magnitude_phase(gx,gy)
mag[mag .< 0.5] = 0.0  # Threshold magnitude image
thinned, subpix =  thin_edges_subpix(mag, grad_angle)
```
"""
thin_edges(img::AbstractArray{T,2}, gradientangles::AbstractArray, border::AbstractString="replicate") where {T} =
    thin_edges_nonmaxsup(img, gradientangles, border)
thin_edges_subpix(img::AbstractArray{T,2}, gradientangles::AbstractArray, border::AbstractString="replicate") where {T} =
    thin_edges_nonmaxsup_subpix(img, gradientangles, border)

# Code below is related to non-maximal suppression, and was ported to Julia from
# http://www.csse.uwa.edu.au/~pk/research/matlabfns/Spatial/nonmaxsup.m
# (Please conserve the original copyright below.)

# NONMAXSUP - Non-maxima suppression
#
# Usage:
#          (im,location) = nonmaxsup(img, gradientangles, radius);
#
# Function for performing non-maxima suppression on an image using
# gradient angles.  Gradient angles are assumed to be in radians.
#
# Input:
#   img - image to be non-maxima suppressed.
#
#   gradientangles - image containing gradient angles around each pixel in radians
#                    (-pi,pi)
#
#   radius  - Distance in pixel units to be looked at on each side of each
#             pixel when determining whether it is a local maxima or not.
#             This value cannot be less than 1.
#             (Suggested value about 1.2 - 1.5)
#
# Returns:
#   im        - Non maximally suppressed image.
#   location  - `Graphics.Point` image holding subpixel locations of edge
#               points.
#
# Notes:
#
# This function uses bilinear interpolation to estimate
# intensity values at ideal, real-valued pixel locations on each side of
# pixels to determine if they are local maxima.

# Copyright (c) 1996-2013 Peter Kovesi
# Centre for Exploration Targeting
# The University of Western Australia
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# The Software is provided "as is", without warranty of any kind.

# December  1996 - Original version
# September 2004 - Subpixel localization added
# August    2005 - Made Octave compatible
# October   2013 - Final thinning applied to binary image for Octave
#                  compatbility (Thanks to Chris Pudney)
# June      2014 - Ported (and modified significantly) to Julia (Kevin Squire)

import .Point

if !applicable(zero, Point)
    import Base.zero
    zero(::Type{Point}) = Point(0.0,0.0)
end

# Used to encode the sign, integral, and fractional components of
# an offset from a coordinate
struct CoordOffset
    s::Int      # sign
    i::Int      # integer part
    f::Float64  # fractional part
end

CoordOffset(x::Float64) = ((frac,i) = modf(x); CoordOffset(sign(frac), round(Int, i), abs(frac)))
(-)(off::CoordOffset) = CoordOffset(-off.s,-off.i, off.f)
(*)(x::Number, off::CoordOffset) = x*(off.i + off.s*off.f)
(*)(off::CoordOffset, x::Number) = x*(off.i + off.s*off.f)
(+)(x::Number, off::CoordOffset) = x + off.i + off.s*off.f
(+)(off::CoordOffset, x::Number) = x + off.i + off.s*off.f

# Precalculate x and y offsets relative to centre pixel for each orientation angle
function _calc_discrete_offsets(θ, radius)

    θ_count = round(Int, 2π/θ)
    θ = 2π/θ_count
    angles = (0:θ_count)*θ

    # x and y offset of points at specified radius and angles
    # from each reference position.

    xoffs = [CoordOffset( x) for x in  radius * cos.(angles)]
    yoffs = [CoordOffset(-y) for y in  radius * sin.(angles)]

    return θ, xoffs, yoffs
end

_discretize_angle(angle::AbstractFloat, invθ) =
    angle < 0 ? round(Int, (angle + 2π)*invθ)+1 : round(Int, angle*invθ)+1

# Interpolate the value of an offset from a particular pixel
#
# Returns (interpolated value, min_value of adjacent pixels in direction of offset)
#
# The second value is made available to eliminate double edges; if the value at
# (x,y) is less than the value or values adjacent to it in the direction of the
# gradient (xoff,yoff), then it is not a local maximum

function _interp_offset(img::AbstractArray, x::Integer, y::Integer, xoff::CoordOffset, yoff::CoordOffset, Ix, Iy, pad)
    fx = Ix[x + xoff.i + pad]
    fy = Iy[y + yoff.i + pad]
    cx = Ix[x + xoff.i + xoff.s + pad]
    cy = Iy[y + yoff.i + yoff.s + pad]

    tl = img[fy,fx]    # Value at bottom left integer pixel location.
    tr = img[fy,cx]    # bottom right
    bl = img[cy,fx]    # top left
    br = img[cy,cx]    # top right

    upperavg = tl + xoff.f * (tr - tl)  # Now use bilinear interpolation to
    loweravg = bl + xoff.f * (br - bl)  # estimate value at x,y

    min_adjacent = (fx == x) & (fy == y) ? min(tr,bl) : tl

    return (upperavg + yoff.f * (loweravg - upperavg), min_adjacent)
end

# Core edge thinning algorithm using nonmaximal suppression
function thin_edges_nonmaxsup_core!(out::AbstractArray{T,2}, location::AbstractArray{Point,2},
                                    img::AbstractArray{T,2}, gradientangles::AbstractMatrix, radius, border, theta) where T
    calc_subpixel = !isempty(location)

    # Error checking
    size(img) == size(gradientangles) == size(out) || error("image, gradient angle, and output image must all be the same size")
    calc_subpixel && size(location) != size(img) && error("subpixel location has a different size than the input image")
    radius < 1.0 && error("radius must be >= 1")

    # Precalculate x and y offsets relative to centre pixel for each orientation angle
    θ, xoffs, yoffs = _calc_discrete_offsets(theta, radius)
    iθ = 1/θ

    # Indexes to use for border handling
    pad = ceil(Int, radius)
    Ix = Images.padindexes(img, 2, pad, pad, border)
    Iy = Images.padindexes(img, 1, pad, pad, border)

    # Now run through the image interpolating grey values on each side
    # of the centre pixel to be used for the non-maximal suppression.

    (height,width) = size(img)

    for x = 1:width, y = 1:height
        (c = img[y,x]) == 0 && continue  # For thresholded images

        or = _discretize_angle(gradientangles[y,x],iθ)   # Disretized orientation
        v1, n1 = _interp_offset(img, x, y, xoffs[or], yoffs[or], Ix, Iy, pad)

        if (c > v1) & (c >= n1) # We need to check the value on the other side...
            v2, n2 = _interp_offset(img, x, y, -xoffs[or], -yoffs[or], Ix, Iy, pad)

            if (c > v2) & (c >= n2)  # This is a local maximum.
                                     # Record value in the output image.
                if calc_subpixel
                    # Solve for coefficients of parabola that passes through
                    # [-1, v2]  [0, img] and [1, v1].
                    # v = a*r^2 + b*r + c

                    # c = img[y,x]
                    a = (v1 + v2)/2 - c
                    b = a + c - v2

                    # location where maxima of fitted parabola occurs
                    r = -b/2a
                    location[y,x] = Point(x + r*xoffs[or], y + r*yoffs[or])

                    if T<:AbstractFloat
                        # Store the interpolated value
                        out[y,x] = a*r^2 + b*r + c
                    else
                        out[y,x] = c
                    end
                else
                    out[y,x] = c
                end
            end
        end
    end

    out
end

# Main function call when subpixel location of edges is not needed
function thin_edges_nonmaxsup!(out, img, gradientangles, border::AbstractString="replicate";
                               radius::Float64=1.35, theta=pi/180)
    thin_edges_nonmaxsup_core!(out, Matrix{Point}(0,0), img, gradientangles, radius, border, theta)
end

function thin_edges_nonmaxsup(img, gradientangles, border::AbstractString="replicate";
                                 radius::Float64=1.35, theta=pi/180)
    (height,width) = size(img)
    out = zeros(eltype(img), height, width)
    thin_edges_nonmaxsup_core!(out, Matrix{Point}(undef,0,0), img, gradientangles, radius, border, theta)
end

# Main function call when subpixel location of edges is desired
function thin_edges_nonmaxsup_subpix!(out, location, img, gradientangles,
                                      border::AbstractString="replicate";
                                      radius::Float64=1.35, theta=pi/180)
    eltype(location) != Point && error("Preallocated subpixel location array/image must have element type Graphics.Point")

    thin_edges_nonmaxsup_core!(out, location, img, gradientangles, radius, border, theta)
    img, location
end

function thin_edges_nonmaxsup_subpix(img, gradientangles,
                                     border::AbstractString="replicate";
                                     radius::Float64=1.35, theta=pi/180)
    (height,width) = size(img)
    out = zeros(eltype(img), height, width)
    location = zeros(Point, height, width)
    thin_edges_nonmaxsup_core!(out, location, img, gradientangles, radius, border, theta)
    out, location
end

"""
```
canny_edges = canny(img, (upper, lower), sigma=1.4)
```

Performs Canny Edge Detection on the input image.

Parameters :

  (upper, lower) :  Bounds for hysteresis thresholding
  sigma :           Specifies the standard deviation of the gaussian filter

# Example

```julia
imgedg = canny(img, (Percentile(80), Percentile(20)))
```
"""
function canny(img_gray::AbstractMatrix{T}, threshold::Tuple{N,N}, sigma::Number = 1.4) where {T<:NumberLike, N<:Union{NumberLike,Percentile{NumberLike}}}
    img_grayf = imfilter(img_gray, KernelFactors.IIRGaussian((sigma,sigma)), NA())
    img_grad_y, img_grad_x = imgradients(img_grayf, KernelFactors.sobel)
    img_mag, img_phase = magnitude_phase(img_grad_x, img_grad_y)
    img_nonMaxSup = thin_edges_nonmaxsup(img_mag, img_phase)
    if N<:Percentile{}
        upperThreshold ,lowerThreshold = StatsBase.percentile(img_nonMaxSup[:], [threshold[i].p for i=1:2])
    else
        upperThreshold, lowerThreshold = threshold
    end
    img_thresholded = hysteresis_thresholding(img_nonMaxSup, upperThreshold, lowerThreshold)
    edges = map(i -> i >= 0.9, img_thresholded)
    edges
end

canny(img::AbstractMatrix, threshold::Tuple{N,N}, args...) where {N<:Union{NumberLike,Percentile{NumberLike}}} =
    canny(convert(Array{Gray}, img), args...)

function hysteresis_thresholding(img_nonMaxSup::AbstractArray{T, 2}, upperThreshold::Number, lowerThreshold::Number) where T
    img_thresholded = map(i -> i > lowerThreshold ? i > upperThreshold ? 1.0 : 0.5 : 0.0, img_nonMaxSup)
    queue = CartesianIndex{2}[]
    R = CartesianIndices(size(img_thresholded))

    I1, Iend = first(R), last(R)
    for I in R
      if img_thresholded[I] == 1.0
        img_thresholded[I] = 0.9
        push!(queue, I)
        while !isempty(queue)
          q_top = popfirst!(queue)
          for J in CartesianIndices(map((f,l)->f:l,(max(I1, q_top - I1)).I, (min(Iend, q_top + I1)).I))
            if img_thresholded[J] == 1.0 || img_thresholded[J] == 0.5
              img_thresholded[J] = 0.9
              push!(queue, J)
            end
          end
        end
      end
    end
    img_thresholded
end

function padindexes(img::AbstractArray{T,n}, dim, prepad, postpad, border::AbstractString) where {T,n}
    M = size(img, dim)
    I = Vector{Int}(undef, M + prepad + postpad)
    I = [(1 - prepad):(M + postpad);]
    if border == "replicate"
        I = min.(max.(I, 1), M)
    elseif border == "circular"
        I = 1 .+ mod.(I .- 1, M)
    elseif border == "symmetric"
        I = [1:M; M:-1:1][1 .+ mod.(I .- 1, 2 * M)]
    elseif border == "reflect"
        I = [1:M; M-1:-1:2][1 .+ mod.(I .- 1, 2 * M - 2)]
    else
        error("unknown border condition")
    end
    I
end

ext: jl


text:
# calling a function in a shared library:
lang = ccall( (:getenv, "libc"),Ptr{Uint8}, (Ptr{Uint8},), "LANGUAGE")
# ok on Linux
# on Windows: ERROR: error compiling anonymous: could not load module libc:
# The specified module could not be found.
bytestring(lang) #> "en_US"

# test existence of library:
find_library(["libc"]) #> "libc"

max_stars_repo_path: Module 1/Chapter09/callc.jl

max_issues_repo_path: Module 1/Chapter09/callc.jl

max_forks_repo_path: Module 1/Chapter09/callc.jl


text:
export ComputeMode, ComputeFast, ComputeAccurate

"""
    ComputeMode
"""
struct ComputeMode{T} end

"""
    ComputeFast

Shorthand for `ComputeMode{:fast}`
"""
const ComputeFast = ComputeMode{:fast}

"""
    ComputeAccurate

Shorthand for `ComputeMode{:ComputeAccurate}`
"""
const ComputeAccurate = ComputeMode{:ComputeAccurate}
Base.promote_rule(::Type{ComputeFast}, ::Type{ComputeAccurate}) = ComputeAccurate


text:
##################################################################
#
# Pointers to Sundials objects
#
##################################################################

"""
    Base type for dummy placeholders that help to
    providing typed pointers for Sundials objects
    (KINSOL, CVODE, IDA).

    See `Handle`.
"""
abstract type AbstractSundialsObject end

struct CVODEMem <: AbstractSundialsObject end
const CVODEMemPtr = Ptr{CVODEMem}

struct ARKODEMem <: AbstractSundialsObject end
const ARKODEMemPtr = Ptr{ARKODEMem}

struct IDAMem <: AbstractSundialsObject end
const IDAMemPtr = Ptr{IDAMem}

struct KINMem <: AbstractSundialsObject end
const KINMemPtr = Ptr{KINMem}

abstract type SundialsHandle end

"""
   Handle for Sundials objects (CVODE, IDA, KIN).

   Wraps the reference to the pointer to the Sundials object.
   Manages automatic destruction of the referenced objects when it is
   no longer in use.
"""
struct Handle{T <: AbstractSundialsObject} <: SundialsHandle
    ptr_ref::Ref{Ptr{T}} # pointer to a pointer

    function Handle(ptr::Ptr{T}) where T <: AbstractSundialsObject
        h = new{T}(Ref{Ptr{T}}(ptr))
        finalizer(release_handle, h.ptr_ref)
        return h
    end
end

mutable struct MatrixHandle{T<:SundialsMatrix} <: SundialsHandle
    ptr::SUNMatrix
    destroyed::Bool
    function MatrixHandle(ptr::SUNMatrix,M::T) where T<:SundialsMatrix
        h = new{T}(ptr,false)
        finalizer(release_handle, h)
        return h
    end
end

mutable struct LinSolHandle{T<:SundialsLinSol} <: SundialsHandle
    ptr::SUNLinearSolver
    destroyed::Bool
    function LinSolHandle(ptr::SUNLinearSolver,M::T) where T<:SundialsLinSol
        h = new{T}(ptr,false)
        finalizer(release_handle, h)
        return h
    end
end

Base.unsafe_convert(::Type{Ptr{T}}, h::Handle{T}) where T = h.ptr_ref[]
Base.convert(::Type{Ptr{T}}, h::Handle{T}) where T = h.ptr_ref[]
Base.convert(::Type{Ptr{Ptr{T}}}, h::Handle{T}) where {T} = convert(Ptr{Ptr{T}}, h.ptr_ref[])

release_handle(ptr_ref::Ref{Ptr{T}}) where {T} = throw(MethodError("Freeing objects of type $T not supported"))
release_handle(ptr_ref::Ref{Ptr{KINMem}}) = ((ptr_ref[] != C_NULL) && KINFree(ptr_ref); nothing)
release_handle(ptr_ref::Ref{Ptr{CVODEMem}}) = ((ptr_ref[] != C_NULL) && CVodeFree(ptr_ref); nothing)
release_handle(ptr_ref::Ref{Ptr{ARKODEMem}}) = ((ptr_ref[] != C_NULL) && ARKodeFree(ptr_ref); nothing)
release_handle(ptr_ref::Ref{Ptr{IDAMem}}) = ((ptr_ref[] != C_NULL) && IDAFree(ptr_ref); nothing)

function release_handle(h::MatrixHandle{DenseMatrix})
    if !isempty(h)
        Sundials.SUNMatDestroy_Dense(h.ptr)
        h.destroyed = true
    end
    nothing
end
function release_handle(h::MatrixHandle{BandMatrix})
    if !isempty(h)
        Sundials.SUNMatDestroy_Band(h.ptr)
        h.destroyed = true
    end
    nothing
end
function release_handle(h::MatrixHandle{SparseMatrix})
    if !isempty(h)
        Sundials.SUNMatDestroy_Sparse(h.ptr)
        h.destroyed = true
    end
    nothing
end

function release_handle(h::LinSolHandle{Dense})
    if !isempty(h)
        Sundials.SUNLinSolFree_Dense(h.ptr)
        h.destroyed = true
    end
    nothing
end

function release_handle(h::LinSolHandle{Band})
    if !isempty(h)
        Sundials.SUNLinSolFree_Band(h.ptr)
        h.destroyed = true
    end
    nothing
end

function release_handle(h::LinSolHandle{SPGMR})
    if !isempty(h)
        Sundials.SUNLinSolFree_SPGMR(h.ptr)
        h.destroyed = true
    end
    nothing
end

function release_handle(h::LinSolHandle{SPFGMR})
    if !isempty(h)
        Sundials.SUNLinSolFree_SPFGMR(h.ptr)
        h.destroyed = true
    end
    nothing
end

function release_handle(h::LinSolHandle{SPBCGS})
    if !isempty(h)
        Sundials.SUNLinSolFree_SPBCGS(h.ptr)
        h.destroyed = true
    end
    nothing
end

function release_handle(h::LinSolHandle{PCG})
    if !isempty(h)
        Sundials.SUNLinSolFree_PCG(h.ptr)
        h.destroyed = true
    end
    nothing
end

function release_handle(h::LinSolHandle{PTFQMR})
    if !isempty(h)
        Sundials.SUNLinSolFree_SPTFQMR(h.ptr)
        h.destroyed = true
    end
    nothing
end

function release_handle(h::LinSolHandle{KLU})
    if !isempty(h)
        Sundials.SUNLinSolFree_KLU(h.ptr)
        h.destroyed = true
    end
    nothing
end

Base.empty!(h::LinSolHandle) = release_handle(h)
Base.empty!(h::MatrixHandle) = release_handle(h)
Base.empty!(h::Handle{T}) where {T} = release_handle(h.ptr_ref)
Base.isempty(h::Handle) = h.ptr_ref[] == C_NULL
Base.isempty(h::MatrixHandle) = h.destroyed
Base.isempty(h::LinSolHandle) = h.destroyed

##################################################################
#
# Convenience typealiases for Sundials handles
#
##################################################################

const CVODEh  =  Handle{CVODEMem}
const ARKODEh =  Handle{ARKODEMem}
const KINh    =  Handle{KINMem}
const IDAh    =  Handle{IDAMem}


text:
import CpuId, LLVM

let features = filter(ext -> (m = match(r"\d", ext); isnothing(m) ? true : m.offset != 2 ) , split(unsafe_string(LLVM.API.LLVMGetHostCPUFeatures()), ','))
    offsetnottwo(::Nothing) = true
    offsetnottwo(m::RegexMatch) = m.offset != 2

    avx512f = any(isequal("+avx512f"), features)
    avx2 = any(isequal("+avx2"), features)
    avx = any(isequal("+avx"), features)

    register_size = avx512f ? 64 : (avx ? 32 : 16)
    register_count = avx512f ? 32 : 16
    cache_size = CpuId.cachesize()
    num_cores = CpuId.cpucores()

    @eval const REGISTER_SIZE = $register_size
    @eval const REGISTER_COUNT = $register_count
    @eval const FP256 = $(CpuId.cpufeature(CpuId.FP256)) # Is AVX2 fast?
    @eval const CACHELINE_SIZE = $(CpuId.cachelinesize())
    @eval const CACHE_SIZE = $cache_size
    @eval const NUM_CORES = $num_cores
    @eval const SIMD_NATIVE_INTEGERS = $(avx2)

    for ext ∈ features
        @eval const $(Symbol(replace(Base.Unicode.uppercase(ext[2:end]), r"\." => "_"))) = $(first(ext) == '+')
    end
end

const FMA3 = FMA


text:
using Test
using CompilerPluginTools
using MLStyle
using Expronicon

function foo(x)
    if x > 0
        2x
    else
        3x
    end
end

@testset "map(f, ::CodeInfo)" begin
    ci = code_lowered(foo, (Float64, ))[1]
    test_ci = map(NewCodeInfo(ci)) do stmt
        @match stmt begin
            Expr(:call, f, args...) => Expr(:call, :overdub, f, args...)
            _ => stmt
        end
    end

    @test_codeinfo test_ci begin
        :(overdub($(GlobalRef(Main, :>)), $(SlotNumber(2)), 0))
        ::GotoIfNot
        :(overdub($(GlobalRef(Main, :*)), 2, $(SlotNumber(2))))
        ReturnNode(SSAValue(3))
        :(overdub($(GlobalRef(Main, :*)), 3, $(SlotNumber(2))))
        ReturnNode(SSAValue(5))
    end
end

@testset "NewCodeInfo" begin
    ci = code_lowered(foo, (Float64, ))[1]
    new = NewCodeInfo(ci)
    for (v, stmt) in new
        @switch stmt begin
            @case Expr(:call, GlobalRef(Main, :>), a, b)
                delete!(new, v)
                insert!(new, v, :(1 + 1))
            @case Expr(:call, f, args...)
                x = insert!(new, v, "variable: %$v")
                new[v] = Expr(:call, :overdub, x, f, args...)
            @case _
                nothing
        end
    end
    test_ci = finish(new)

    @test_codeinfo test_ci begin
        :(1 + 1)
        ::GotoIfNot

        "variable: %3"

        :(overdub($(SSAValue(3)), $(GlobalRef(Main, :*)), 2, $(SlotNumber(2))))
        ReturnNode(SSAValue(4))

        "variable: %5"
    end

    @testset "slot insertion" begin
        ci = code_lowered(foo, (Float64, ))[1]
        new = NewCodeInfo(ci)
        insert!(new.slots, 1, Symbol("#a#"))
        insert!(new.slots, 1, Symbol("#b#"))
        nci = finish(new)
        @test nci.code[1] == Expr(:call, GlobalRef(Main, :>), SlotNumber(4), 0)
        @test nci.code[2] == GotoIfNot(SSAValue(1), 5)
        @test nci.slotnames == [Symbol("#b#"), Symbol("#a#"), Symbol("#self#"), :x]
    end

    @testset "pc=1 multi push" begin
        ci = code_lowered(foo, (Float64, ))[1]
        new = NewCodeInfo(ci)
        push!(new, :(1 + 1))
        push!(new, :(1 + 2))
        push!(new, :(1 + 3))

        for (v, stmt) in new
            if v == 3
                x = insert!(new, v, :(1 + 2))
                new[v] = :(1 + $x)
            end
        end

        test_ci = finish(new)
        @test test_ci.code[1] == :(1 + 1)
        @test test_ci.code[2] == :(1 + 2)
        @test test_ci.code[3] == :(1 + 3)
        @test test_ci.code[4] == Expr(:call, GlobalRef(Main, :>), SlotNumber(2), 0)
        @test test_ci.code[5] == GotoIfNot(SSAValue(4), 9)
        @test test_ci.code[6] == :(1 + 2)
        @test test_ci.code[7] == :(1 + $(SSAValue(6)))
    end

    @testset "multiple insert!(new, 2, ...)" begin
        ci = code_lowered(foo, (Float64, ))[1]
        new = NewCodeInfo(ci)
        insert!(new, 2, :(1 + 1))
        insert!(new, 2, :(1 + 2))
        insert!(new, 2, :(1 + 3))
        test_ci = finish(new)

        @test test_ci.code[4] == :(1 + 1)
        @test test_ci.code[3] == :(1 + 2)
        @test test_ci.code[2] == :(1 + 3)
    end

    @testset "multiple push!(new, ...) with setindex" begin
        ci = code_lowered(foo, (Float64, ))[1]
        new = NewCodeInfo(ci)
        push!(new, :(1 + 1))
        push!(new, :(1 + 2))
        push!(new, :(1 + 3))

        for (v, stmt) in new
            if v == 3
                x = push!(new, :(1 + 2))
                x = push!(new, :(1 + 3))
                x = push!(new, :(1 + $x))
                new[v] = :(1 + $x)
            end
        end

        test_ci = finish(new)
        @test test_ci.code[9] == :(1 + $(SSAValue(8)))
    end
end

@testset "rm_code_coverage_effect" begin
    ci = @make_codeinfo begin
        #=%1 =# Expr(:code_coverage_effect)::Nothing
        #=%2 =# QuoteNode(1.0)::Float64
        #=%3 =# Expr(:call, sin, SSAValue(2))::Float64
        #=%4 =# Expr(:code_coverage_effect)::Nothing
        #=%5 =# QuoteNode(2.0)::Float64
        #=%6 =# Expr(:call, sin, SSAValue(5))::Float64
        #=%7 =# Expr(:code_coverage_effect)::Nothing
        #=%8 =# Expr(:call, <, SSAValue(6), QuoteNode(1.0))::Bool
        #=%9 =# GotoIfNot(SSAValue(8), 10)
        #=%10=# Expr(:call, +, SSAValue(3), SSAValue(6))::Float64
        #=%11=# ReturnNode(SSAValue(10))::Float64
    end

    test_ci = CompilerPluginTools.rm_code_coverage_effect(ci)

    @test_codeinfo test_ci begin
        #=%1=# QuoteNode(1.0)::Float64
        #=%2=# Expr(:call, sin, SSAValue(1))::Float64
        #=%3=# QuoteNode(2.0)::Float64
        #=%4=# Expr(:call, sin, SSAValue(3))::Float64
        #=%5=# Expr(:call, <, SSAValue(4), QuoteNode(1.0))::Bool
        #=%6=# GotoIfNot(SSAValue(5), 7)
        #=%7=# Expr(:call, +, SSAValue(2), SSAValue(4))::Float64
        #=%8=# ReturnNode(SSAValue(7))::Float64
    end
end

@testset "test code_coverage_effect bb start" begin
    ci = @make_codeinfo begin
        Expr(:code_coverage_effect)::Nothing
        Expr(:code_coverage_effect)::Nothing
        Expr(:code_coverage_effect)::Nothing
        Expr(:code_coverage_effect)::Nothing
        Expr(:call, GlobalRef(Main, :foo), 2)::Nothing
        Expr(:call, GlobalRef(Main, :measure), 2)::Int
        Expr(:code_coverage_effect)::Nothing
        Expr(:code_coverage_effect)::Nothing
        Expr(:call, GlobalRef(Main, :measure_cmp), SSAValue(6), 1)::Bool
        GotoIfNot(SSAValue(9), 14)
        Expr(:code_coverage_effect)::Nothing
        Expr(:code_coverage_effect)::Nothing
        Expr(:call, GlobalRef(Main, :foo), 2)::Nothing
        Expr(:code_coverage_effect)::Nothing
        ReturnNode(nothing)::Nothing
    end

    stmt = CompilerPluginTools.rm_code_coverage_effect(ci).code[4]
    @test stmt.cond == SSAValue(3)
    @test stmt.dest == 6
end


text:
module ParametersEstimator

using NamedArrays
using LinearAlgebra, ForwardDiff, FiniteDiff
using Statistics
using DataFormatter
using DataFrames
using GLMakie

include("functions.jl")

end

max_stars_repo_name: Leooop/ParametersEstimator

max_issues_repo_name: Leooop/ParametersEstimator

max_forks_repo_name: Leooop/ParametersEstimator


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#     Define constants.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

################################################################################
#                          Types for C compatibility
################################################################################

const __u8  = UInt8
const __u16 = UInt16
const __u32 = UInt32
const __u64 = UInt64

################################################################################
#                                     GPIO
################################################################################

const GPIOMEM_SIZE = 1024

const GPIO_MODE_GET = Dict(0b000 => :in,
                           0b001 => :out,
                           0b100 => :alt0,
                           0b101 => :alt1,
                           0b110 => :alt2,
                           0b111 => :alt3,
                           0b011 => :alt4,
                           0b010 => :alt5)

const GPIO_MODE_SET = Dict(:in   => 0b000,
                           :out  => 0b001,
                           :alt0 => 0b100,
                           :alt1 => 0b101,
                           :alt2 => 0b110,
                           :alt3 => 0b111,
                           :alt4 => 0b011,
                           :alt5 => 0b010)

################################################################################
#                                     SPI
################################################################################

const _SPI_BUFFER_SIZE = 16


text:
struct CatType{T} end

function make_url(::Type{CatType{:_cat}}, info::Esinfo, paths::AbstractString, name::AbstractString)
	"$(info.url)/_cat/$paths/$name"
end

function make_url(::Type{CatType{:_cat}}, info::Esinfo,paths::AbstractString)
	"$(info.url)/_cat/$paths"
end

macro catgenfun(methods, paths, names )
	funname = esc(Symbol(string("cat_", paths )) )
	name    = esc(Symbol(names) )

	Expr(:function ,
		Expr(:call ,
			funname ,
			Expr(:parameters,
				Expr(:(...), :kw)),
			Expr(:(::) , :info , :Esinfo),
			Expr(:(::) , name , :AbstractString)
			) ,
		Expr(:block ,
			quote
				querys = Dict(kw...)
				url    = make_url(CatType{:_cat}, info, $paths, name)
				@catexport info $methods  url  querys
			end )
	)
end

macro catgenfun(methods, paths )
	funname = esc(Symbol(string("cat_", paths )) )

	Expr(:function ,
		Expr(:call ,
			funname ,
			Expr(:parameters,
				Expr(:(...), :kw)),
			Expr(:(::) , :info , :Esinfo)
			) ,
		Expr(:block ,
			quote
				querys = Dict(kw...)
				url    = make_url(CatType{:_cat}, info, $paths)
				@catexport info $methods  url  querys
			end )
	)
end

##
@catgenfun "GET" "aliases"  "name"
@catgenfun "GET" "aliases"
@catgenfun "GET" "fielddata" "fields"
@catgenfun "GET" "fielddata"
@catgenfun "GET" "allocation"  "node_id"
@catgenfun "GET" "allocation"
@catgenfun "GET" "count"  "index"
@catgenfun "GET" "count"
@catgenfun "GET" "health"
@catgenfun "GET" "master"
@catgenfun "GET" "nodeattrs"
@catgenfun "GET" "pending_tasks"
@catgenfun "GET" "plugins"
@catgenfun "GET" "recovery"  "index"
@catgenfun "GET" "recovery"
@catgenfun "GET" "repositories"
@catgenfun "GET" "segments"  "index"
@catgenfun "GET" "segments"
@catgenfun "GET" "shards"  "index"
@catgenfun "GET" "shards"
@catgenfun "GET" "snapshots"  "repository"
@catgenfun "GET" "snapshots"
@catgenfun "GET" "tasks"
@catgenfun "GET" "templates"  "name"
@catgenfun "GET" "templates"
@catgenfun "GET" "thread_pool"  "thread_pool_patterns"
@catgenfun "GET" "thread_pool"

max_stars_repo_name: frankill/ES

max_issues_repo_name: frankill/ES

max_forks_repo_name: frankill/ES


text:
module FESpaceWithConstantFixedTests
using Gridap
using Gridap.Geometry
using Gridap.FESpaces
using Gridap.ReferenceFEs
using Test

domain = (0,1,0,1)
partition = (4,4)
model = CartesianDiscreteModel(domain,partition)

trian = get_triangulation(model)

order = 2

V = FESpace(model,ReferenceFE(lagrangian,Float64,order);conformity=:L2)

V0 = FESpaceWithConstantFixed(V,true,rand(1:num_free_dofs(V)))
test_single_field_fe_space(V0)

@test Gridap.FESpaces.ConstantApproach(V0) == Gridap.FESpaces.FixConstant()

uh0 = interpolate(V0) do x
    sin(4*pi*(x[1]+x[2]^2)) + 3
end
using Gridap.Visualization

#writevtk(trian,"trian",nsubcells=20,cellfields=["uh0"=>uh0])

end # module


text:
using BinaryBuilder

# Collection of sources required to build Python
sources = [
    "https://www.python.org/ftp/python/3.6.5/Python-3.6.5.tar.xz" =>
    "f434053ba1b5c8a5cc597e966ead3c5143012af827fd3f0697d21450bb8d87a6"
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/Python-3.6.5/
make distclean
./configure --prefix=${prefix} --host=${target} --build=x86_64-linux-gnu
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Windows(:i686),
    Windows(:x86_64),
    Linux(:i686, :glibc),
    Linux(:x86_64, :glibc),
    Linux(:aarch64, :glibc),
    Linux(:armv7l, :glibc),
    Linux(:powerpc64le, :glibc),
    MacOS()
]

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libpython", :libpython),
    ExecutableProduct(prefix, "python", :python)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    "https://github.com/staticfloat/ZlibBuilder/releases/download/v1.2.11-3/build.jl",
    "https://github.com/staticfloat/Bzip2Builder/releases/download/v1.0.6-0/build.jl",
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, "Python", sources, script, platforms, products, dependencies)

max_stars_repo_name: staticfloat/PythonBuilder

max_issues_repo_name: staticfloat/PythonBuilder

max_forks_repo_name: staticfloat/PythonBuilder


text:
function generate_2D_dataset(n::Int, x1_scope::Tuple, x2_scope::Tuple, accurancy::Number, equidistant::Bool=false)
# n: Number of generated points in 2D (i.e. points in the form [x1, x2])
# x1_scope: Definition domain of the 2D function if x=x1,y=x2.
# x2_scope: Function range of the 2D function if x=x1,y=x2.
# accurancy: Decimal number representing valid decimal places in function values.
# equidistant: Whether x1 points should be equidistant or not(defaut), i.e. generated randomly as x2 points.
    x1_first, x1_last = x1_scope
    x1_min, x1_max = x1_scope
    x2_min, x2_max = x2_scope

    if equidistant == true
        step = x1_scope/n
        x1_set = collect(x1_min : step : x1_max)
    else
        x1_set = rand(x1_min : accurancy : x1_max, n)
    end
    x2_set = rand(x2_min : accurancy : x2_max, n)
    return [x1_set x2_set]
end

function create_2D_linear_function(omega,b)
    # returns anonymous function which accepts single number 'x1' as input
    # a*x + b*y +c = 0 --> y = -(c+ax)/b
    # omega1*x1 + omega2*x2 + b = 0 --> x2 = -(b+w1*x1)/w2
    w1,w2 = omega[1],omega[2]
    f(x1) = -(b+w1*x1)/w2
    return f
end

function get_minmax_vals(dataset::Array, two_dimensions=true)
    if two_dimensions != true
        println(stderr, "Only 2D supported, exiting...\n")
        exit(1)
    end
    x1 = dataset[:,1]
    x2 = dataset[:,2]
    x1_bounds=[minimum(x1),maximum(x1)]
    x2_bounds=[minimum(x2),maximum(x2)]
    return x1_bounds, x2_bounds
end

max_stars_repo_name: Dano-drevo/sop-julia-demo

max_issues_repo_name: Dano-drevo/sop-julia-demo

max_forks_repo_name: Dano-drevo/sop-julia-demo


text:
function runabcCancer(ABCsetup::ABCRejectionModel, targetdata; progress = false)

  ABCsetup.nmodels > 1 || error("Only 1 model specified, use ABCRejection method to estimate parameters for a single model")

  #initalize array of particles
  particles = Array{ApproxBayes.ParticleRejectionModel}(undef, ABCsetup.Models[1].nparticles)

  i = 1 #set particle indicator to 1
  its = 0 #keep track of number of iterations
  distvec = zeros(Float64, ABCsetup.Models[1].nparticles) #store distances in an array

  if progress == true
    p = Progress(ABCsetup.Models[1].nparticles, 1, "Running ABC rejection algorithm...", 30)
  end

  dist, out, newparams = 0.0,0.0,0.0

  while (i < (ABCsetup.Models[1].nparticles + 1)) & (its < ABCsetup.Models[1].maxiterations)

    its += 1
    #sample uniformly from models
    model = rand(1:ABCsetup.nmodels)
    correctmodel = false
    while correctmodel == false
      #get new proposal parameters
      newparams = ApproxBayes.getproposal(ABCsetup.Models[model].prior, ABCsetup.Models[model].nparams)
      #simulate with new parameters
      dist, out, cm = ABCsetup.Models[model].simfunc(newparams, ABCsetup.Models[model].constants, targetdata)
      correctmodel = cm
    end

    #if simulated data is less than target tolerance accept particle
    if dist < ABCsetup.Models[1].ϵ
      particles[i] = ApproxBayes.ParticleRejectionModel(newparams, model, dist, out)
      distvec[i] = dist
      i +=1
      if progress == true
        next!(p)
      end
    end
  end

  i > ABCsetup.Models[1].nparticles || error("Only accepted $(i-1) particles with ϵ < $(ABCsetup.Models[1].ϵ). \n\tDecrease ϵ or increase maxiterations ")

  out = ApproxBayes.ABCrejectionmodelresults(particles, its, ABCsetup, distvec)
  return out
end

function runabcCancer(ABCsetup::ABCSMCModel, targetdata; verbose = false, progress = false, savepopulations = false, resultsdirectory = "", Nmaxinf = 10^10, sname = "cancerabc", VAF = [])

  ABCsetup.nmodels > 1 || error("Only 1 model specified, use ABCSMC method to estimate parameters for a single model")

  #run first population with parameters sampled from prior
  if verbose
    println("##################################################")
    println("Use ABC rejection to get first population")
  end
  ABCrejresults = runabcCancer(ApproxBayes.ABCRejectionModel(
            map(x -> x.simfunc, ABCsetup.Models),
            map(x -> x.nparams, ABCsetup.Models),
            ABCsetup.Models[1].ϵ1,
            map(x -> x.prior, ABCsetup.Models),
            constants = map(x -> x.constants, ABCsetup.Models),
            nparticles = ABCsetup.Models[1].nparticles,
            maxiterations = ABCsetup.Models[1].maxiterations),
            targetdata, progress = progress);

  oldparticles, weights = ApproxBayes.setupSMCparticles(ABCrejresults, ABCsetup)
  ϵ = quantile(ABCrejresults.dist, ABCsetup.α) # set new ϵ to αth quantile
  ϵvec = [ϵ] #store epsilon values
  numsims = [ABCrejresults.numsims] #keep track of number of simualtions
  particles = Array{ApproxBayes.ParticleSMCModel}(undef, ABCsetup.nparticles) #define particles array
  weights, modelprob = ApproxBayes.getparticleweights(oldparticles, ABCsetup)
  ABCsetup = ApproxBayes.modelselection_kernel(ABCsetup, oldparticles)

  modelprob = ABCrejresults.modelfreq

  if verbose
    println("Run ABC SMC \n")
  end

  popnum = 1
  finalpop = false

  if verbose
    show(ApproxBayes.ABCSMCmodelresults(oldparticles, numsims, ABCsetup, ϵvec))
  end

  newparticle, dist, out, priorp = 0.0,0.0,0.0,0.0

  while (ϵ >= ABCsetup.ϵT) & (sum(numsims) <= ABCsetup.maxiterations)
    if verbose
      println("######################################## \n")
      println("########################################")
      println("Population number: $(popnum) \n")
    end

    i = 1 #set particle indicator to 1
    particles = Array{ApproxBayes.ParticleSMCModel}(undef, ABCsetup.nparticles)
    distvec = zeros(Float64, ABCsetup.nparticles)
    its = 1

    if progress == true
      p = Progress(ABCsetup.nparticles, 1, "ABC SMC population $(popnum), new ϵ: $(round(ϵ, digits = 2))...", 30)
    end
    while i < ABCsetup.nparticles + 1

      #draw model from previous model probabilities
      mstar = wsample(1:ABCsetup.nmodels, modelprob)

      #perturb model
      mdoublestar = ApproxBayes.perturbmodel(ABCsetup, mstar, modelprob)
      correctmodel = false
      #sometimes even though the input is for multiple clones,
      #simulations will be returned that don't have the requested number of clones,
      #this can happen if for example the simulation finishes before t1 is reached.
      #To overcome this we continually resample until we get simulations with the correct number of clones
      while correctmodel == false
        # sample particle with correct model
        j = wsample(1:ABCsetup.nparticles, weights[mdoublestar, :])
        particletemp = oldparticles[j]
        #perturb particle
        newparticle = ApproxBayes.perturbparticle(particletemp, ABCsetup.Models[mdoublestar].kernel)
        #calculate priorprob
        priorp = ApproxBayes.priorprob(newparticle.params, ABCsetup.Models[mdoublestar].prior)

        if priorp == 0.0 #return to beginning of loop if prior probability is 0
          break
        end

        #simulate with new parameters
        dist, out, cm = ABCsetup.Models[mdoublestar].simfunc(newparticle.params, ABCsetup.Models[mdoublestar].constants, targetdata)
        correctmodel = cm
      end

      if priorp == 0.0 #return to beginning of loop if prior probability is 0
        continue
      end

      #if simulated data is less than target tolerance accept particle
      if dist < ϵ
        particles[i] = newparticle
        particles[i].other = out
        particles[i].distance = dist
        distvec[i] = dist
        i += 1
        if progress == true
          next!(p)
        end
      end
      its += 1
    end

    particles, weights = ApproxBayes.smcweightsmodel(particles, oldparticles, ABCsetup, modelprob)
    weights, modelprob = ApproxBayes.getparticleweights(particles, ABCsetup)
    oldparticles = particles
    ABCsetup = ApproxBayes.modelselection_kernel(ABCsetup, particles)

    if finalpop == true
      break
    end

    ϵ = quantile(distvec, ABCsetup.α)

    if ϵ < ABCsetup.ϵT
      ϵ = ABCsetup.ϵT
      push!(ϵvec, ϵ)
      push!(numsims, its)
      popnum = popnum + 1
      finalpop = true
      continue
    end

    push!(ϵvec, ϵ)
    push!(numsims, its)

    if ((( abs(ϵvec[end - 1] - ϵ )) / ϵvec[end - 1]) < ABCsetup.convergence) == true
      println("New ϵ is within $(round(ABCsetup.convergence * 100, digits = 2))% of previous population, stop ABC SMC")
      break
    end

    if savepopulations == true
        println()
        println("Saving population $popnum data and plots\n")
        abcres = ApproxBayes.ABCSMCmodelresults(particles, numsims, ABCsetup, ϵvec)
        makedirectories(joinpath(resultsdirectory, sname, "populations", "population_$popnum"))
        posteriors, DFmp = getresults(abcres, resultsdirectory, sname, VAF, save = true, Nmaxinf = Nmaxinf, savepopulations = savepopulations, popnum = popnum);
        populationresults = Results(ABCsetup, abcres, VAF, posteriors, DFmp, sname)
        saveallplots(populationresults; resultsdirectory = resultsdirectory, outputformat = ".pdf", savepopulations = savepopulations, popnum = popnum)
    end

    popnum = popnum + 1
    if verbose
      show(ApproxBayes.ABCSMCmodelresults(oldparticles, numsims, ABCsetup, ϵvec))
    end

    #sometimes models die out early in the inference, restart the process if this happens
    numdeadmodels = sum(modelprob.==0.0)
    deadmodels = collect(0:(ABCsetup.nmodels-1))[modelprob.==0.0]
    if (popnum <= 4) & (numdeadmodels > 0)
      warn("One of the models died out, restarting inference...")
      out = runabcCancer(ABCsetup, targetdata; verbose = verbose, progress = progress)
      return out
    end
  end

  out = ApproxBayes.ABCSMCmodelresults(particles, numsims, ABCsetup, ϵvec)
  return out
end

function runabcCancer(ABCsetup::ABCSMC, targetdata; verbose = false, progress = false)

  #run first population with parameters sampled from prior
  if verbose
    println("##################################################")
    println("Use ABC rejection to get first population")
  end
  ABCrejresults = ApproxBayes.runabc(ApproxBayes.ABCRejection(ABCsetup.simfunc, ABCsetup.nparams,
                  ABCsetup.ϵ1, ABCsetup.prior; nparticles = ABCsetup.nparticles,
                  maxiterations = ABCsetup.maxiterations, constants = ABCsetup.constants), targetdata);

  oldparticles, weights = ApproxBayes.setupSMCparticles(ABCrejresults, ABCsetup)
  ϵ = quantile(ABCrejresults.dist, ABCsetup.α) # set new ϵ to αth quantile
  ϵvec = [ϵ] #store epsilon values
  numsims = [ABCrejresults.numsims] #keep track of number of simualtions
  particles = Array{ApproxBayes.ParticleSMC}(ABCsetup.nparticles) #define particles array
  ABCsetup = ApproxBayes.modelselection_kernel(ABCsetup, oldparticles)

  if verbose
    println("Run ABC SMC \n")
  end

  popnum = 1
  finalpop = false
  newparticle, dist, out, priorp = 0.0,0.0,0.0,0.0

  while (ϵ > ABCsetup.ϵT) & (sum(numsims) < ABCsetup.maxiterations)

    i = 1 #set particle indicator to 1
    particles = Array{ApproxBayes.ParticleSMC}(ABCsetup.nparticles)
    distvec = zeros(Float64, ABCsetup.nparticles)
    its = 1
    if progress == true
      p = Progress(ABCsetup.nparticles, 1, "ABC SMC population $(popnum), new ϵ: $(round(ϵ, digits = 2))...", 30)
    end
    while i < ABCsetup.nparticles + 1

      correctmodel = false
      while correctmodel == false
        j = wsample(1:ABCsetup.nparticles, weights)
        particle = oldparticles[j]
        newparticle = ApproxBayes.perturbparticle(particle, ABCsetup.kernel)
        priorp = ApproxBayes.priorprob(newparticle.params, ABCsetup.prior)
        if priorp == 0.0 #return to beginning of loop if prior probability is 0
          break
        end
        #simulate with new parameters
        dist, out, cm = ABCsetup.simfunc(newparticle.params, ABCsetup.constants, targetdata)
        correctmodel = cm
      end

    if priorp == 0.0 #return to beginning of loop if prior probability is 0
      continue
    end

      #if simulated data is less than target tolerance accept particle
      if dist < ϵ
        particles[i] = newparticle
        particles[i].other = out
        particles[i].distance = dist
        distvec[i] = dist
        i += 1
        if progress == true
          next!(p)
        end
      end

      its += 1
    end

    particles, weights = ApproxBayes.smcweights(particles, oldparticles, ABCsetup.prior)
    oldparticles = particles
    ABCsetup = ApproxBayes.modelselection_kernel(ABCsetup, oldparticles)

    if finalpop == true
      break
    end

    ϵ = quantile(distvec, ABCsetup.α)

    if ϵ < ABCsetup.ϵT
      ϵ = ABCsetup.ϵT
      push!(ϵvec, ϵ)
      push!(numsims, its)
      popnum = popnum + 1
      finalpop = true
      continue
    end

    push!(ϵvec, ϵ)
    push!(numsims, its)

    if ((( abs(ϵvec[end - 1] - ϵ )) / ϵvec[end - 1]) < ABCsetup.convergence) == true
      if verbose
        println("New ϵ is within $(round(ABCsetup.convergence * 100, digits = 2))% of previous population, stop ABC SMC")
    end
      break
    end

    popnum = popnum + 1
  end

  out = ApproxBayes.ABCSMCresults(particles, numsims, ABCsetup, ϵvec)
  return out
end


text:
"""
```
ct_kalman_filter(y, T, R, C, Q, Z, D, E, tspan, s_0 = Vector(), P_0 = Matrix();
    outputs = [:loglh, :pred, :filt], Nt0 = 0)
tspan - Number of time units betwen each observation
"""
function ct_kalman_filter(y::Matrix{S},
    T::Matrix{S}, R::Matrix{S}, C::Vector{S},
    Q::Matrix{S}, Z::Matrix{S}, D::Vector{S}, E::Matrix{S}, tspan::Float64,
    s_0::Vector{S} = Vector{S}(undef, 0), P_0::Matrix{S} = Matrix{S}(undef, 0, 0);
    outputs::Vector{Symbol} = [:loglh, :pred, :filt],
    Nt0::Int = 0, method = Tsit5(),
    reltol::Float64 = 1e-8, abstol::Float64 = 1e-8) where {S<:AbstractFloat}

    # Determine outputs
    return_loglh = :loglh in outputs
    return_pred  = :pred  in outputs
    return_filt  = :filt  in outputs

    # Dimensions
    Ns = size(T, 1) # number of states
    Nt = size(y, 2) # number of periods of data

    # Initialize inputs and outputs
    k = StateSpaceRoutines.KalmanFilter(T, R, C, Q, Z, D, E, s_0, P_0)

    mynan = convert(S, NaN)
    loglh  = return_loglh ? fill(mynan, Nt)         : Vector{S}(0)
    s_pred = return_pred  ? fill(mynan, Ns, Nt)     : Matrix{S}(0, 0)
    P_pred = return_pred  ? fill(mynan, Ns, Ns, Nt) : Array{S, 3}(0, 0, 0)
    s_filt = return_filt  ? fill(mynan, Ns, Nt)     : Matrix{S}(0, 0)
    P_filt = return_filt  ? fill(mynan, Ns, Ns, Nt) : Array{S, 3}(0, 0, 0)

    # Populate initial states
    s_0 = k.s_t
    P_0 = k.P_t

    # Loop through periods t
    for t = 1:Nt
        # Forecast
        forecast!(k, tspan; method = method, reltol = reltol, abstol = abstol)
        if return_pred
            s_pred[:,    t] = k.s_t
            P_pred[:, :, t] = k.P_t
        end

        # Update and compute log-likelihood
        StateSpaceRoutines.update!(k, y[:, t]; return_loglh = return_loglh)
        if return_filt
            s_filt[:,    t] = k.s_t
            P_filt[:, :, t] = k.P_t
        end
        if return_loglh
            loglh[t]        = k.loglh_t
        end

        # Update s_0 and P_0 if Nt0 > 0
        if t == Nt0
            s_0 = k.s_t
            P_0 = k.P_t
        end
    end

    # Populate final states
    s_T = k.s_t
    P_T = k.P_t

    # Remove presample periods
    loglh, s_pred, P_pred, s_filt, P_filt =
        StateSpaceRoutines.remove_presample!(Nt0, loglh, s_pred, P_pred, s_filt, P_filt; outputs = outputs)

    return loglh, s_pred, P_pred, s_filt, P_filt, s_0, P_0, s_T, P_T
end

"""
```
ct_kalman_filter(y, T, R, C, Q, Z, D, E, tspan, s_0 = Vector(), P_0 = Matrix();
    outputs = [:loglh, :pred, :filt], Nt0 = 0)
tspan - Number of time units betwen each observation
n_subinterval - Number of subintervals to use ODE integration with
"""
function ct_kalman_filter(y::Matrix{S},
                          T::Matrix{S}, R::Matrix{S}, C::Vector{S}, Q::Matrix{S},
                          Z::Matrix{S}, D::Vector{S}, E::Matrix{S}, tspan::Float64, n_subinterval::Int64,
                          s_0::Vector{S} = Vector{S}(0), P_0::Matrix{S} = Matrix{S}(0, 0);
                          outputs::Vector{Symbol} = [:loglh, :pred, :filt],
                          Nt0::Int = 0, method = Tsit5(),
                          reltol::Float64 = 1e-8, abstol::Float64 = 1e-8) where {S<:AbstractFloat}

    # Determine outputs
    return_loglh = :loglh in outputs
    return_pred  = :pred  in outputs
    return_filt  = :filt  in outputs

    # Dimensions
    Ns = size(T, 1) # number of states
    Nt = size(y, 2) # number of periods of data

    # Initialize inputs and outputs
    k = StateSpaceRoutines.KalmanFilter(T, R, C, Q, Z, D, E, s_0, P_0)

    mynan = convert(S, NaN)
    loglh  = return_loglh ? fill(mynan, Nt)         : Vector{S}(0)
    s_pred = return_pred  ? fill(mynan, Ns, Nt)     : Matrix{S}(0, 0)
    P_pred = return_pred  ? fill(mynan, Ns, Ns, Nt) : Array{S, 3}(0, 0, 0)
    s_filt = return_filt  ? fill(mynan, Ns, Nt)     : Matrix{S}(0, 0)
    P_filt = return_filt  ? fill(mynan, Ns, Ns, Nt) : Array{S, 3}(0, 0, 0)

    # Populate initial states
    s_0 = k.s_t
    P_0 = k.P_t

    # Loop through periods t
    for t = 1:Nt
        # Forecast
        forecast!(k, tspan, n_subinterval; method = method, reltol = reltol, abstol = abstol)
        if return_pred
            s_pred[:,    t] = k.s_t
            P_pred[:, :, t] = k.P_t
        end

        # Update and compute log-likelihood
        StateSpaceRoutines.update!(k, y[:, t]; return_loglh = return_loglh)
        if return_filt
            s_filt[:,    t] = k.s_t
            P_filt[:, :, t] = k.P_t
        end
        if return_loglh
            loglh[t]        = k.loglh_t
        end

        # Update s_0 and P_0 if Nt0 > 0
        if t == Nt0
            s_0 = k.s_t
            P_0 = k.P_t
        end
    end

    # Populate final states
    s_T = k.s_t
    P_T = k.P_t

    # Remove presample periods
    loglh, s_pred, P_pred, s_filt, P_filt =
        remove_presample!(Nt0, loglh, s_pred, P_pred, s_filt, P_filt; outputs = outputs)

    return loglh, s_pred, P_pred, s_filt, P_filt, s_0, P_0, s_T, P_T
end

"""
Forecast functions
"""
function forecast!(k::StateSpaceRoutines.KalmanFilter, tspan::Float64; method = Tsit5(),
                   reltol::Float64 = 1e-8, abstol::Float64 = 1e-8) where {S<:AbstractFloat}
    T, R, Q = k.T, k.R, k.Q
    s0, P0 = k.s_t, k.P_t
    f(u,p,t) = T*u
    g(u,p,t) = T*u + u*T' + R*Q*R'
    tspan_int = (0., tspan)
    s_prob = ODEProblem(f, s0, tspan_int)
    P_prob = ODEProblem(g, P0, tspan_int)
    s_sol = DifferentialEquations.solve(s_prob, method, reltol = reltol, abstol = abstol)
    P_sol = DifferentialEquations.solve(P_prob, method, reltol = reltol, abstol = abstol)

    k.s_t = s_sol.u[end]
    k.P_t = P_sol.u[end]
end

# INCOMPLETE
# Not sure how to write down the ODE system correctly since I don't have the same
# kind of recursion that I can do otherwise. Moreover, I effectively will get the
# same predicted states regardless, so I'm not sure if this will improve the
# measurement equation or accuracy.
function forecast!(k::StateSpaceRoutines.KalmanFilter, tspan::Float64, n_subinterval::Int64; method = Tsit5(),
                   reltol::Float64 = 1e-8, abstol::Float64 = 1e-8) where {S<:AbstractFloat}
    T, R, Q = k.T, k.R, k.Q
    s_filt, P0 = k.s_t, k.P_t
    indiv_dim = size(T, 1)
    s0 = s_filt[1:indiv_dim]
    f(u,p,t) = T*u
    g(u,p,t) = T*u + u*T' + R*Q*R'
    tspan_int = (0., span)
    s_prob = ODEProblem(f, s0, tspan_int)
    P_prob = ODEProblem(g, P0, tspan_int)
    s_sol = solvect(s_prob, method, reltol = reltol, abstol = abstol)
    P_sol = solvect(P_prob, method, reltol = reltol, abstol = abstol)
    k.s_t = s_sol.u
    k.P_t = P_sol.u
end


text:
# Scale from/to [0,1].
scalefrom(a,b,t) = @. (t-a)/(b-a)
scaleto(a,b,t) = @. a + t*(b-a)

# Unique real roots of a quadratic.
function realroots(a,b,c)
	a==0 && return [-c/b]
	bh = -b/2
	d = bh^2 - a*c
	if d < 0
		[]
	elseif d==0
		[bh/a]
	else
		r = (bh + sign(bh)*sqrt(d)) / a
		[r,c/(r*a)]
	end
end

# Are three given points arranged in counterclockwise order?
function isccw(a::Number,b::Number,c::Number)
	v(z) = SVector(1,real(z),imag(z))
	det([v(a) v(b) v(c)]) > 0
end

# Use 2nd order finite differences to approximate a tangent.
function fdtangent(z,t::Real)
	ϵ = eps(typeof(float(t)))
	ϵ3 = 0.5*ϵ^(1/3)
	if t < ϵ3
		τ = (-1.5*z(t) + 2*z(t+ϵ3) - 0.5*z(t+2ϵ3)) / ϵ3
	elseif t > 1-ϵ3
		τ = (1.5*z(t) - 2*z(t-ϵ3) + 0.5*z(t-2ϵ3)) / ϵ3
	else
		τ = (z(t+ϵ3) - z(t-ϵ3)) / (2ϵ3)
	end
	return τ
end

# Select points adaptively to make a smooth-appearing curve.
function adaptpoints(point,utangent,a,b;depth=6,curvemax=0.05)
	function refine(tl,tr,zl,zr,τl,τr,maxdz,d=depth)
		dzkap = dist(τr,τl)  # approximately, the stepsize over radius of curvature
		tm = (tl+tr)/2
		zm = point(tm)
		τm = utangent(tm)
		if d > 0 && (dzkap > curvemax || dist(zr,zl) > maxdz )
			zl = refine(tl,tm,zl,zm,τl,τm,maxdz,d-1)
			zr = refine(tm,tr,zm,zr,τm,τr,maxdz,d-1)
			return [zl;zm;zr]
		else
			return zm
		end
	end

	d = (b-a)/4
	tt = d*[0,0.196,0.41,0.592,0.806]   # avoid common symmetry points
	t = [a .+ tt; a + d .+ tt; a + 2d .+ tt; a + 3d .+ tt; b]
	z = point.(t)
	τ = utangent.(t)

	# on the Riemann sphere, use distance in R^3
	if z[1] isa Spherical
		dist = (u,v) -> norm(S2coord(u)-S2coord(v))
	else
		dist = (u,v) -> abs(u-v)
	end
	m = length(t)
	scale = maximum(dist(z[i],z[j]) for i=2:m-1, j=2:m-1 if j > i)
	zfinal = z[[1]]
	for j = 1:length(t)-1
		znew = refine(t[j],t[j+1],z[j],z[j+1],τ[j],τ[j+1],scale/25)
		append!(zfinal,znew)
		push!(zfinal,z[j+1])
	end
	return zfinal
end

# Do adaptive integration to estimate the integral of `f` over [`a`,`b`] to desired
# error tolerance `tol`.
function intadapt(f,a,b,tol)
    # Use error estimation and recursive bisection.
    function do_integral(a,fa,b,fb,m,fm,tol,depth)
        # These are the two new nodes and their f-values.
        xl = (a+m)/2;  fl = f(xl);
        xr = (m+b)/2;  fr = f(xr);
        t = [a,xl,m,xr,b]              # all 5 nodes at this level

        # Compute the trapezoid values iteratively.
        h = (b-a)
        T = [0.,0.,0.]
        T[1] = h*(fa+fb)/2
        T[2] = T[1]/2 + (h/2)*fm
        T[3] = T[2]/2 + (h/4)*(fl+fr)

        S = (4*T[2:3]-T[1:2]) / 3      # Simpson values
        E = (S[2]-S[1]) / 15           # error estimate

        if abs(E) < tol*(1+abs(S[2]))  # acceptable error?
            Q = S[2]                   # yes--done
		else
			if depth==0
				@warn "Too many recursions to determine integral"
				Q = S[2]
			else
    	        # Error is too large--bisect and recurse.
        	    QL = do_integral(a,fa,m,fm,xl,fl,tol,depth-1)
            	QR = do_integral(m,fm,b,fb,xr,fr,tol,depth-1)
				Q = QL + QR
			end
        end
        return Q
    end

    m = (b+a)/2
    Q = do_integral(a,f(a),b,f(b),m,f(m),tol,50)
    return Q
end

function enclosing_circle(z::AbstractVector,expansion=2)
	xa,xb = extrema(real(z))
	ya,yb = extrema(imag(z))
	zc = complex((xa+xb)/2,(ya+yb)/2)
	R = length(z) > 1 ? maximum(@. abs(z - zc)) : max(1,abs(zc))
	return zc,expansion*R
end

function enclosing_box(z::AbstractVector,expansion=2)
	zc = sum(z)/length(z)
	dz = z .- zc
	rx = length(z) > 1 ? maximum(@. abs(real(dz))) : max(1,abs(real(zc)))
	ry = length(z) > 1 ? maximum(@. abs(imag(dz))) : max(1,abs(imag(zc)))
	return real(zc).+expansion*[-rx,rx],imag(zc).+expansion*[-ry,ry]
end

# indices of the closest pair of points from two lists
function argclosest(z1,z2)
	i1 = [argmin(abs.(z1.-z)) for z in z2]
	i2 = argmin(abs.(z2.-z1[i1]))
	return i1[i2],i2
end


text:
"""
    @repeat times begin
        body
    end

Repeatedly evaluate body for the specified number of times.
"""
macro repeat(times, body)
    quote
        for i = 1:$times
            $body
        end
    end |> esc
end


text:
# Julia wrapper for header: /usr/local/include/libavutil/mastering_display_metadata.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    av_mastering_display_metadata_alloc,
    av_mastering_display_metadata_create_side_data

function av_mastering_display_metadata_alloc()
    ccall((:av_mastering_display_metadata_alloc,libavutil),Ptr{AVMasteringDisplayMetadata},())
end

function av_mastering_display_metadata_create_side_data(frame)
    ccall((:av_mastering_display_metadata_create_side_data,libavutil),Ptr{AVMasteringDisplayMetadata},(Ptr{AVFrame},),frame)
end

max_stars_repo_name: asafasban/AAVideoIO

max_issues_repo_name: asafasban/AAVideoIO

max_forks_repo_name: asafasban/AAVideoIO


text:
@testset "Testing FindNextTree" begin
    # test ceiling integer division
    @test CachedArrays.cdiv(64, 64) == 1
    @test CachedArrays.cdiv(63, 64) == 1
    @test CachedArrays.cdiv(65, 64) == 2

    @test CachedArrays.cdiv(128, 64) == 2
    @test CachedArrays.cdiv(127, 64) == 2
    @test CachedArrays.cdiv(129, 64) == 3

    @testset "Mask" begin
        m = CachedArrays.Mask(0)
        @test CachedArrays.hasentry(m) == false
        for i in 1:64
            @test CachedArrays.hasentry(m) == false
        end

        for i in 1:64
            # Test with single bits set.
            m = CachedArrays.Mask(one(UInt64) << (i-1))
            @test CachedArrays.hasentry(m)
            @test CachedArrays.hasentryat(m, i)
            for j in 1:i-1
                @test CachedArrays.hasentry(m, j) == true
            end
            for j in i:64
                @test CachedArrays.hasentry(m, j) == false
            end

            # Now, set all kinds of bits
            m = CachedArrays.Mask((one(UInt64) << (i - 1)) - 1)
            # If i == 1, then the mask value is, in fact, zero
            if i > 1
                @test CachedArrays.hasentry(m)
            end

            for j in 1:i-2
                @test CachedArrays.hasentry(m, j) == true
            end
            for j in i-1:64
                @test CachedArrays.hasentry(m, j) == false
            end
        end
    end

    @testset "FindNextTree" begin
        # Start simple with just a 1-level tree
        tree = CachedArrays.FindNextTree(64)
        @test length(tree.runs) == 1

        for i in 1:64
            CachedArrays.setentry!(tree, i)
            @test CachedArrays.hasentryat(tree.runs[end][1], i) == true
            CachedArrays.clearentry!(tree, i)
            @test CachedArrays.hasentryat(tree.runs[end][1], i) == false
        end

        # For now, test corner cases of interest.
        # Later, we'll move onto a longer stress test.
        setentry!(x...) = CachedArrays.setentry!(x...)
        clearentry!(x...) = CachedArrays.clearentry!(x...)
        numentries(x...) = CachedArrays.numentries(x...)

        # Okay, now we get a little wilder
        # Use a bit-array to provide a test against.
        len = 64^3
        tree = CachedArrays.FindNextTree(len)
        @test length(tree.runs) == 3

        # Right now, the top most level of the tree should have no entries.
        @test all(x -> iszero(numentries(x)), tree.runs[1])
        setentry!(tree, 1)

        # Make sure the top level of the tree got set.
        @test numentries(tree.runs[1][1]) == 1
        @test numentries(tree.runs[2][1]) == 1

        # If we set another bit in the same bin, it shouldn't change anything.
        setentry!(tree, 2)
        @test numentries(tree.runs[2][1]) == 1
        @test findnext(tree, 1) == 2
        @test findnext(tree, 2) == nothing

        setentry!(tree, 64)
        @test findnext(tree, 1) == 2
        @test findnext(tree, 2) == 64
        clearentry!(tree, 64)
        @test findnext(tree, 1) == 2
        @test findnext(tree, 2) == nothing

        setentry!(tree, 120)
        @test findnext(tree, 2) == 120
        clearentry!(tree, 120)
        @test findnext(tree, 2) == nothing

        # Set the entry at the end of the mask
        setentry!(tree, 64^3)
        @test findnext(tree, 2) == 64^3
        @test findnext(tree, 1) == 2
        clearentry!(tree, 64^3)
        @test findnext(tree, 2) == nothing
        clearentry!(tree, 1)
        @test findnext(tree, 1) == 2
        clearentry!(tree, 2)
        @test findnext(tree, 1) == nothing

        # Start stress testing.
        Random.seed!(1234)

        reference = falses(len)
        setentries = ones(Int64, 10000)
        search_indices = ones(Int64, 100000)

        @time for iter in 1:10
            # Populate indices to set
            for i in eachindex(setentries)
                setentries[i] = rand(1:len)
            end

            # Set entries
            for entry in setentries
                reference[entry] = true
                setentry!(tree, entry)
            end

            # Populate search indices
            for i in eachindex(search_indices)
                search_indices[i] = rand(1:len)
            end

            # Fuzz test
            for i in search_indices
                # There's a slight difference in behavior between the tree and the reference.
                # If an entry is set, `findnext(::BitArray, i)` will return `i`, where-as
                # the tree will still search for the next.
                #
                # So, if this entry is set, we expect the behavior to be different and skip
                # this test.
                reference[i] && continue

                if findnext(reference, i) != findnext(tree, i)
                    message = """
                    Failed on iteration $iter and index $i
                    Reference: $(findnext(reference, i))
                    Tree: $(findnext(tree, i))
                    """
                    error(message)
                end
            end

            # Clear entries for next round.
            for entry in setentries
                reference[entry] = false
                clearentry!(tree, entry)
            end
        end
    end
end


text:
using TimeZones: optional

function strip(ex::Expr)
    args = []
    for arg in ex.args
        if isa(arg, Expr)
            arg = strip(arg)
        end

        if !isa(arg, LineNumberNode)
            push!(args, arg)
        end
    end

    return Expr(ex.head, args...)
end

Base.isequal(a::Array{Expr}, b::Array{Expr}) = map(strip, a) == map(strip, b)

# We can include comments in our comparison here since the expressions are on the same line
@test optional(:(foo(a, b, c=3; d=4) = nothing)) == [:(foo(a, b, c=3; d=4) = nothing)]

@test isequal(
    optional(
        :(foo(a, b=2, c) = nothing)
    ),
    [
        :(foo(a, b, c) = nothing),
        :(foo(a, c) = foo(a, 2, c)),
    ],
)

@test isequal(
    optional(
        :(foo(a, b=2, c=3, d, e=5; f=6) = nothing)
    ),
    [
        :(foo(a, b, c, d, e=5; f=6) = nothing),
        :(foo(a, b, d, e=5; f=6) = foo(a, b, 3, d, e; f=f)),
        :(foo(a, d, e=5; f=6) = foo(a, 2, 3, d, e; f=f)),
    ],
)

const I = Integer
@test isequal(
    optional(
        :(function ZonedDateTime(y::I, m::I=1, d::I=1, h::I=0, mi::I=0, s::I=0, ms::I=0, tz::TimeZone)
            ZonedDateTime(DateTime(y,m,d,h,mi,s,ms), tz)
        end)
    ),
    [
        :(function ZonedDateTime(y::I,m::I,d::I,h::I,mi::I,s::I,ms::I,tz::TimeZone)
            ZonedDateTime(DateTime(y,m,d,h,mi,s,ms), tz)
        end),
        :(ZonedDateTime(y::I,m::I,d::I,h::I,mi::I,s::I,tz::TimeZone) = ZonedDateTime(y,m,d,h,mi,s,0,tz)),
        :(ZonedDateTime(y::I,m::I,d::I,h::I,mi::I,tz::TimeZone) = ZonedDateTime(y,m,d,h,mi,0,0,tz)),
        :(ZonedDateTime(y::I,m::I,d::I,h::I,tz::TimeZone) = ZonedDateTime(y,m,d,h,0,0,0,tz)),
        :(ZonedDateTime(y::I,m::I,d::I,tz::TimeZone) = ZonedDateTime(y,m,d,0,0,0,0,tz)),
        :(ZonedDateTime(y::I,m::I,tz::TimeZone) = ZonedDateTime(y,m,1,0,0,0,0,tz)),
        :(ZonedDateTime(y::I,tz::TimeZone) = ZonedDateTime(y,1,1,0,0,0,0,tz)),
    ],
)

# Currently demonstrates an issue when the type given doesn't match value
# @test isequal(
#     optional(
#         :(f(a::Float64=1, b) = nothing)
#     ),
#     [
#         :(f(a::Float64, b) = nothing),
#         :(f(b) = f(1.0, b)),
#     ]
# )


text:
# Populate a SeisData structure with random channels
"""
    randSeisData()

Generate 8 to 24 channels of random seismic data as a SeisData object.

    randSeisData(N)

Generate `N` channels of random seismic data as a SeisData object.

### Keywords
| KW      | Default | Type      | Meaning                                 |
|-------- |:--------|:----------|:-------------------------------         |
| s       | 0.6     | Float64   | chance each channel is seismic data     |
| c       | 0.2     | Float64   | chance each channel is irregular (fs=0) |
| nx      | 0       | Int64     | number of samples in each channel       |
| fs_min  | 0.0     | Float64   | regular channels have fs ≥ fs_min       |
| a0      | false   | Bool      | recalculate :resp constants?            |

### Behavior
* If `nx` ≤ 0, the number of samples is determined randomly for each channel
* `100*s` is the minimum percentage of channels with guaranteed seismic data
* `100*c` is the maximum percentage of channels with irregularly-sampled data
* `s` takes precedence over `c`; they are not renormalized. Thus, for example:
  + `randSeisData(12, c=1.0, s=1.0)` works like `randSeisData(12, c=0.0, s=1.0)`
  + `randSeisData(c=1.0)` works like `randSeisData(s=0.0)`
* For `N` channels, `s*N` is rounded up; `c*N` is rounded down. Thus:
  + `randSeisData(10, c=0.28)` and `randSeisData(10, c=0.2)` are equivalent
  + `randSeisData(10, s=0.28)` and `randSeisData(10, s=0.3)` are equivalent
* By default, `:resp.a0 = 1.0` for all channels; `a0=true` calls `resp_a0!(S)`
  + `:resp.p` and `:resp.z` are random; `a0=true` isn't guaranteed to work

See also: `randSeisChannel`, `resp_a0!`
"""
function randSeisData(N::Int64; c::Float64=0.2, s::Float64=0.6, nx::Integer=0, fs_min::Float64=0.0, a0::Bool=false)
  S = SeisData(N)

  # Evaluate probabilities, with s taking precedence
  if s != 0.6
    c = max(0.0, 1.0-s)
  end
  if c != 0.2
    s = max(0.0, 1.0-c)
  end

  # determine number of channels of each
  n_seis = max(min(ceil(Int, s*S.n), S.n), 0)
  n_irr = max(min(floor(Int, c*S.n), S.n-n_seis), 0)

  # fill in and shuffle data_spec
  data_spec = zeros(UInt8, S.n)
  data_spec[1:n_seis] .= 0x01
  data_spec[n_seis+1:n_seis+n_irr] .= 0x02
  shuffle!(data_spec)

  # populate all channels according to data_spec
  for i = 1:S.n
    if data_spec[i] == 0x01
      S[i] = randSeisChannel(s=true, nx=nx, fs_min=fs_min)
    elseif data_spec[i] == 0x02
      S[i] = randSeisChannel(c=true, nx=nx, fs_min=fs_min)
    else
      S[i] = randSeisChannel(nx=nx, fs_min=fs_min)
    end
  end
  a0 && resp_a0!(S)
  return S
end
randSeisData(; c::Float64=0.2, s::Float64=0.6, nx::Int64=0, fs_min::Float64=0.0, a0::Bool=false) = randSeisData(rand(8:24), c=c, s=s, nx=nx, fs_min=fs_min, a0=a0)


text:
function replace_invalid_character(string::AbstractString)
    substitutions = ["\\" => "/",
                     "\"" => "'",
                     "\n" => "\\n",
                     "\r" => "\\r"]
    for sub in substitutions
        string = replace(string, sub)
    end
    return string
end

"""
    stringify(; kwargs...)

Convert keywords arguments into a string that conforms the log event message structure used in the `Seq` logger.
"""
function stringify(; kwargs...)
    event_properties = Vector{String}()
    for kwarg in kwargs
        # kwarg.second could have invalid characters, kwarg.first not
        clean_kwarg_second = replace_invalid_character("$(kwarg.second)")
        property = "\"$(kwarg.first)\":\"$(clean_kwarg_second)\""
        push!(event_properties, property)
    end
    return join(event_properties, ",")
end

"""
    to_seq_level(log_level::Base.CoreLogging.LogLevel)

Given a `Logging.LogLevel` return the string identifier for the log level used on the `Seq` server.
"""
function to_seq_level(log_level::Base.CoreLogging.LogLevel)
    if log_level == Logging.Debug
        return "Debug"
    elseif log_level == Logging.Info
        return "Info"
    elseif log_level == Logging.Warn
        return "Warning"
    elseif log_level == Logging.Error
        return "Error"
    end
end

"""
    joinurl(left::AbstractString, right::AbstractString)

Join the left and right part of a URL, by removing trailing frontslashes and add
an additional frontslashes if required.
"""
function joinurl(left::AbstractString, right::AbstractString)
    leftStripped = rstrip(left, '/')
    rightStripped = lstrip(right, '/')
    return join([leftStripped, rightStripped], '/')
end


text:
#Is this usable data?
isconcrete(x::T) where {T<:Number} = !ismissing(x) && isfinite(x)
isconcrete(x::(Irrational)) = true
isconcrete(x) = !ismissing(x)

function isallconcrete(xs)
    ans = true
    for x in xs
        if !isconcrete(x)
            ans = false
            break
        end
    end
    return ans
end

function concretize(xss::AbstractVector...)
    if all(map(isallconcrete, xss))
        return xss
    end

    @info "concretizing" length.(xss) xss...

    count = 0
    for j in 1:length(xss[1])
        for xs in xss
            if !isconcrete(xs[j])
                @goto next_j1
            end
        end

        count += 1
        @label next_j1
    end

    yss = Vector{AbstractVector}(undef, length(xss))
    for (i, xs) in enumerate(xss)
        yss[i] = Vector{eltype(xs)}(undef, count)
    end

    k = 1
    for j in 1:length(xss[1])
        for xs in xss
            if !isconcrete(xs[j])
                @goto next_j2
            end
        end

        for (i, xs) in enumerate(xss)
            yss[i][k] = xs[j]
        end
        k += 1

        @label next_j2
    end

    @info "returning" yss...

    return tuple(yss...)
end

# How many concrete elements in an iterable
function concrete_length(xs)
    return sum(isconcrete.(skipmissing(xs)))
end

function concrete_length(xs::Iterators.Flatten)
    n = 0
    for obj in xs.it
        n += concrete_length(obj)
    end
    n
end

function concrete_minimum(xs)
    if isempty(xs)
        error("argument must not be empty")
    end

    x_min = first(xs)
    for x in xs
        if Gadfly.isconcrete(x) && isfinite(x)
            x_min = x
            break
        end
    end

    for x in xs
        if Gadfly.isconcrete(x) && isfinite(x) && x < x_min
            x_min = x
        end
    end
    return x_min
end

function concrete_maximum(xs)
    if isempty(xs)
        error("argument must not be empty")
    end

    x_max = first(xs)
    for x in xs
        if Gadfly.isconcrete(x) && isfinite(x)
            x_max = x
            break
        end
    end

    for x in xs
        if Gadfly.isconcrete(x) && isfinite(x) && x > x_max
            x_max = x
        end
    end
    return x_max
end

function concrete_minmax(xs, xmin::T, xmax::T) where T<:Real
    if eltype(xs) <: Base.Callable
        return xmin, xmax
    end

    for x in xs
        if !ismissing(x) && isconcrete(x)
            xT = convert(T, x)
            if isnan(xmin) || xT < xmin
                xmin = xT
            end
            if isnan(xmax) || xT > xmax
                xmax = xT
            end
        end
    end
    xmin, xmax
end

function concrete_minmax(xs, xmin::T, xmax::T) where T
    for x in xs
        if !ismissing(x) && isconcrete(x)
            xT = convert(T, x)
            if xT < xmin
                xmin = xT
            end
            if xT > xmax
                xmax = xT
            end
        end
    end
    xmin, xmax
end

function concrete_minmax(xs::Iterators.Flatten, xmin::T, xmax::T) where T<:Real
    for obj in xs.it
        xmin, xmax = concrete_minmax(obj, xmin, xmax)
    end
    xmin, xmax
end

# Create a new object of type T from a with missing values (i.e., those set to
# nothing) inherited from b.
function inherit(a::T, b::T) where T
    c = copy(a)
    inherit!(c, b)
    c
end

function inherit!(a::T, b::T) where T
    for field in fieldnames(T)
        aval = getfield(a, field)
        bval = getfield(b, field)
        # TODO: this is a hack to let non-default labelers overide the defaults
        if aval === nothing || aval === string || aval == showoff
            setfield!(a, field, bval)
        elseif typeof(aval) <: Dict && typeof(bval) <: Dict
            merge!(aval, getfield(b, field))
        end
    end
    nothing
end

isnothing(u) = u === nothing
issomething(u) = !isnothing(u)

negate(f) = x -> !f(x)

function has(xs::AbstractArray{T,N}, y::T) where {T,N}
    for x in xs
        if x == y
            return true
        end
    end
    return false
end

Maybe(T) = Union{T, (Nothing)}

lerp(x::Float64, a, b) = a + (b - a) * max(min(x, 1.0), 0.0)

# Generate a unique id, primarily for assigning IDs to SVG elements.
let next_unique_svg_id_num = 0
    global unique_svg_id
    function unique_svg_id()
        uid = @sprintf("id%d", next_unique_svg_id_num)
        next_unique_svg_id_num += 1
        uid
    end
end

# Remove any markup or whitespace from a string.
function escape_id(s::AbstractString)
    s = replace(s, r"<[^>]*>" => "")
    s = replace(s, r"\s" => "_")
    s
end

# Can a numerical value be treated as an integer
is_int_compatable(::Integer) = true
is_int_compatable(x::T) where {T <: AbstractFloat} = abs(x) < maxintfloat(T) && float(int(x)) == x
is_int_compatable(::Any) = false

# Construct a jscall to store arbitrary data in the element
function jsdata(key::AbstractString, value::AbstractString, arg::Vector{Measure}=Measure[])
    return jscall(
        """
        data("$(escape_string(key))", $(value))
        """, arg)
end

# Construct jscall properties to store arbitrary data in plotroot elements.
function jsplotdata(key::AbstractString, value::AbstractString, arg::Vector{Measure}=Measure[])
    return jscall(
        """
        plotroot().data("$(escape_string(key))", $(value))
        """, arg)
end

svg_color_class_from_label(label::AbstractString) = @sprintf("color_%s", escape_id(label))

using Dates

# Arbitrarily order colors
color_isless(a::Color, b::Color) =
        color_isless(convert(RGB{Float32}, a), convert(RGB{Float32}, b))
color_isless(a::TransparentColor, b::TransparentColor) =
        color_isless(convert(RGBA{Float32}, a), convert(RGBA{Float32}, b))

function color_isless(a::RGB{Float32}, b::RGB{Float32})
    if a.r < b.r
        return true
    elseif a.r == b.r
        if a.g < b.g
            return true
        elseif a.g == b.g
            return a.b < b.b
        else
            return false
        end
    else
        return false
    end
end

function color_isless(a::RGBA{Float32}, b::RGBA{Float32})
    if color_isless(color(a), color(b))
        return true
    elseif color(a) == color(b)
        return a.alpha < b.alpha
    else
        return false
    end
end

function group_color_isless(a::(Tuple{S, T}),
                            b::(Tuple{S, T})) where {S, T <: Colorant}
    if a[1] < b[1]
        return true
    elseif a[1] == b[1]
        return color_isless(a[2], b[2])
    else
        return false
    end
end

#trim longer arrays to the size of the smallest one
#and zip the arrays.
function trim_zip(xs...)
    mx = max(map(length, xs)...)
    mn = min(map(length, xs)...)
    if mx == mn
        zip(xs...)
    else
        zip([length(x) == mn ? x : x[1:mn] for x in xs]...)
    end
end

# Convenience constructors of IndirectArrays
function discretize_make_ia(values::AbstractVector, levels)
    index = something.(indexin(values, levels), 0)
    any(iszero, index) && throw(ArgumentError("values not in levels encountered"))
    IndirectArray(index, levels)
end
discretize_make_ia(values::AbstractVector)         = discretize_make_ia(values, unique(values))
discretize_make_ia(values::AbstractVector, ::Nothing) = discretize_make_ia(values)

discretize_make_ia(values::IndirectArray)         = values
discretize_make_ia(values::IndirectArray, ::Nothing) = values

discretize_make_ia(values::CategoricalArray) =
    discretize_make_ia(values, intersect(push!(levels(values), missing), unique(values)))
discretize_make_ia(values::CategoricalArray, ::Nothing) = discretize_make_ia(values)
function discretize_make_ia(values::CategoricalArray{T}, levels::Vector) where {T}
    mapping = something.(indexin(CategoricalArrays.index(values.pool), levels), 0)
    pushfirst!(mapping, something(findfirst(ismissing, levels), 0))
    index = [mapping[x+1] for x in values.refs]
    any(iszero, index) && throw(ArgumentError("values not in levels encountered"))
    return IndirectArray(index, convert(Vector{T},levels))
end
function discretize_make_ia(values::CategoricalArray{T}, levels::CategoricalVector{T}) where T
    _levels = map!(t -> ismissing(t) ? t : get(t), Vector{T}(length(levels)), levels)
    discretize_make_ia(values, _levels)
end


text:
#!/usr/bin/env julia
#
# This is an adapted version of format.jl from JuliaFormatter with the
# following license:
#
#    MIT License
#    Copyright (c) 2019 Dominique Luna
#
#    Permission is hereby granted, free of charge, to any person obtaining a
#    copy of this software and associated documentation files (the
#    "Software"), to deal in the Software without restriction, including
#    without limitation the rights to use, copy, modify, merge, publish,
#    distribute, sublicense, and/or sell copies of the Software, and to permit
#    persons to whom the Software is furnished to do so, subject to the
#    following conditions:
#
#    The above copyright notice and this permission notice shall be included
#    in all copies or substantial portions of the Software.
#
#    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
#    OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
#    MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN
#    NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,
#    DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
#    OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE
#    USE OR OTHER DEALINGS IN THE SOFTWARE.
#
using Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()

using JuliaFormatter

include("clima_formatter_options.jl")

help = """
Usage: climaformat.jl [flags] [FILE/PATH]...

Formats the given julia files using the CLIMA formatting options.  If paths
are given it will format the julia files in the paths. Otherwise, it will
format all changed julia files.

    -v, --verbose
        Print the name of the files being formatted with relevant details.

    -h, --help
        Print this message.
"""

function parse_opts!(args::Vector{String})
    i = 1
    opts = Dict{Symbol, Union{Int, Bool}}()
    while i ≤ length(args)
        arg = args[i]
        if arg[1] != '-'
            i += 1
            continue
        end
        if arg == "-v" || arg == "--verbose"
            opt = :verbose
        elseif arg == "-h" || arg == "--help"
            opt = :help
        else
            error("invalid option $arg")
        end
        if opt in (:verbose, :help)
            opts[opt] = true
            deleteat!(args, i)
        end
    end
    return opts
end

opts = parse_opts!(ARGS)
if haskey(opts, :help)
    write(stdout, help)
    exit(0)
end
if isempty(ARGS)
    filenames = readlines(`git diff --name-only --diff-filter=AM HEAD`)
    filter!(f -> endswith(f, ".jl"), filenames)
else
    filenames = ARGS
end

format(filenames; clima_formatter_options..., opts...)


text:
##
#
# Twindragon fractal
#
# Tomáš Kalvoda
# KAM FIT ČVUT, 2015
#

# helper type and functions
type SimpleSet
  container::Array{T, 2}
  maxlen::Int
  length::Int

  SimpleSet() = new(zeros(T, 2, array_length), array_length, 0)
end

@inline function empty!(set::SimpleSet)
  set.length = 0
end

@inline function append!(set::SimpleSet, x::Vector{T})
  set.length += 1
  @inbounds set.container[1, set.length] = x[1] # devectorization
  @inbounds set.container[2, set.length] = x[2]
end

@inline function initiate!(set::SimpleSet, x::Vector{T})
  empty!(set)
  append!(set, x)
end

@inline function multiply!(set::SimpleSet, m::Array{T}) #m is expected to by 2x2
  @fastmath @inbounds @simd for j in 1:length(set)
    a = set.container[1, j]; b = set.container[2, j];
    set.container[1, j] = m[1,1] * a + m[1,2] * b
    set.container[2, j] = m[2,1] * a + m[2,2] * b
  end
end

@inline length(set::SimpleSet) = set.length

@inline ball(z::Vector{T}) = z[1]^2 + z[2]^2 < radius

# twindragon iteration
function twindragon(z::Vector{T}, oldset::SimpleSet, newset::SimpleSet, in_domain::Function, jmax::Int=60)
  initiate!(oldset, z)
  empty!(newset)
  y = T[0.0, 0.0]
  j = 0

  while j < jmax && length(oldset) > 0
    multiply!(oldset, mat)
    @fastmath @inbounds for d in 1:number_of_digits
      @simd for k in 1:length(oldset)
        y[1] = mdigits[1, d] + oldset.container[1, k]
        y[2] = mdigits[2, d] + oldset.container[2, k]

        if in_domain(y)
          append!(newset, y)

          if length(newset) == newset.maxlen #TO FIX
            return T(1)
          end
        end
      end
    end

    j += 1
    temp = oldset
    oldset = newset
    newset = temp
    empty!(newset)
  end

  if j < jmax
    return convert(T, j/jmax)
  else
    return T(1)
  end
end

max_stars_repo_name: kalvotom/julia-examples

max_issues_repo_name: kalvotom/julia-examples

max_forks_repo_name: kalvotom/julia-examples


text:
using OrdinaryDiffEq, Test, Random, LinearAlgebra, SparseArrays

# Parameters
Nc = 22
η = 1.0
κ = 1.0
T = (0.0,100.0)

# Matrix definitions
A = sparse(diagm(1 => sqrt.(Complex[1:Nc;])))
H = η*(A + A') - 1.0im*κ*A'*A

u0 = zeros(ComplexF64, Nc+1)
u0[1] = 1.0
function f_psos(du, u, t, p)
    du .= -1.0im*H*u
end

# Callback
rng = MersenneTwister(rand(UInt))
jumpnorm = Ref(rand(rng))
djumpnorm(x::Vector{ComplexF64}, t, integrator) = norm(x)^2 - (1-jumpnorm[])
function dojump(integrator)
    x = integrator.u
    t = integrator.t

    x .= normalize(A*x)
    jumpnorm[] = rand(rng)
end

cb = ContinuousCallback(djumpnorm,dojump)

prob = ODEProblem{true}(f_psos,u0,T)

sol_tot = []
Ntraj = 100
for i=1:Ntraj
    rng = MersenneTwister(rand(UInt))
    # Tweaking tolerances and dtmax also is not reliable
    sol = solve(prob,DP5(),save_everystep=true,callback=cb,
            abstol=1e-8,reltol=1e-6,dtmax=10)
    push!(sol_tot, sol)
end

# This number has to be η^2/κ^2 in steady-state; all trajectories should converge there
n = [[norm(A * normalize(s.u[j]))^2 for j=1:length(s.t)] for s=sol_tot]

@test all(η^2/κ^2 .≈ [k[end] for k in n])

#=
using Plots
gr()

p1 = plot(sol_tot[1].t, n[1], lw = 2)
for i=2:Ntraj
    plot!(p1,sol_tot[i].t, n[i])
end

p2 = plot(sol_tot[1].t, norm.(sol_tot[1].u).^2)
for i=2:Ntraj
    plot!(p2,sol_tot[i].t, norm.(sol_tot[i].u).^2)
end
=#

using OrdinaryDiffEq, DiffEqCallbacks, Test

# Initial state
u0=[0, -0.25, 0.42081, 0]

function hheom!(du,u,p,t)
    du[1] = u[3]
    du[2] = u[4]
    du[3] = -u[1] - 2u[1]*u[2]
    du[4] = -u[2] - (u[1]^2 - u[2]^2)
    return nothing
end

@inline Vhh(q1, q2) = 1//2 * (q1^2 + q2^2 + 2q1^2 * q2 - 2//3 * q2^3)
@inline Thh(p1, p2) = 1//2 * (p1^2 + p2^2)
@inline Hhh(q1, q2, p1, p2) = Thh(p1, p2) + Vhh(q1, q2)
@inline Hhh(u::AbstractVector) = Hhh(u...)

# Energy
const E = Hhh(u0)

function ghh(resid, u, p)
    resid[1] = Hhh(u[1],u[2],u[3],u[4]) - E
    resid[2:4] .= 0
end

# energy conserving callback:
# important to use save = false, I dont want rescaling points
cb = ManifoldProjection(ghh, nlopts=Dict(:ftol=>1e-13), save = false)

# Callback for Poincare surface of section
function psos_callback(j, direction = +1, offset::Real = 0,
    callback_kwargs = Dict{Symbol, Any}(:abstol=>1e-9))

    # Prepare callback:
    s = sign(direction)
    cond = (u,t,integrator) -> s*(u - offset)
    affect! = (integrator) -> nothing

    cb = DiffEqBase.ContinuousCallback(cond, nothing, affect!; callback_kwargs...,
    save_positions = (true,false), idxs = j)
end

# with this callback, the saved values of variable 1 should be zero
poincarecb = psos_callback(1)

totalcb = CallbackSet(poincarecb, cb)

prob = ODEProblem(hheom!, u0, (0.0, 100.0), callback = totalcb)

extra_kw = Dict(:save_start=>false, :save_end=>false)
DEFAULT_DIFFEQ_KWARGS = Dict{Symbol, Any}(:abstol => 1e-9, :reltol => 1e-9)

sol = solve(prob, Vern9(); extra_kw..., DEFAULT_DIFFEQ_KWARGS..., save_everystep = false)

Es = [Hhh(sol[:, i]) for i in 1:length(sol)]
Eerror = maximum(@. abs(E - Es))

a = sol[1, :]

@test Eerror < 1e-10
for el in a
    @test abs(el) < 1e-10
end


text:
# Animate a 3D Trajectory

getbounds(list) = [minimum(list), maximum(list)]

function animTraj3D(pt::parseTrajectory, lag = 5, land = [0; 0; 0])
    xyzList = splitDimensions(pt.sList)

    x = xyzList[1]
    y = xyzList[2]
    z = xyzList[3]

    x_b = getbounds(x)
    y_b = getbounds(y)
    z_b = getbounds(z)

    cube = getbounds([x_b..., y_b...])

    N = length(x)

    my_anim = @animate for i in 1:(N + lag)
        j = min(i, N)
        plot3d(x[1:j], y[1:j], z[1:j], markershape = :circle,
                        zlabel = "Z (m)", ylabel = "Y (m)", xlabel = "X(m)",
                        legend = :none, zlims = (z_b[1], z_b[2]))
        scatter3d!([land[1]], [land[2]], [land[3]],
                        markershape = :xcross, markercolor = :red,
                        markersize = 8)
        title!("Trajectory over time")

        xlims!(cube...)
        ylims!(cube...)
    end

    return my_anim
end

max_stars_repo_name: danineamati/TrajOptSOCPs


text:
function convert_cells_to_pop(cells, cell_phenotypes, guides)
    cells_to_phenotypes = [DefaultDict{Float64, Int}(0) for _ in 1:length(guides)]

    @inbounds for i in eachindex(cells)
        cells_to_phenotypes[cells[i]][cell_phenotypes[i]] += 1
    end

    cells_to_phenotypes
end

function test_crispri_construction()
    # Test the two different CRISPR types
    lib = Library(CRISPRi())
    setup = FacsScreen()
    guides, guide_freqs_dist = construct_library(setup, lib)
    cells, cell_phenotypes = build_cells(lib.cas9_behavior, guides, guide_freqs_dist, 10^6)

    # In a CRISPRi screen all cells should have the same phenotype
    all(map(length, convert_cells_to_pop(cells, cell_phenotypes, guides)) .== 1)
end
@test test_crispri_construction()

function test_crisprko_construction()
    lib = Library(CRISPRn())
    setup = FacsScreen()
    guides, guide_freqs_dist = construct_library(setup, lib)
    cells, cell_phenotypes = build_cells(lib.cas9_behavior, guides, guide_freqs_dist, 10^7)

    arr = convert_cells_to_pop(cells, cell_phenotypes, guides)
    nonzeros = arr[find(x->length(x) != 1, arr)]
    results = zeros(CRISPRn().knockout_dist.K, length(nonzeros))
    for (idx, guide) in enumerate(nonzeros)
        if -0.0 in keys(guide)
            @assert guide[0.0] == 0
            delete!(guide, 0.0)
            ks = sort(collect(keys(guide)))
            tot = sum(values(guide))
            results[:, idx] = [guide[key]/tot for key in ks]
        else
            @assert guide[0.0] != 0
            ks = sort(collect(keys(guide)), rev=true)
            tot = sum(values(guide))
            results[:, idx] = [guide[key]/tot for key in ks]
        end
    end

    isapprox(mean(results[1, :] ./ results[2, :]), 1, atol=0.25) &&
    isapprox(mean(results[1, :] ./ results[3, :]), 4, atol=0.25)
end

@test test_crisprko_construction()


text:
module Helper

params = Dict(
    :default =>
    Dict(:x0 => 0.1, :dim => 1000, :steps => 10000, :r => 4.0, :epsilon => 0.1),
)

function init(name)
    param = copy(params[name])
    x = zeros(Float64, (param[:dim], param[:steps]))
    x[:, 1] = param[:x0]
    param[:x] = x
    param
end

end

max_stars_repo_name: tkf/comparatist

max_issues_repo_name: tkf/comparatist

max_forks_repo_name: tkf/comparatist


text:
# dependencies:
# - none

module util

function matlab_display(obj)
    class_name = string( typeof(obj) )
    props = propertynames(obj)
    sp = map(string, props)
    L = map(length, sp)
    Lmax = Base.maximum(L)

    props_vec = Vector{String}(undef, length(sp))
    for kk = 1:length(sp)
        props_vec[kk] = "\t" * lpad(sp[kk], Lmax) * ": " * string( getproperty(obj, props[kk]) ) * "\n"
    end

    str = class_name * " with properties:\n"
    str_end = reduce(*, props_vec)
    return str * str_end * "\n"
end

function swap!(array::AbstractArray, i1::Int, i2::Int)
    tmp = array[i1];
    array[i1] = array[i2];
    array[i2] = tmp;
    return array
end

function mid_point(a::T, b::T)::T where {T}
    return (b - a) / 2 + a
end

function mid_point(a::T, b::T)::T where {T<:Integer}
    return floor( (b - a) / 2 + a )
end

function swappop!(array::AbstractArray, idx::Int)
	tmp = array[end]
	array[end] = array[idx]
	array[idx] = tmp
	return pop!(array)
end

function modulo_index(n, L)
	if n <= 0
		n = n % L
		n = L + n
		#  # L in domain 0.. -L+1
	end

	return  (n-1) % L + 1
end

function next_index(index, L)
    # ring indexing
    if index >= L
        index = one(typeof(index))
    else
        index = index + one(typeof(index))
    end
    return index
end

function prev_index(index, L)
    # ring indexing
    if index <= 1
        index = L
    else
        index = index - one(typeof(index))
    end
    return index
end

function is_vector(A)
    return length(size(A)) == 1
end

function is_matrix(A)
    return length(size(A)) == 2
end

function maximum(a::AbstractArray{T}) where {T}
    n = length(a)
    if n < 1
        return (nothing, zero(Int64))
    end

    max_elem = a[1]
    max_idx = one(Int64)

    @inbounds @simd for kk = 2:n
        @inbounds a_ = a[kk]
        if isless(max_elem, a_)
            max_elem = a_
            max_idx = kk
        end
    end

    return (max_elem, max_idx)
end

function maximum(f, a::AbstractArray{T}) where {T}
    n = length(a)
    if n < 1
        return (nothing, zero(Int64))
    end

    max_elem = f(a[1])
    max_idx = one(Int64)

    @inbounds @simd for kk = 2:n
        @inbounds fa = f(a[kk])
        if isless( max_elem, fa )
            max_elem = fa
            max_idx = kk
        end
    end

    return (max_elem, max_idx)
end

function linspace(x_start, x_end, x_length::Int)
    if x_length <= 1
        return [x_end]
    end
    x_delta = (x_end - x_start) / (x_length - 1)
    x_ = collect( x_start:x_delta:x_end )
    return x_
end

function logspace(a, b, n::Int = 50)
	retval = 10 .^ collect( LinRange(a, b, n) )
end

function mesh_grid(x, y)
    X = [i for i in x, j in 1:length(y)]
    Y = [j for i in 1:length(x), j in y]
    return (X, Y)
end

function linear_search(
    x::Vector{T},
    xq::T;
    left_init_index = 1) where {T}

    # it is assumed that the x is sorted in increasing order
    idx = 0 # if zero -> no valid index found

    L = length(x)

    # handle extrapolation
    if xq <= x[1]
        return 0
    elseif xq >= x[L]
        return L+1
    end

    # handle interpolation
    ii = left_init_index
    while ii <= L
        if xq <= x[ii]
            idx = ii-1
            break
        end
        ii += 1
    end

    return idx

end

function binary_search(
    x,
    xq;
    left_init_index = 1,
    right_init_index = length(x))

    # it is assumed that the x is sorted in increasing order
    idx = 0 # if zero -> no valid index found
    # if -1     left extrapolated
    # if L+1    right extrapolated

    L = length(x)

    # handle extrapolation
    if xq < x[1]
        return 0
    elseif xq > x[L]
        return L+1
    end

    left_idx  = left_init_index
    right_idx = right_init_index
    mid_idx   = div(right_idx + left_idx, 2)

    while true

        if xq <= x[mid_idx]
            right_idx = mid_idx
        else # xq > curve.x[mid_idx]
            left_idx = mid_idx
        end

        mid_idx = div(right_idx + left_idx, 2)
        if mid_idx == left_idx
            break
        end
    end
    idx = left_idx

    return idx

end

end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia