text:
## Add code
#push!(LOAD_PATH, "/Users/dchilder/Dropbox/DSGE project/matlab/Huggett/Julia/clusterversion")
#to be able to load module

# A bunch of functions for working with Chebyshev polynomials

function dctTypeI(coeffs)
    # Shameless copy/translation of function described below in package Chebfun, ripped from context
    #COEFFS2VALS   Convert Chebyshev coefficients to values at Chebyshev points
    #of the 2nd kind.
    #   V = COEFFS2VALS(C) returns the values of the polynomial V(i,1) = P(x_i) =
    #   C(1,1)*T_{0}(x_i) + ... + C(N,1)*T_{N-1}(x_i), where the x_i are
    #   2nd-kind Chebyshev nodes.
    #
    #   If the input C is an (N+1)xM matrix then V = COEFFS2VALS(C) returns the
    #   (N+1)xM matrix of values V such that V(i,j) = P_j(x_i) = C(1,j)*T_{0}(x_i)
    #   + C(2,j)*T_{1}(x_i) + ... + C(N,j)*T_{N-1}(x_i).
    #
    # See also VALS2COEFFS, CHEBPTS.

    # Copyright 2017 by The University of Oxford and The Chebfun Developers.
    # See http://www.chebfun.org/ for Chebfun information.

    ################################################################################
    # [Developer Note]: This is equivalent to Discrete Cosine Transform of Type I.
    #
    # [Mathematical reference]: Sections 4.7 and 6.3 Mason & Handscomb, "Chebyshev
    # Polynomials". Chapman & Hall/CRC (2003).
    ################################################################################

    # *Note about symmetries* The code below takes steps to
    # ensure that the following symmetries are enforced:
    # even Chebyshev COEFFS exactly zero ==> VALUES are exactly odd
    # odd Chebychev COEFFS exactly zero ==> VALUES are exactly even
    # These corrections are required because the MATLAB FFT does not
    # guarantee that these symmetries are enforced.

    #Convert to floating point type
    coeffs=float(coeffs)

    # Get the length of the input:
    n = size(coeffs, 1);

    # Trivial case (constant or empty):
    if ( n <= 1 )
        fvalues = coeffs;
        return
    end

    # check for symmetry
    isEvencol = maximum(abs.(coeffs[2:2:end,:]),1).== 0;
    isOddcol = maximum(abs.(coeffs[1:2:end,:]),1).== 0;

    # Scale them by 1/2:
    coeffs[2:n-1,:] = coeffs[2:n-1,:]/2;

    # Mirror the coefficients (to fake a DCT using an FFT):
    tmp = [ coeffs ; coeffs[n-1:-1:2,:] ];

    if ( isreal(coeffs) )
        # Real-valued case:
        fvalues = real(fft(tmp,1));
    elseif ( isreal(1i*coeffs) )
        # Imaginary-valued case:
        fvalues = 1i*real(fft(imag(tmp),1));
    else
        # General case:
        fvalues = fft(tmp,1);
    end

    # Flip and truncate:
    fvalues = fvalues[n:-1:1,:];

    # enforce symmetry
    for j=1:size(fvalues,2)
        if (isEvencol[j])
            fvalues[:,j] = (fvalues[:,j]+flipdim(fvalues[:,j],1))/2;
        elseif (isOddcol[j])
            fvalues[:,j] = (fvalues[:,j]-flipdim(fvalues[:,j],1))/2;
        end
    end

    return fvalues

end #dctTypeI

function c2v1(coeffs)
    #c2v1 is a shameless copy of chebfun function described below to get values at type I Chebyshev points

    #COEFFS2VALS   Convert Chebyshev coefficients to values at Chebyshev points
    #of the 1st kind.
    #   V = COEFFS2VALS(C) returns the values of the polynomial V(i,1) = P(x_i) =
    #   C(1,1)*T_{N-1}(x_i) + C(2,1)*T_{N-2}(x_i) + ... + C(N,1), where the x_i are
    #   1st-kind Chebyshev nodes.
    #
    #   If the input C is an (N+1)xM matrix then V = COEFFS2VALS(C) returns the
    #   (N+1)xM matrix of values V such that V(i,j) = P_j(x_i) = C(1,j)*T_{N-1}(x_i)
    #   + C(2,j)*T_{N-2}(x_i) + ... + C(N,j)
    #
    # See also VALS2COEFFS, CHEBPTS.

    # Copyright 2017 by The University of Oxford and The Chebfun Developers.
    # See http://www.chebfun.org/ for Chebfun information.

    #Convert to floating point type
    coeffs=float(coeffs)

    # Get the length of the input:
    n = size(coeffs, 1);
    m = size(coeffs,2);

    # Trivial case (constant):
    if ( n <= 1 )
        return fvalues = coeffs;
    end

    # check for symmetry
    isEvencol = maximum(abs.(coeffs[2:2:end,:]),1).== 0;
    isOddcol = maximum(abs.(coeffs[1:2:end,:]),1).== 0;

    #Precompute weight vector
    w = (exp.(-1im*[0:2*n-1...]*pi/(2*n))/2);
    w[1] = 2*w[1];
    w[n+1] = 0;
    w[n+2:end] = -w[n+2:end];

    # Mirror the values for FFT:
    c_mirror = [coeffs ; ones(1, m) ; coeffs[end:-1:2,:]];

    # Apply the weight vector:
    c_weight = broadcast(*, c_mirror, w);
    fvaltemp = fft(c_weight,1);

    # Truncate and flip the order:
    fvalues = fvaltemp[n:-1:1,:];

    # Post-process:
    if ( isreal(coeffs) )
        # Real-valued case:
        fvalues = real(fvalues);
    elseif ( isreal(1im*coeffs) )
        # Imaginary-valued case:
        fvalues = 1im*imag(fvalues);
    end

    # enforce symmetry
    for j=1:size(fvalues,2)
        if (isEvencol[j])
            fvalues[:,j] = (fvalues[:,j]+flipdim(fvalues[:,j],1))/2;
        elseif (isOddcol[j])
            fvalues[:,j] = (fvalues[:,j]-flipdim(fvalues[:,j],1))/2;
        end
    end

    return fvalues

end #c2v1

function ccquadwts(n,chebkind)
    #CCQUADWTS   Quadrature weights for Chebyshev points of 2nd kind.
    #   CCQUADWTS(N) returns the N weights for Clenshaw-Curtis quadrature on 2nd-kind
    #   Chebyshev points.
    #
    # See also CHEBPTS, BARYWTS.

    # Copyright 2017 by The University of Oxford and The Chebfun Developers.
    # See http://www.chebfun.org/ for Chebfun information.

    ################################################################################
    # DEVELOPER NOTE:
    # We use a variant of Waldvogel's algorithm [1], due to Nick Hale. (See below)
    # We note this is similar to Greg Von Winkel's approach, which can be found on
    # the MathWorks File Exchange.
    #
    # Let $f(x) = \sum_{k=0}^nc_kT_k(x)$, then\vspace*{-3pt} }
    #   I(f) = v.'*c
    # where
    #   v = \int_{-1}^1T_k(x)dx = { 2/(1-k^2) : k even
    #                             { 0         : k odd
    #     = v'*inv(TT)*f(x) where TT_{j,k} = T_k(x_j)
    #     = (inv(TT)'*v)'*f(x)
    # Therefore
    #   I(f) = w.'f(x) => w = inv(TT).'*v;
    # Here inv(TT).' = inv(TT) is an inverse discrete cosine transform of Type I.
    #
    # Furthermore, since odd entries in v are zero, can compute via FFT without
    # doubling up from N to 2N (though we still need to double up from N/2 to N to
    # facilitate the use of ifft).
    #
    # References:
    #   [1] Joerg Waldvogel, "Fast construction of the Fejer and Clenshaw-Curtis
    #       quadrature rules", BIT Numerical Mathematics 46 (2006), pp 195-202.
    #   [2] Greg von Winckel, "Fast Clenshaw-Curtis Quadrature",
    #       http://www.mathworks.com/matlabcentral/fileexchange/6911, (2005)
    ################################################################################

    if ( n == 0 )                      # Special case (no points!)
        w = [];
    elseif ( n == 1 )                  # Special case (single point)
        w = 2;
    elseif (chebkind==2)                                 # General case, Polynomials of 2nd kind
        tvec = 1 .- [2.0:2.0:(n-1)...].^2
        tk=zeros(size(tvec,1)+1,1);
        tk[1,1]=1;
        tk[2:end,1]=tvec
        c = 2 ./ tk';  # Exact integrals of T_k (even)
        cvec=c[floor(Int,n/2):-1:2]';
        ck=zeros(1,size(c,2)+size(cvec,2))
        ck[1,1:size(c,2)] = c;
        ck[1,size(c,2)+1:end] = cvec; # Mirror for DCT via FFT
        intw = real(ifft(ck,2));                   # Interior weights
        zz=intw[1]/2;                      # Boundary weights
        w=zeros(1,n);
        w[1:size(intw,2)]=intw;
        w[1] = zz;
        w[n] = zz;
    elseif (chebkind==1)                                 # General case, Polynomials of 1st kind
        tvec = 1 .- [2.0:2.0:(n-1)...].^2
        tk=zeros(size(tvec,1)+1,1);
        tk[1,1]=1;
        tk[2:end,1]=tvec
        m = 2 ./ tk';  # Exact integrals of T_k (even)
        # Mirror the vector for the use of ifft:
        if mod(n,2)==1 #n is odd
            cvec=-m[floor(Int,(n+1)/2):-1:2]';
            ck=zeros(1,size(m,2)+size(cvec,2))
            ck[1,1:size(m,2)] = m;
            ck[1,size(m,2)+1:end] = cvec; # Mirror for DCT via FFT
        else #n is even
            cvec=-m[floor(Int,n/2):-1:2]';
            ck=zeros(1,size(m,2)+size(cvec,2)+1)
            ck[1,1:size(m,2)] = m;
            ck[1,size(m,2)+2:end] = cvec;
        end
        v = exp.(1im*[0:n-1...]'*pi/n);      # weight (rotation) vector
        c = ck .* v;                      # Apply the weight vector
        w = real(ifft(c,2));             # Call ifft
    else
        println("Error: Chebkind must be 1 or 2")
    end

    return w

end #ccquadwts

function chebpts(n,lo,hi,chebkind)

    #Function to return grid points and associated Clenshaw Curtis quadrature weights for
    #Chebyshev polynomials
    if ( n == 0 )     # Special case (no points)
        x = [];
    elseif n==1
        x = 0;
    else
        if chebkind==2
            #second kind (the ones with points at the boundary)
            m = n - 1;
            # Points for default domain of [-1,1]
            x = copy(transpose(sin.(pi*(-m:2:m)/(2*m))));  # (Use of sine enforces symmetry.)
        elseif chebkind==1
            x = copy(transpose(sin.(pi*((-n+1:2:n-1)/(2*n))))); # (Use of sine enforces symmetry.)
        else
            println("Error: Chebtype Must be 1 or 2")
        end
    end

    xr=((hi .- lo)/2) * x .+ (hi .+ lo)/2; # Rescaled to new grid

    xwts=((hi-lo)/2)*ccquadwts(n,chebkind); #Rescale Clenshaw Curtis Quadrature Weights

    (cgrid,cwts)=[xr,xwts];

end #chebpts

# cheb2leg
# A function to convert Chebyshev coefficients into orthonormal Legendre polynomial coefficients
# Based on cheb2leg.m in Chebfun for Matlab

function cheb2leg(c_cheb,normalize)
    # Input is array of Chebyshev coefficients,
    # set normalize to true to output coefficients on orthonormal legendre basis

    chebsize=size(c_cheb)
    if chebsize[1]==1
        c_leg=c_cheb
    elseif  chebsize[1] > 1
        #     # Use direct approach:
        c_leg = cheb2leg_direct(c_cheb,normalize);
    end

    ## Matlab version gives option to call approximate algorithm when n large
    ## Not included here: maybe useful feature to add
    # else
    #     # Use fast algorithm:
    #     c_leg = cheb2leg_fast( c_cheb, normalize );
    # end

return c_leg

end #cheb2leg

function cheb2leg_direct(c_cheb, normalize)
    #CHEB2LEG_DIRECT   Convert Cheb to Leg coeffs using the 3-term recurrence.
    N = size(c_cheb,1)                     # Number of columns.
    m = size(c_cheb,2)                   # Number of rows.
    N = N - 1;                          # Degree of polynomial.
    x = cos.(.5*pi*transpose((0:2*N))/N);         # 2*N+1 Chebyshev grid (reversed order).
    temparg=zeros(2*N+1,m);
    temparg[1:N+1,1:m]=c_cheb;
    ff = dctTypeI(temparg);                           #Ordering and normalization differ from Matlab version
    f=flipdim(ff,1);                                   # Values on 2*N+1 Chebyshev grid.
    w = copy(transpose(ccquadwts(2*N+1,2)));     # Clenshaw-Curtis quadrature weights.
    Pm2 = 1;
    Pm1 = x;                   # Initialize.
    L = zeros(2*N+1, N+1);              # Vandermonde matrix.
    L[:,1] = 1+0*x;          #P_0
    L[:,2] = x;              #P_1.
    for k = 1:N-1 # Recurrence relation:
        P = (2-1/(k+1))*Pm1 .* x - (1-1/(k+1))*Pm2;
        Pm2 = Pm1;
        Pm1 = P;
        L[:,2+k] = P;
    end
    scale = (2*[0:N...]+1)/2;            # Scaling in coefficients [NIST, (18.3.1)].
    c_leg = broadcast(*, transpose(L)*(broadcast(*, f ,w)), scale); # Legendre coeffs.

    # Normalize:
    if ( normalize )
        c_leg  = broadcast(*, c_leg, 1 ./ sqrt.([0:N...]+.5) );
    end

    return c_leg

end #cheb2leg_direct

function bary(x, fvals, xk)
    #BARY   Barycentric interpolation formula.

    # xk are (2nd kind Chebyshev) grid points, fvals are values at grid points,
    # x are points where you want to interpolate function (using Chebyshev polynomial)

    #   BARY(X, FVALS, XK, VK) uses the 2nd form barycentric formula with weights VK
    #   to evaluate an interpolant of the data {XK, FVALS(:,k)} at the points X.
    #   Note that XK and VK should be column vectors, and FVALS, XK, and VK should
    #   have the same length.
    #
    #   BARY(X, FVALS) assumes XK are the 2nd-kind Chebyshev points and VK are the
    #   corresponding barycentric weights.
    #
    #   If size(FVALS, 2) > 1 then BARY(X, FVALS) returns values in the form
    #   [F_1(X), F_2(X), ...], where size(F_k(X)) = size(X) for each k.
    #

    # Copyright 2017 by The University of Oxford and The Chebfun Developers.
    # See http://www.chebfun.org/ for Chebfun information.

    # Parse inputs:
    x=x[:,:]; # Represent as matrix

    n = size(fvals,1);
    m = size(fvals,2);
    sizex = size(x);
    ndimsx = ndims(x);

    if ( (m > 1) && (ndimsx > 2) )
        error("Too many Dimensions")
    end

    # Omit default behavior: always ask for grids and weights
    # # Default to Chebyshev nodes and barycentric weights:
    # if ( nargin < 3 )
    #     xk = chebtech2.chebpts(n);
    # end
#
# if ( nargin < 4 )
#     vk = chebtech2.barywts(n);
# end

#Use Barycentric Weights corresponding to second kind Chebyshev

vk = barywts(n);

# Check that input is a column vector:
if ( (ndimsx > 2) || (sizex[2] > 1) )
    x = x[:];
end

# The function is a constant.
if ( n == 1 )
    return fx = repeat(fvals, size(x,1), 1);
end

# The function is NaN.
if ( any(isnan.(fvals)) )
    return    fx = NaN(size(x,1), m);
end

# The main loop:
if ( prod(size(x)) < 4*size(xk,1) )  # Loop over evaluation points
    # Note: The value "4" here was detemined experimentally.

    # Initialise return value:
    fx = zeros(size(x, 1), m);

    # Loop:
    for j = 1:prod(size(x))
        xvar = vk ./ (x[j]*ones(size(xk)) - xk);
        fx[j,:] = (xvar'*fvals) ./ sum(xvar);
    end
else                            # Loop over barycentric nodes
    # Initialise:
    numerat = zeros(size(x, 1), m);
    denom = numerat;

    # Loop:
    for j = 1:size(xk,1)
        tmp = (vk[j] ./ (x - xk[j]));
        numerat = numerat + broadcast(*, tmp, fvals[j,:]);
        denom = broadcast(+, denom, tmp);
    end
    fx = numerat ./ denom;
end

# Try to clean up NaNs:

nanind=find(isnan,fx[:,1]);

for i = 1:size(nanind,1)
    k=nanind[i]; #Place where there is a NaN
    index = findfirst(xk,x[k]);    # Find the corresponding node
    if  index!=0
        fx[k,:] = fvals[index,:];   # Set entry/entries to the correct value
    end
end

# Reshape the output if possible:
if ( (m == 1) && ( (ndimsx > 2) || (sizex[2] > 1) ) )
    fx = reshape(fx, sizex);
elseif ( (m > 1) && ( (ndimsx == 2) || (sizex[2] > 1) ) )
    fx = reshape(fx, sizex[1], m*prod(size(x))/sizex[1]);
end

return fx

end #bary

function barywts(n)
    #returns Barycentric weights for interpolation using Chebyshev polynomials of the second kind

    #BARYWTS   Barycentric weights for Chebyshev points of 2nd kind.
    #   BARYWTS(N) returns the N barycentric weights for polynomial interpolation on
    #   a Chebyshev grid of the 2nd kind. The weights are normalised so that they
    #   have infinity norm equal to 1 and the final entry is positive.
    #
    # See also BARY, CHEBPTS.

    # Copyright 2017 by The University of Oxford and The Chebfun Developers.
    # See http://www.chebfun.org/ for Chebfun information.

    ################################################################################
    # See Thm. 5.2 of Trefethen, Approximation Theory and Approximation Practice,
    # SIAM, 2013 for more information.
    ################################################################################

    if ( n == 0 )                      # Special case (no points)
        v = [];
    elseif ( n == 1 )                  # Special case (single point)
        v = 1;
    else                               # General case
        v = [ones(n-1,1) ; .5];        # Note v(end) is positive.
        v[end-1:-2:1] = -1;
        v[1] = .5*v[1];
    end

    return v

end #barywts

# ChebDerivMat
# Create matrix corresponding to collocation representation of differentiation
# based on recurrence relation for derivatives of Chebyshev polynomials at second kind chebyshev points (Tn(x))
# Returns matrix with proper sclaing, but

function ChebDerivMat(n,lo,hi)

    D = zeros(n,n) #Initialize recurrence relation
    # D is coefficients to coefficients transform
    if n>1
        D[1,2]=1.0 #Polynomial of order 1 has derivative
        if n>2
            D[2,3]=4
            if n>3
                for j=4:n
                    D[j-1,j]=2*(j-1)
                    D[:,j]+=((j-1)/(j-3))*D[:,j-2]
                end
            end
        end
    end

    Minv=dctTypeI(eye(n))
    M=pinv(Minv)

    Deriv=(2/(hi-lo))*Minv*D*M #Grid values to grid values transform

    return Deriv, D

end #ChebDerivMat


text:
module Day09

using AdventOfCode2019

function day09(input::String = readInput(joinpath(@__DIR__, "..", "data", "day09.txt")))
    program = parse.(Int, split(input, ","))
    channelIn = Channel{Int}(1)
    put!(channelIn, Int(1))
    out1 = AdventOfCode2019.IntCode.run_program!(program, channelIn, nothing, nothing)
    close(channelIn)

    channelIn = Channel{Int}(2)
    put!(channelIn, Int(2))
    out2 = AdventOfCode2019.IntCode.run_program!(program, channelIn, nothing, nothing)
    close(channelIn)
    return [out1[end], out2[end]]
end

end # module

ext: jl


text:
function _prepare_approximation(nodes_b::Matrix{T},
                                kernel::RK
                               ) where {T <: AbstractFloat, RK <: ReproducingKernel_0}

    spline = _prepare_approximation(nothing, nodes_b, kernel)
    return spline
end

function _prepare_approximation(nodes::Union{Matrix{T}, Nothing},
                                nodes_b::Matrix{T},
                                kernel::RK
                               ) where {T <: AbstractFloat, RK <: ReproducingKernel_0}

     if isnothing(nodes)
         n_1 = 0
         n = size(nodes_b, 1)
     else
        n = size(nodes, 1)
        n_1 = size(nodes, 2)
        if n != size(nodes_b, 1)
            error("Matrices 'nodes' and 'nodes_b' have mismatched dimensions.")
        end
     end

     n_b = size(nodes_b, 1)
     n_1_b = size(nodes_b, 2)

     min_bound = Vector{T}(undef, n)
     compression::T = T(0.)
     maxx::T = T(0.)
     @inbounds for i = 1:n
         if !isnothing(nodes)
             min_bound[i] = nodes[i,1]
             maxx = nodes[i,1]
         else
             min_bound[i] = nodes_b[i,1]
             maxx = nodes_b[i,1]
         end
         for j = 1:n_1
             min_bound[i] = min(min_bound[i], nodes[i,j])
             maxx = max(maxx, nodes[i,j])
         end
         for j = 1:n_1_b
             min_bound[i] = min(min_bound[i], nodes_b[i,j])
             maxx = max(maxx, nodes_b[i,j])
         end
         compression = max(compression, maxx - min_bound[i])
     end

     if compression <= eps(T(1.0))
         error("Cannot prepare the spline: `nodes` data are not correct.")
     end

     if !isnothing(nodes)
         t_nodes = similar(nodes)
     else
         t_nodes = nothing
     end
     @inbounds for j = 1:n_1
         for i = 1:n
             t_nodes[i,j] = (nodes[i,j] - min_bound[i]) / compression
         end
     end
     t_nodes_b = similar(nodes_b)
     @inbounds for j = 1:n_1_b
         for i = 1:n
             t_nodes_b[i,j] = (nodes_b[i,j] - min_bound[i]) / compression
         end
     end

     if !isnothing(nodes)
         t_nodes_all = [t_nodes t_nodes_b]
     else
         t_nodes_all = t_nodes_b
     end

     if T(kernel.ε) == T(0.0)
         ε = _estimate_ε(t_nodes_all)
         if isa(kernel, RK_H0)
             kernel = RK_H0(ε)
         elseif isa(kernel, RK_H1)
             ε *= T(1.5)
             kernel = RK_H1(ε)
         elseif isa(kernel, RK_H2)
             ε *= T(2.0)
             kernel = RK_H2(ε)
         else
             error("Incorrect `kernel` type.")
         end
     end

     gram = _gram(t_nodes_all, kernel)
     chol = nothing
     try
         chol = cholesky(gram)
     catch
         error("Cannot prepare the spline: Gram matrix is degenerate.")
     end

     cond = _estimate_cond(gram, chol)

     spline = NormalSpline(kernel,
                           compression,
                           t_nodes,
                           t_nodes_b,
                           nothing,
                           nothing,
                           nothing,
                           nothing,
                           nothing,
                           nothing,
                           nothing,
                           nothing,
                           nothing,
                           nothing,
                           min_bound,
                           gram,
                           chol,
                           nothing,
                           nothing,
                           cond,
                           -1
                          )
     return spline
end

function _construct_approximation(spline::NormalSpline{T, RK},
                                  values_lb::Vector{T},
                                  values_ub::Vector{T},
                                  maxiter::Int,
                                  ftol::T,
                                  cleanup::Bool = false
                                 ) where {T <: AbstractFloat, RK <: ReproducingKernel_0}

    if !isnothing(spline._nodes)
         error("Spline object is not correctly prepared ('spline._nodes' is not empty).")
    end

    spline, nit_done = _construct_approximation(spline, T[], values_lb, values_ub, maxiter, ftol, cleanup)
    return spline, nit_done
end

function _construct_approximation(spline::NormalSpline{T, RK},
                                  values::Vector{T},
                                  values_lb::Vector{T},
                                  values_ub::Vector{T},
                                  maxiter::Int,
                                  ftol::T,
                                  cleanup::Bool = false
                                 ) where {T <: AbstractFloat, RK <: ReproducingKernel_0}

    if isnothing(spline._nodes) && isnothing(spline._nodes_b)
        error("Spline object is not correctly prepared ('spline._nodes' is empty and 'spline._nodes_b' is empty).")
    end
    if isnothing(spline._nodes_b)
        error("Spline object is not correctly prepared ('spline._nodes_b' is empty).")
    end

    m1 = 0
    if !isnothing(spline._nodes)
        m1 = size(spline._nodes, 2)
    end
    m2 = size(spline._nodes_b, 2)

    if length(values) != m1
        error("Number of 'values' does not correspond to the number of interpolating nodes ('spline._nodes').")
    end
    if length(values_lb) != m2
        error("Number of 'values_lb' does not correspond to the number of approximating nodes ('spline._nodes_b').")
    end
    if length(values_ub) != m2
        error("Number of 'values_ub' does not correspond to the number of approximating nodes.")
    end

    @inbounds for i = 1:m2
        if values_ub[i] < values_lb[i]
            error("Incorrect bounds: 'values_ub[$i]' is less than 'values_lb[$i]'.")
        end
        if (values_ub[i] - values_lb[i]) <= sqrt(eps(T(1.0)))
            error("Incorrect bounds: 'values_ub[$i]' is almost equal to 'values_lb[$i]'.")
        end
     end

    if isnothing(spline._chol)
        error("Gram matrix was not factorized.")
    end

    umax = T(0.)
    @inbounds for i = 1:m2
        val = abs(values_lb[i])
        if val != Inf && val > umax
            umax = val
        end
        val = abs(values_ub[i])
        if val != Inf && val > umax
            umax = val
        end
    end

    ineq_active = zeros(Int, m2)

    values_b  = Vector{T}(undef, m2)
    incorrect_bounds::Bool = false
    @inbounds for i = 1:m2
        if abs(values_lb[i]) == Inf && abs(values_ub[i]) == Inf
            incorrect_bounds = true
            break
        end
        if abs(values_lb[i]) != Inf && abs(values_ub[i]) != Inf
            values_b[i] = (values_lb[i] + values_ub[i]) / T(2.0)
        end
        if abs(values_lb[i]) == Inf
            values_b[i] = values_ub[i] - umax
        end
        if abs(values_ub[i]) == Inf
            values_b[i] = values_lb[i] + umax
        end
    end
    if incorrect_bounds
        error("Incorrect_bounds. Both 'values_ub[i]' and 'values_lb[i]' are set to Inf.")
    end

    if m1 > 0
        values_all = [values; values_b]
    else
        values_all = values_b
    end

    mu = Vector{T}(undef, size(spline._gram, 1))
    ldiv!(mu, spline._chol, values_all)

    spline = NormalSpline(spline._kernel,
                          spline._compression,
                          spline._nodes,
                          spline._nodes_b,
                          values,
                          values_lb,
                          values_ub,
                          nothing,
                          nothing,
                          nothing,
                          nothing,
                          nothing,
                          nothing,
                          nothing,
                          spline._min_bound,
                          spline._gram,
                          spline._chol,
                          mu,
                          nothing,
                          spline._cond,
                          -1
                         )

    precision = T(1.e-10) # TODO change

    spline, nit_done = _qp(spline, ineq_active, maxiter, ftol, precision, cleanup)
    return spline, nit_done
end

###################

function _evaluate(spline::NormalSpline{T, RK}, points::Matrix{T},
                  ) where {T <: AbstractFloat, RK <: ReproducingKernel_0}
    if !isnothing(spline._nodes_b) || !isnothing(spline._d_nodes_b)
        _evaluate_approximation(spline, points)
    else
        _evaluate_interpolation(spline, points)
    end
end

function _evaluate_gradient(spline::NormalSpline{T, RK},
                            point::Vector{T}
                           ) where {T <: AbstractFloat, RK <: ReproducingKernel_0}
    if !isnothing(spline._nodes_b) || !isnothing(spline._d_nodes_b)
        _evaluate_approximation_gradient(spline, point)
    else
        _evaluate_interpolation_gradient(spline, point)
    end
end

function _evaluate_approximation_gradient(spline::NormalSpline{T, RK},
                                          point::Vector{T}
                                         ) where {T <: AbstractFloat, RK <: ReproducingKernel_0}
    # TODO - implement
    error("_evaluate_approximation_gradient: Not implemented.")
end

function _evaluate_approximation(spline::NormalSpline{T, RK}, points::Matrix{T},
                                ) where {T <: AbstractFloat, RK <: ReproducingKernel_0}
    if isnothing(spline)
        error("Spline was not prepared.")
    end
    if isnothing(spline._mu)
        error("Spline was not constructed.")
    end

    if isnothing(spline._ier < 0)
        error("Spline coefficients are not calculated.")
    end

    if size(points, 1) != size(spline._nodes_b, 1)
        if size(points, 1) == 1 && size(points, 2) > 1
            error("Incorrect first dimension of the `points` parameter (use 'evaluate_at' function for evaluating the spline at one point).")
        else
            error("Incorrect first dimension of the `points` parameter (the spline was built in the space of different dimension).")
        end
    end

    if isnothing(spline._nodes)
       n_1 = 0
    else
       n_1 = size(spline._nodes, 2)
    end

    n = size(spline._nodes_b, 1)
    n_1_b = size(spline._nodes_b, 2)
    m = size(points, 2)

    pts = similar(points)
    @inbounds for j = 1:m
        for i = 1:n
            pts[i,j] = (points[i,j] - spline._min_bound[i]) / spline._compression
        end
    end

    spline_values = Vector{T}(undef, m)
    h_values = Vector{T}(undef, n_1)

    if spline._ier != 0 #  optimal solution not found

        h_values_b = Vector{T}(undef, n_1_b)
        @inbounds for p = 1:m
            for i = 1:n_1_b
                h_values_b[i] = _rk(spline._kernel, pts[:,p], spline._nodes_b[:,i])
            end
            if n_1 > 0
                for i = 1:n_1
                    h_values[i] = _rk(spline._kernel, pts[:,p], spline._nodes[:,i])
                end
                spline_values[p] = sum(spline._mu .* [h_values; h_values_b; -h_values_b])
            else
                spline_values[p] = sum(spline._mu .* [h_values_b; -h_values_b])
            end
        end
    else # optimal solution found
        nak = length(spline._active)
        h_values_b = Vector{T}(undef, nak)
        mu_b = Vector{T}(undef, nak)
        @inbounds for p = 1:m
            spline_values[p] = T(0.)
            if n_1 > 0
                for i = 1:n_1
                    h_values[i] = _rk(spline._kernel, pts[:,p], spline._nodes[:,i])
                end
                spline_values[p] += sum(spline._mu[1:n_1] .* h_values)
            end

            for i = 1:nak
                ii = spline._active[i]
                if ii > 0
                    mu_b[i] = spline._mu[ii+n_1]
                    h_values_b[i] = _rk(spline._kernel, pts[:,p], spline._nodes_b[:,ii])
                    s = 0
                else
                    ii = -ii
                    mu_b[i] = spline._mu[ii+n_1+n_1_b]
                    h_values_b[i] = -_rk(spline._kernel, pts[:,p], spline._nodes_b[:,ii])
                    s = 0
                end
            end
            spline_values[p] += sum(mu_b .* h_values_b)
        end
    end # if spline._ier != 0

    return spline_values
end


text:
using SimJulia

# minimalistic simjulia evactuation example
# - random ambulances are generated, they only do one trip
# - random patients generated
# - simulation stops when all processes are terminated
const hospitals = ["St. Luc", "Erasme", "UZ Jette"]

struct Victim
    name::String
    state::Int
end
Base.show(io::IO, v::Victim) = print(io, "victim $(v.name) in state $(v.state)")

mutable struct Ambulance
    number::Int
    hospital::String
    traveltime::Int
    p::Process
    function Ambulance(sim::Environment, number::Int, hospital::String, traveltime::Int, victims::Array{Victim, 1})
        a = new()
        a.number = number
        a.hospital = hospital
        a.traveltime = traveltime
        a.p = @process ambulance(sim, a, victims)
    end
end

Base.show(io::IO, a::Ambulance) = print(io, "ambulance $(a.number) from $(a.hospital)")

@resumable function ambulance(sim::Environment, a::Ambulance, victims::Array{Victim, 1})
    @info "$(now(sim)) - $(a) moving"
    @yield timeout(sim, a.traveltime)
    @info "$(now(sim)) - $(a) on site"
    @yield timeout(sim, rand(5:10))
    # pick random victim
    v = popat!(victims, rand(1:length(victims)))
    @info "$(now(sim)) - $(a) leaving site with $(v)"
    @yield timeout(sim, a.traveltime)
    @info "$(now(sim)) - $(a) at hospital with $(v)"
end

function disaster(n_v=10, n_a=5)
    @info "\n\t\tEvaction demo\n\n"
    sim = Simulation()
    victims = [Victim("person $(i)", rand(1:4)) for i in 1:n_v]
    ambs = [Ambulance(sim, i, rand(hospitals), 10 + rand(1:10), victims) for i in 1:n_a]
    run(sim)
    @info "$(length(victims))/$(n_v) victims remaining on site"
    println()
end

disaster()

# minimalistic object transfering
# - only a single jump is done
# - not all steps are present
# - chutes are transferred from one store to another
# - chutes is warehouse "age" in their ready state

mutable struct Chute
    id::Int
    age::Int
    ready::Bool
end

Base.show(io::IO, c::Chute) = print(io, """chute $(c.id) (age $(c.age), $(c.ready ? "" : "not") ready)""")
age!(c::Chute) = c.age += 1

@resumable function fillmag(sim::Environment, mag::Store, n::Int=5)
    chutes = [Chute(i, 0, true) for i in 1:n]
    for chute in chutes
        @yield put(mag, chute)
    end
    @info "$(now(sim)) - current warehouse: $(mag.items)"
end

# cycle for single jump
@resumable function jump(sim::Environment, n::Int,  mag::Store, folders::Store)
    @info "$(now(sim)) - dropping $(n) paratroopers"
    for _ in 1:n
        chute = @yield get(mag)
        chute.ready = false
        @info "got $(chute)"
        @process chutecycle(sim, chute, folders)
    end
end

# aging process
@resumable function aging(sim::Environment, mag::Store)
    while true
        @yield timeout(sim, 1)
        age!.(mag.items)
    end
end

@resumable function chutecycle(sim::Environment, c::Chute, folders::Store)
    @info "$(now(sim)) - $(c) is drying"
    @yield timeout(sim, rand(1:3))
    @info "$(now(sim)) - $(c) is dry, going to fold"
    # store chute with the folders
    @yield put(folders, c)
    @info "$(now(sim)) - waiting list folding: $(folders.items)"
end

@resumable function fold(sim::Environment, mag::Store, folders::Store)
    while true
        @yield timeout(sim, 10)
        for _ in 1:37 # daily folds per person
            if length(folders.items) > 0
                chute = @yield get(folders)
                chute.ready = true
                @info "$(now(sim)) - folded $(chute)"
                @yield put(mag, chute)
            end
        end
    end
end

function airborne(nchutes=3)
    @info "\n\t\tParatrooper demo\n\n"
    sim = Simulation()
    mag = Store{Chute}(sim)
    folders = Store{Chute}(sim)
    # fill storage with chutes
    @process fillmag(sim, mag, nchutes)
    # add single jump to simulation
    @process jump(sim, 2, mag, folders)
    # add folding process
    @process fold(sim, mag, folders)
    # add aging process
    @process aging(sim, mag)
    run(sim, 20)
    @info "$(length(mag.items)) chutes in warehouse, $(length(folders.items)) waiting to be folded"
    @info mag.items

    println()
end

airborne()

########################################################
# Simple and empty process

@resumable function filler(sim::Environment, s::Store)
    i = 0
    while true
        @yield put(s, i)
        @info "$(now(sim)) - current store: $(s.items)"
        i += 1
        @yield timeout(sim, 1)
    end
end

@resumable function taker(sim::Environment, s::Store)
    while true
        res = @yield get(s)
        @info "$(now(sim)) - got $(res)"
        @yield timeout(sim, rand(1:4))
    end
end

function testfun()
    @info "\n\t\tStore usage demo\n\n"
    sim = Simulation()
    s = Store{Int}(sim)
    @process filler(sim, s)
    @process taker(sim, s)
    run(sim, 10)
end

testfun()

# worker process with prioties
# - each job has an ID and a priority
# - items are removed from store according to their priority
# - task can be split over multiple days

mutable struct job
    number::Int           # job identification
    prio::Int             # priority (0: normal, 1: high, 2: super high (i.e. unfinished))
    duration # standard time delta when working with datetime
    start       # moment task started
    stop        # moment when finished
    finished::Bool        # task terminated
    # only intialise required fields
    function job(number; prio::Int=0, duration)
        j = new()
        j.number = number
        j.prio = prio
        j.duration = duration
        return j
    end
end

Base.show(io::IO, j::job) = print(io, "job $(j.number) (prio: $(j.prio), duration: $(j.duration))")

@resumable function work(sim::Environment, tasks::Store{job}, completed::Store{job})
    while true
        # check if (super) high prio jobs exist and get one
        if length([t for t in tasks.items if t.prio == 2]) > 0
            @info "$(now(sim)) - VERY HIGH prio task(s) found"
            current = @yield get(tasks, x-> x.prio == 2)
        elseif length([t for t in tasks.items if t.prio == 1]) > 0
            @info "$(now(sim)) - HIGH prio task(s) found"
            current = @yield get(tasks, x-> x.prio == 1)
        elseif length([t for t in tasks.items if t.prio == 0]) > 0
            @info "$(now(sim)) - default prio task(s) found"
            current = @yield get(tasks, x-> x.prio == 0)
        else
            @info "$(now(sim)) - Nothing to do!"
            current = nothing
        end

        # next day if no work
        if isnothing(current)
            @yield timeout(sim, 24)
        # proces the task
        else
            @info "$(now(sim)) - working on task $(current)"
            current.start = now(sim)
            if (current.start + current.duration) % 24 > 16
                @info "$(now(sim)) - not enough time for task $(current)"
                # calculate remaining time
                dt =  current.duration - ((current.start + current.duration) % 24 - 16)
                # add remaing part of task on next day with super high prio
                current.duration -= dt
                current.prio = 2
                @yield put(tasks, current)
                # move on to end of day
                @yield timeout(sim, 16 - now(sim))
                @info "$(now(sim)) - end of day!"
                @yield timeout(sim, 16)
            else
                @yield timeout(sim, current.duration)
                current.stop = now(sim)
                current.finished = true
                @info "$(now(sim)) - task $(current) is done!"
                @yield put(completed, current)
            end
        end
    end
end

function jobman()
    @info "\n\t\tJobman demo\n\n"
    sim = Simulation(8)
    jobs = Store{job}(sim)
    completedjobs = Store{job}(sim)
    # add some jobs
    for i in 1:10
        push!(jobs.items, job(i, prio=rand([0,1]), duration = 10))
    end
    @process work(sim, jobs, completedjobs)
    run(sim, 50)
    @info "$(length(completedjobs.items)) completed tasks: \n\n$(completedjobs.items)\n\n"
    @info "$(length(jobs.items)) planned tasks: \n\n$(jobs.items)\n"
end

jobman()


text:
export
    device, setdevice,
    synchronize

const contexts = Ptr{Void}[]

function initctx()
    cuInit(0)
    for dev = 0:devcount()-1
        ref = Ptr{Void}[0]
        cuCtxCreate(ref, 0, dev)
        push!(contexts, ref[1])
    end
    setdevice(0)
end

function device()
    ref = Ptr{Void}[0]
    cuCtxGetDevice(ref)
    Int(ref[1])
end

setdevice(dev::Int) = cuCtxSetCurrent(contexts[dev+1])

synchronize() = cuCtxSynchronize()


text:
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the compressible Euler equations
gamma = 1.4
equations = CompressibleEulerEquations1D(gamma)

initial_condition = initial_condition_density_wave

surface_flux = flux_central
polydeg = 5
solver = DGSEM(polydeg, surface_flux)

coordinates_min = (-1.0,)
coordinates_max = ( 1.0,)
mesh = TreeMesh(coordinates_min, coordinates_max,
                initial_refinement_level=2,
                n_cells_max=30_000)

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)

###############################################################################
# ODE solvers, callbacks etc.

tspan = (0.0, 2.0)
ode = semidiscretize(semi, tspan)

###############################################################################
# run the simulation

sol = solve(ode, BS3(), save_everystep=false);


text:
import LDLFactorizations
const LDLF = LDLFactorizations

"""
    LDLFact <: LSBackend

Use LDLFactorizations backend.

Options available:
* Any numerical type `T`
* Augmented system with LDLᵀ factorization
"""
struct LDLFact <: LSBackend end

# ==============================================================================
#   LDLFLinearSolver
# ==============================================================================

"""
    LDLFLinearSolver{Tv}

Linear solver for the 2x2 augmented system with ``A`` sparse.

Uses LDLFactorizations.jl to compute an LDLᵀ factorization of the quasi-definite augmented system.
"""
mutable struct LDLFLinearSolver{Tv<:Real} <: AbstractLinearSolver{Tv}
    m::Int  # Number of rows
    n::Int  # Number of columns

    # TODO: allow for user-provided ordering,
    #   and add flag (default to false) to know whether user ordering should be used

    # Problem data
    A::SparseMatrixCSC{Tv, Int}
    θ::Vector{Tv}
    regP::Vector{Tv}  # primal regularization
    regD::Vector{Tv}  # dual regularization

    # Factorization
    F::LDLF.LDLFactorization{Tv}

    # TODO: constructor with initial memory allocation
    function LDLFLinearSolver(A::SparseMatrixCSC{Tv, Int}) where{Tv<:Real}
        m, n = size(A)
        θ = ones(Tv, n)

        S = [
            spdiagm(0 => -θ)  A';
            A spdiagm(0 => ones(m))
        ]

        # TODO: PSD-ness checks
        F = LDLF.ldl(S)
        return new{Tv}(m, n, A, θ, ones(Tv, n), ones(Tv, m), F)
    end

end

AbstractLinearSolver(
    ::LDLFact,
    ::AugmentedSystem,
    A::AbstractMatrix{Tv}
) where{Tv<:Real} = LDLFLinearSolver(sparse(A))

backend(::LDLFLinearSolver) = "LDLFactorizations.jl"
linear_system(::LDLFLinearSolver) = "Augmented system"

"""
    update_linear_solver!(ls, θ, regP, regD)

Update LDLᵀ factorization of the augmented system.

Update diagonal scaling ``\\theta``, primal-dual regularizations, and re-compute
    the factorization.
Throws a `PosDefException` if matrix is not quasi-definite.
"""
function update_linear_solver!(
    ls::LDLFLinearSolver{Tv},
    θ::AbstractVector{Tv},
    regP::AbstractVector{Tv},
    regD::AbstractVector{Tv}
) where{Tv<:Real}
    # Sanity checks
    length(θ)  == ls.n || throw(DimensionMismatch(
        "θ has length $(length(θ)) but linear solver is for n=$(ls.n)."
    ))
    length(regP) == ls.n || throw(DimensionMismatch(
        "regP has length $(length(regP)) but linear solver has n=$(ls.n)"
    ))
    length(regD) == ls.m || throw(DimensionMismatch(
        "regD has length $(length(regD)) but linear solver has m=$(ls.m)"
    ))

    # Update diagonal scaling
    ls.θ .= θ
    # Update regularizers
    ls.regP .= regP
    ls.regD .= regD

    # Re-compute factorization
    # TODO: Keep S in memory, only change diagonal
    S = [
        spdiagm(0 => -ls.θ .- regP)  ls.A';
        ls.A spdiagm(0 => regD)
    ]

    # TODO: PSD-ness checks
    try
        ls.F = LDLF.ldl(S)
    catch err
        isa(err, LDLF.SQDException) && throw(PosDefException(-1))
        rethrow(err)
    end

    return nothing
end

"""
    solve_augmented_system!(dx, dy, ls, ξp, ξd)

Solve the augmented system, overwriting `dx, dy` with the result.
"""
function solve_augmented_system!(
    dx::Vector{Tv}, dy::Vector{Tv},
    ls::LDLFLinearSolver{Tv},
    ξp::Vector{Tv}, ξd::Vector{Tv}
) where{Tv<:Real}
    m, n = ls.m, ls.n

    # Set-up right-hand side
    ξ = [ξd; ξp]

    # Solve augmented system
    d = ls.F \ ξ

    # Recover dx, dy
    @views dx .= d[1:n]
    @views dy .= d[(n+1):(m+n)]

    # TODO: Iterative refinement
    return nothing
end


text:
@with_kw struct recipe
	epochs::Array{Int64,1} = [1000]
	N::Array{Int64,1} = [8000]

	θ::Float64=2000
	h::Float64=0.5

	s₋::Float64=-457.0
	s₊::Float64=500.0

	dfe::String="point"
	param_one::Float64=1.0
	param_two::Float64=1.0
	s::Array{Float64,1}=-[0.0]
	s_mult::Array{Float64,1}=[1.0]
	prob::Array{Float64,1}=[0.0]

	n_anc::Int64=1000
	burnin_period::Bool=false

	relax::Bool=false
	epoch_relaxation::Array{Bool}=fill(false,length(epochs))
	s_relaxation::Float64=0.0
	s_relaxation_threshold::Float64=0.0

	F::Array{Float64,1}=zeros(length(N))
	seed::Int64=rand(1:10^8)

	trajectories::Array{Int64,1} = Int64[]
	trajectories_output::String = ""

	@assert length(N)==length(epochs)  "N and epochs must be equal in length";
	@assert length(s)==length(prob)    "s and probs must be equal in length";
end

@with_kw mutable struct prf_output
	fixations::Int64 = 0;
	loss::Int64=0;
	total_mut::Int64 = 0;
	count_mut::Int64=0;
end

@with_kw mutable struct mutation
	# The frequency of the mutation
	frequency::Float64=0;
	# Count of the mutation
	count_samp::Int64=1;
	# Selection coefficient
	s::Float64=0.0;
	# Dominance factor
	h::Float64=0.5;
	# The generation that mutation has arise
	age::Int64=1;
	# Number of mutation. Tag to avoid similar nodes
	num::Int64=1;
	# Mutation type. 1=neutral; 2=deleterious; 3=strong advantegeous; 4=weakly advantegeous
	type::Int8=0
end

function relax_selection!(mutation_list::LinkedList{mutation}, new_s, sel_threshold::Float64, relax_type::Int64)

	if(!isempty(mutation_list))

		# Index to delete node
		node=mutation_list.node.next;
		while (node.next != mutation_list.node.next)
			item        = node.data;

			freq            = item.frequency;
			s               = item.s;
			h               = item.h;
			freq            = freq_inbreed(s,freq,h);
			num       = item.num;

			if ((s/2) <= sel_threshold && relax_tye == 0)
				item.s = new_s * 2;
			elseif ((s/2) <= sel_threshold && relax_tye == 1)
				item.s = s * new_s * 2;
			end
			node  = node.next
		end
	end
end

function freq_inbreed(sel::Float64,freq::Float64,h::Float64,F::Float64)
	#The fitnesses of the genotypes are A1A1 = 1; A1A2 = 1-h*(2*s); A2A2 = 1-(2*s), where A2 is the derived allele. Under any conditions that the program is run, the value of s for any segregating site cannot exceed 0.5.
	# parameters struct assume positive values of s to concieve positive selection. Changing to negative value
	s::Float64=-sel
	return (((1.0 - s) * (freq * freq + F * freq * (1.0 - freq))) + ((1.0 - h * s) * freq * (1.0 - freq) * (1.0 - F))) / (((1.0 - freq) * (1.0 - freq) + (1.0 - freq) * F * freq) + ((1.0 - h * s) * 2.0 * freq * (1.0 - freq) * (1.0 - F)) + ((1.0 - s) * (freq * freq + F * freq * (1.0 - freq))))
end

function add_mutation!(mutation_list::LinkedList{mutation},r::Ptr{gsl_rng},N::Float64,h::Float64,s::Float64,θ::Float64,freq::Float64,dfe::String,param_one::Float64,param_two::Float64,s_mult::Float64,n_anc::Int64,age::Int64,trajectories::Array{Int64,1},relax::Bool=false)

	num_mut::Float64 = ran_poisson(r,θ/2.0)
	count_mut::Int64 = length(mutation_list)

	for x::Int64=1:num_mut

		if dfe == "point"
			s = s;
		elseif dfe == "gamma"
			##/Use this for gamma distribution in Boyko 2008:
			gamma = ran_gamma(r, param_one, param_two * s_mult);
			# note, this is scale for a Nanc=1000, using the boyko params
			s = - gamma / (n_anc * 2);
		end

        count_mut += 1;
        mut        = mutation(frequency = freq, h = h, s = s*2, count_samp = 0.0, age = age, num = count_mut, type=1)
		push!(mutation_list,mut);

		if(!isempty(trajectories) && count_mut in trajectories)
			i = findfirst(isequal(count_mut),trajectories);
			@printf(io, "%d\t%lf\n", trajectories[i], (1.0/N));
		end
	end
end

function drift_sel!(mutation_list::LinkedList{mutation},r::Ptr{gsl_rng},N::Float64,F::Float64,trajectories::Array{Int64,1})

	l = 0;
	f = 0;

	if(!isempty(mutation_list))

		# Index to delete node
		node=mutation_list.node.next;
		while (node.next != mutation_list.node.next)
			item        = node.data;

            freq            = item.frequency;
            s::Float64      = item.s;
            h::Float64      = item.h;
            freq::Float64   = freq_inbreed(s,freq,h,F);

            num       = item.num;

            #count::Int64    = rand(Binomial(N,freq));
            count::Int64    = ran_binomial(r,freq,N)
            item.count_samp = count;
            freq            = count / N;
            item.frequency  = freq;

			if (freq > 0.0 && freq < 1.0)
				if(!isempty(trajectories) && num in trajectories)
					i = findfirst(isequal(num),trajectories);
					@printf(io, "%d\t%lf\n", trajectories[i], freq);
				end
				# Move node
				node = node.next
			else
				if(freq == 0)
					l +=1;
				elseif(freq == 1)
					f +=1;
				end
				# Delete node
				deleteat!(mutation_list,node);
				# Move node
				node = node.next
			end
		end
		return(l,f)
	else
		return(l,f)
	end
end

function burnin(param::recipe,r::Ptr{gsl_rng})

	@unpack N,θ,s,h,dfe,param_one,param_two,s_mult,n_anc,trajectories,relax = param;

	n_size = N[1]-1;
	s_size = N[1]-1;
	sel    = s[1];

	theoretical_sfs     = zeros(n_size,2);
	number_of_mutations = 0;

	mutation_list = LinkedList{mutation}()

	f(x=Float64,j=Int64) = f_of_q_lambda(x,j,n_size,s_size,sel,θ*2)
	@time @inbounds for j::Int64=2:n_size
		m = LinkedList{mutation}()

		res,err = quadgk(x->f(x,j),0,1)
		age = n_size * 10;
		add_mutation!(mutation_list,r,Float64(n_size),h,s,θ,j/n_size,dfe,param_one,param_two,sel,n_anc,age,trajectories,relax);
		# append!(mutation_list,m);
	end

	return(mutation_list);
end

function f_of_q_lambda(x::Float64,j::Int64,N_burnin::Int64,sample_size::Int64,point_sel::Float64,theta::Float64)

	gamma = N_burnin * point_sel;

	if (abs(gamma) > 1.0e-7)

		sfs_function_term_one = (1-exp(-2*N_burnin*(-point_sel)*(1-x)))/(1-exp(-2*N_burnin*(-point_sel)));
		sfs_function_term_two::Float64 = (2/(x*(1-x)));

		binomial_value::Float64 = ran_binomial_pdf(j,x,sample_size);

		f = theta/2 * sfs_function_term_one * sfs_function_term_two * binomial_value;

	else
		f = theta/2 * 2/x * ran_binomial_pdf(j,x,sample_size);
	end
	return  f;
end

function sfs(mutation_list::LinkedList{mutation},r::Ptr{gsl_rng},sample_size::Int64)

	if(!isempty(mutation_list))

		mut  = zeros(length(mutation_list));
		#=neut = zeros(length(mutation_list));
		del  = zeros(length(mutation_list));
		adv  = zeros(length(mutation_list));=#

		i    = 1;
		# Index to delete node
		node = mutation_list.node.next;

		while (node.next != mutation_list.node.next)
			item        = node.data;

			freq        = item.frequency;
			s           = item.s;

			h           = item.h;
			#samp_count = rand(Binomial(sample_size,freq));
			samp_count = ran_binomial(r,freq,sample_size);
			samp_freq = samp_count/sample_size;
			if (samp_freq > 0 && samp_freq < 1.0)
				mut[i] = samp_freq
				#=if(item.type == 1)
					neut[i] = samp_freq
				elseif(item.type == 2)
					del[i] = samp_freq
				elseif(item.type == 3)
					adv[i] = samp_freq
				end=#
			end
			node = node.next;
			i+=1
		end

		mut     = mut[mut .!= 0]

		b        = collect(1/sample_size:1/sample_size:1)
		v(x,b=b) = searchsortedfirst.(Ref(b), x)
		n    = v(mut)

		out      = zeros((sample_size-1,2))
		out[:,1] = collect(1:(sample_size-1))

		@inbounds for i::Int64 in out[:,1]
			out[i,2] = length(n[n .== i])
		end

		out[:,1] = round.(out[:,1]/sample_size,digits=3);

		return out
	else
		out = zeros((sample_size-1,3));
		out[:,1] = round.(out[:,1]/sample_size,digits=3);
		return out
	end
end

function simulate(param::recipe, sample_size::Int64)

	@unpack epochs,N,θ,h,s₋,s₊,dfe,param_one,param_two,s,s_mult,prob,n_anc,burnin_period,relax,epoch_relaxation,s_relaxation,s_relaxation_threshold,F,trajectories,trajectories_output,seed = param;

	if !isempty(trajectories)
		@assert length(unique(trajectories .> 0)) == 1  "ID index must be greater than 0";
		io = open(trajectories_output,"w+")
	end

	# set seed before start
	T = gsl_rng_default;
	r = rng_alloc(T);

	rng_set(r, seed)

	epochs = SVector{length(epochs)}(epochs)
	s_mult = SVector{length(s_mult)}(s_mult)
	s = SVector{length(s)}(s)
	F = SVector{length(F)}(F)

	events = length(epochs);

	if(length(s)==1)
		s      = s[1];
		s_mult = s_mult[1]
	else
		s = s;
		s_mult = s_mult;
	end

	@printf "Demographic History (%i epochs)\n\n" events;
	for e=1:events
		@printf "Ne = %i\tgenerations = %i\tF = %lf\n" N[e] epochs[e] param.F[e];
	end

	mutation_list = LinkedList{mutation}();

	if !burnin_period
		# mutation_list = burnin(mutation_list,param);
		mutation_list = burnin(param,r);
	else
		mutation_list = LinkedList{mutation}();
	end;

	l = 0;
	f = 0;
	N_F_1=0;
	@printf "\n"
	@inbounds for e=1:events
		# Inbreeding Ne

		N_F = N[e] / (1.0 + F[e]);
		@printf "Currently in epoch = %i ; Mutations before epoch's beginning = %i\n"  e length(mutation_list);

		if epoch_relaxation[e]
			printf("Relaxation in Epoch %i\n", e);
			relax_selection!(mutation_list, s_relaxation, s_relaxation_threshold, relaxation_type);
		end

		@inbounds for g=1:epochs[e]
			# Mutation age
			if e == 1
				age = g;
			else
				# Add age to previous epochs generations
				age = g + epochs[e];
			end
			x, y    = drift_sel!(mutation_list,r,N_F,F[e],trajectories);
			add_mutation!(mutation_list,r,N_F,h,s,θ,1.0/N[e],dfe,param_one,param_two,s_mult,n_anc,age,trajectories,relax);
			l = l + x;
			f = f + y;
		end

		if (e < events)
			N_F_1 = N[e + 1] / (1.0 + param.F[e + 1]);
			θ = θ * N_F_1 / N_F;
		end
	end

	if @isdefined io
		close(io)
	end

	out = sfs(mutation_list,r,sample_size);

	return(out,hcat(l,f))
end

function simulate_batch(param::recipe,sample_size::Int64,replicas::Int64)

	mutation_list = LinkedList{mutation}();

	n_params = [deepcopy(param) for i::Int64=1:replicas];
	n_mutations_list = [deepcopy(mutation_list) for i::Int64=1:replicas];
	n_sample_size = [sample_size for i::Int64=1:replicas]

	mut = pmapbatch(simulate,n_params,n_sample_size);
	s = sum(x->x[1][:,2],mut);
	f = sum(x->x[2],mut);

	s = hcat(round.(collect(1:(sample_size-1))/sample_size,digits=3),s);

	return s,f
end

function two_epochs(ne_prior::Union{Int64,Nothing},θ_prior::Union{Float64,Nothing},sample_size::Int64,replicas::Int64;expansion=false,bottleneck=false)

	if(isnothing(θ_prior))
		θ_prior = rand(50:1000,replicas)
	else
		θ_prior = fill(θ_prior,replicas)
	end

	if(isnothing(ne_prior))
		ne_prior = rand(1000:10000,replicas)
	else
		ne_prior = fill(ne_prior,replicas)
	end

	n_sample_size = fill(sample_size,replicas)
	n_params = recipe[];

	if expansion
		ne_rand = rand(2:50,replicas)
	elseif bottleneck
		ne_rand = rand(0.01:0.01:1,replicas)
	else
		ne_rand = rand(0.01:0.01:50,replicas);
	end

	ne_two      =  Int64.(trunc.(ne_rand.*ne_prior,digits=0));
	epoch_two   = rand(100:2000,replicas)

	for i::Int64=1:replicas
		push!(n_params,recipe(N=[ne_prior[i],ne_two[i]],epochs=[1,epoch_two[i]],θ=θ_prior[i],burnin_period=false))
	end

	mut = pmapbatch(simulate,n_params,n_sample_size);

	f(x) = map(y->y[1][:,2],x)
	ps   = hcat(f(mut)...)

	return(ps,hcat(ne_prior,ne_two,epoch_two))
end

function three_epochs(ne_prior::Union{Int64,Nothing},θ_prior::Union{Float64,Nothing},sample_size::Int64,replicas::Int64;expansion=false,bottleneck=false)

	if(isnothing(θ_prior))
		θ_prior = rand(50:1000,replicas)
	else
		θ_prior = fill(θ_prior,replicas)
	end

	if(isnothing(ne_prior))
		ne_prior = rand(1000:10000,replicas)
	else
		ne_prior = fill(ne_prior,replicas)
	end

	if expansion
		ne_rand       = rand(2:50,replicas)
		ne_rand_three = rand(2:50,replicas)
	elseif bottleneck
		ne_rand       = rand(0.01:0.01:1,replicas)
		ne_rand_three = rand(0.01:0.01:1,replicas)
	else
		ne_rand       = rand(0.01:0.01:50,replicas);
		ne_rand_three = rand(0.01:0.01:50,replicas);
	end

    ne_two        = Int64.(trunc.(ne_rand.*ne_prior,digits=0));
    ne_three      = Int64.(trunc.(ne_rand_three.*ne_two,digits=0));

    epoch_two     = rand(100:1000,replicas)
    epoch_three   = rand(100:1000,replicas)

    n_sample_size = fill(sample_size,replicas)
    n_params      = recipe[];
	for i::Int64=1:replicas
		push!(n_params,recipe(N=[ne_prior[i],ne_two[i],ne_three[i]],epochs=[ne_prior[i]*10,epoch_two[i],epoch_three[i]],θ=θ_prior[i]))
	end

	mut = pmapbatch(simulate,n_params,n_sample_size);

	f(x) = map(y->y[1][:,2],x)
	ps   = hcat(f(mut)...)

	return(ps)
end

function exponential(ne_prior::Union{Int64,Nothing},θ_prior::Union{Float64,Nothing},rate::Union{Float64,Nothing},sample_size::Int64,replicas::Int64;expansion=false,bottleneck=false)

	if(isnothing(θ_prior))
		θ_prior = rand(50:1000,replicas)
	else
		θ_prior = fill(θ_prior,replicas)
	end

	if(isnothing(ne_prior))
		ne_prior = rand(1000:10000,replicas)
	else
		ne_prior = fill(ne_prior,replicas)
	end

	if(isnothing(rate))
		rate = rand(0.00005:0.0001:0.02,replicas)
	else
		rate = fill(rate,replicas)
	end

	n_sample_size = fill(sample_size,replicas)
	n_params = recipe[];

	if expansion
		ne_rand = rand(2:10,replicas)
	elseif bottleneck
		ne_rand = rand(0.01:0.01:1,replicas)
	else
		ne_rand = rand(0.01:0.01:50,replicas);
	end

	ne_two      =  trunc.(Int64,ne_rand.*ne_prior);
	epoch_two   = rand(100:10000,replicas)
	epoch_exp   = rand(100:1000,replicas)

	g(p::Int64,r::Float64,t::Int64) = trunc(Int64,p*(1 + r)^t)
	for i::Int64=1:replicas
		tmp_ne    = g.(ne_two[i],rate[i],1:epoch_exp[i])
		tmp_epoch = fill(1,epoch_exp[i])

		push!(n_params,recipe(N=vcat(ne_prior[i],ne_two[i],tmp_ne),epochs=vcat(1,epoch_two[i],tmp_epoch),θ=θ_prior[i],burnin_period=false))
	end

	mut = pmapbatch(simulate,n_params,n_sample_size);

	f(x) = map(y->y[1][:,2],x)
	ps   = hcat(f(mut)...)

	return(ps,hcat(ne_prior,ne_two,epoch_two,epoch_exp,rate))
end


text:
#------------------
# aggregate

"""
    aggregate(object, grid [, fun])

Aggregate `object` (either an `Assemblage` or `Locations` type) to `grid`.
If `object` is an `Assemblage{PointData}` this will grid all points and return an
`Assemblage{GridData}`. `grid` can be a `GridTopology` or a single `Integer`
signifying the aggregation factor for already gridded data, the cellsize for
point data. `fun` is an optional function specifying how to lump occurrences. If
not specified the default function is `any` for Boolean Assemblages and `sum` for
Integer ones.
"""
aggregate(lo::SELocations, factor) = Locations{GridData}(aggregate(lo.coords, factor))
aggregate(gr::GridTopology, factor) = _grid_from_factor(gr, factor)
aggregate(gr::SEGrid, factor::Integer) = (g = _grid_from_factor(gr.grid, factor); aggregate(gr, g))
aggregate(gr::SEGrid, newgrid::GridTopology) = (ind = _apply_grid(gr, newgrid); GridData(sortslices(unique(ind, dims = 1), dims = 1), newgrid))

function aggregate(asm::SEAssemblage{D, T, P}, factor::Union{Integer, Tuple{Integer, Integer}}, fun = _default_fun(D); xmin = nothing, ymin = nothing) where D where T where P <: SELocations

    gt = _grid_from_factor(asm.site.coords, factor; xmin = xmin, ymin = ymin)
    aggregate(asm, gt, fun)
end

function aggregate(asm::SEAssemblage{D}, gt::GridTopology, fun = _default_fun(D)) where D where T where P <: SELocations

    # create the new grid
    tmpgrid = _apply_grid(asm.site.coords, gt)
    tmpsites = sortslices(unique(tmpgrid, dims = 1), dims = 1)

    # fill a new occurrence matrix
    retmat = spzeros(D, nspecies(asm), size(tmpsites, 1))

    # get a dictionary of the indices of each new group
    groupinds = Dict{typeof(first(eachrow(tmpgrid))), Vector{Int}}()
    for (i, row) in enumerate(eachrow(tmpgrid))
        if haskey(groupinds, row)
            push!(groupinds[row], i)
        else
            groupinds[row] = [i]
        end
    end

    for newcell in axes(tmpsites, 1)
        inds = groupinds[tmpsites[newcell, :]]
        if fun === Base.any
            retmat[nzrows(view(occurrences(asm), :, inds)), newcell] .= true
        elseif fun === Base.sum
            retmat[:, newcell] .= vec(colsum(view(occurrences(asm), :, inds)))
        else
            retmat[:, newcell] .= vec(mapslices(fun, occurrences(asm)[:, inds], dims = 2))
        end
    end

    #build the basic objects
    newsitenames = string.(["agg_"], 1:size(retmat, 2))
    sit = Locations{GridData}(GridData(tmpsites, gt), DataFrame(sites = newsitenames))
    retcommat = ComMatrix{D}(retmat, speciesnames(asm), newsitenames)
    sp = SpeciesData(retcommat, traits(asm))
    Assemblage(sit,sp)
end

#-----------
# `aggregate` helper functions

_lowerleft(gr) = (xmin(gr), ymin(gr)) .- (0.5 .* cellsize(gr))

function _apply_grid(gr::SEGrid, newgrid::GridTopology)
    factor = [cellsize(newgrid)...] ./ cellsize(gr)
    shift = (_lowerleft(gr) .- _lowerleft(newgrid)) .% factor ./cellsize(gr)
    newind = ceil.(Int, (shift' .+ gr.indices) ./ factor')
    newind
end
_apply_grid(pt::SEPoints, newgrid::GridTopology) = getindices(coordinates(pt), newgrid)

function _range_from_factor(mi, ma, factor; inmin = nothing)
    newmin = isnothing(inmin) ? (floor(mi / factor) + 0.5) * factor : inmin
    newmax = isnothing(inmin) ? (ceil(ma/factor) - 0.5) * factor : floor((ma - inmin) / factor) * factor + inmin
    range(newmin, stop = newmax, step = factor)
end

function _grid_from_factor(pt::SEPoints, factor; xmin = nothing, ymin = nothing)
    (xmi, xma), (ymi, yma) = mapslices(extrema, coordinates(pt), dims = 1)
    GridTopology(_range_from_factor(xmi, xma, first(factor); inmin = xmin), _range_from_factor(ymi, yma, last(factor); inmin = ymin))
end
_grid_from_factor(gd::SEGrid, factor; xmin = nothing, ymin = nothing) = _grid_from_factor(gd.grid, factor; xmin = xmin, ymin = ymin)
_grid_from_factor(gd::GridTopology, factor; xmin = nothing, ymin = nothing) = GridTopology(_range_from_factor(extrema(gd.xs)..., first(factor); inmin = xmin), _range_from_factor(extrema(gd.ys)..., last(factor); inmin = ymin))

_default_fun(::Type{Bool}) = Base.any
_default_fun(::Type{Integer}) = Base.sum
_default_fun(::Any) = error("Default aggregation functions are only defined for Assemblage{Bool} (presence-absence) and Assemblage{Int}")


text:
function calculate(f, g)
    local a = 0;
    for i = 2:length(f)
        a = a + f[i] * g[i];
    end
    return a;
end

function zeckendorf(j, f, g, y)
    if calculate(f, g) == y
        return (g, true);
    else
        for i = j:length(f)
            g[i] = 1;
            if zeckendorf(i + 2, f, g, y)[2]
                return (g, true);
            end
            g[i] = 0;
        end
    end
    return (g, false);
end

function printer(a)
    deleteat!(a, 1);
    if a[length(a)] == 0
        deleteat!(a, length(a));
    end
    println(join(a[end:-1:1], ""));
end

global x = parse(Int64, readline());
global f = Int64[1, 1];
while f[length(f)] < x
    resize!(f, length(f) + 1);
    f[length(f)] = f[length(f) - 1] + f[length(f) - 2];
end
g = zeros(Int64, length(f));
a = zeckendorf(2, f, g, x);
b = zeckendorf(3, f, g, x);
if a[2]
    printer(a[1]);
else
    printer(b[2]);
end

max_stars_repo_name: MediocreBoris/gitgud

max_issues_repo_name: MediocreBoris/gitgud

max_forks_repo_name: MediocreBoris/gitgud


text:
module PlotlyJSSchemaDocsGenerator

using Base.Markdown: MD
using JSON

# methods to re-construct a plot from JSON
_symbol_dict(x) = x
_symbol_dict(d::Associative) =
    Dict{Symbol,Any}([(Symbol(k), _symbol_dict(v)) for (k, v) in d])

struct SchemaAttribute
    description::Nullable{MD}
    valtype::Nullable{String}
    flags::Nullable{Vector{Any}}
    vals::Nullable{Vector{Any}}
    children::Nullable{Dict{Symbol,SchemaAttribute}}
end

function SchemaAttribute(d::Associative)
    role = pop!(d, :role, "")

    if role == "object" || any(_->isa(_, Associative), values(d))
        # description and valtype ire empty, but children is not
        _desc = pop!(d, :description, "")
        desc = isempty(_desc) ? Nullable{MD}() : Nullable{MD}(MD(_desc))
        valtype = Nullable{String}()
        filter!(d) do k, v
            !(k in (:_deprecated, :tracerefminus)) && !(startswith(string(k), "_"))
        end
        kids = Dict{Symbol,SchemaAttribute}()
        for (a_k, v) in d
            isa(v, Associative) || continue
            kids[a_k] = SchemaAttribute(v)
        end
        children = Nullable{Dict{Symbol,SchemaAttribute}}(kids)
    else
        valtype = Nullable{String}(get(d, :valType, Nullable{String}()))
        desc = Nullable{MD}(
            Base.Markdown.parse(get(d, :description, ""))
        )

        # children is none
        children = Nullable{Dict{Symbol,SchemaAttribute}}()
    end

    flags = valtype == "flaglist" ? Nullable{Vector{Any}}(d["flags"]) :
                                    Nullable{Vector{Any}}()

    vals = valtype == "enumerated" ? Nullable{Vector{Any}}(d["values"]) :
                                       Nullable{Vector{Any}}()

    SchemaAttribute(desc, valtype, flags, vals, children)
end

struct TraceSchema
    name::Symbol
    description::Nullable{MD}
    attributes::Dict{Symbol,SchemaAttribute}
end

function TraceSchema(nm::Symbol, d::Associative, attrs_key=:attributes)
    _attrs = filter!((k, v) -> k != :uid && k != :type, d[attrs_key])
    attrs = Dict{Symbol,SchemaAttribute}()
    for (k, v) in _attrs
        isa(v, Associative) || continue
        attrs[k] = SchemaAttribute(v)
    end
    desc = get(d, :description, "")
    description = isempty(desc) ? Nullable{MD}(MD()) :
                                  Nullable{MD}(Base.Markdown.parse(desc))
    TraceSchema(nm, description, attrs)
end

struct Schema
    traces::Dict{Symbol,TraceSchema}
    layout::TraceSchema

    function Schema()
        _path = joinpath(dirname(@__FILE__), "plotschema.json")
        schema = _symbol_dict(JSON.parse(readstring(_path)))

        traces = Dict{Symbol,TraceSchema}()
        for (k, v) in schema[:schema][:traces]
            traces[k] = TraceSchema(k, v)
        end

        layout = TraceSchema(:layout, schema[:schema][:layout], :layoutAttributes)
        new(traces, layout)
    end
end

doc_html(parent::Symbol, name::Symbol, x::Union{SchemaAttribute,TraceSchema}) =
    takebuf_string(doc_html!(IOBuffer(), parent, name, x))

function doc_html!(buf::IO, parent::Symbol, name::Symbol, sa::SchemaAttribute)
    data_parent = "$(parent)_attributes"
    id = "$(parent)_$(name)"
    print(buf,"
")
    print(buf, """

         ")
     print(buf, name)
     print(buf, """

     """)
     println(buf, "
")

     # if we have a description or children, we need a panel-body
     has_body = !isnull(sa.description) || !isnull(sa.description)

     # add in description if we have one
     if has_body
         print(buf, "
")
     end

     if !isnull(sa.description)
         println(buf, Base.Markdown.html(get(sa.description)))
     end

     if !isnull(sa.children)
         new_parent = Symbol(string(parent), "_", string(name))
         # need to add panel-group
         print(buf, "
")

         # recursively add in children, if any
         kids = get(sa.children)
         for kid_name in sort!(collect(keys(kids)))
             doc_html!(buf, new_parent, kid_name, kids[kid_name])
         end

         # close panel group
         println(buf, "
 ")
     end

     # close body div
     has_body && println(buf, "
  ")

     # close outer divs
     println(buf, "
 ")  #
     println(buf, "
 ")
     buf
end

function doc_html!(buf::IO, data_parent::Symbol, name::Symbol, ts::TraceSchema)
    id = "$(data_parent)_$(name)"
    print(buf, """

       ")
    println(buf, name, "

")

    println(buf, "
")
    println(buf, "
")
    if !isnull(ts.description)
        println(buf, Base.Markdown.html(get(ts.description)))
    end

    print(buf, "
")

    # now add all attributes
    kids = sort!(collect(keys(ts.attributes)))
    for kid_name in kids
        doc_html!(buf, name, kid_name, ts.attributes[kid_name])
    end

    println(buf, "
 ")
    println(buf, "
 ")
    println(buf, "
 ")
    println(buf, "
 ")
    return buf
end

doc_html(s::Schema) = takebuf_string(doc_html!(IOBuffer(), s))

function doc_html!(buf::IO, s::Schema)
    println(buf,
    """<!DOCTYPE html>

      <meta name="viewport" content="width=device-width, initial-scale=1">
      <link rel="stylesheet" href="http://maxcdn.bootstrapcdn.com/bootstrap/3.3.7/css/bootstrap.min.css">

    """)

    println(buf,
    """

Trace attributes

    """)

    # add traces
    for kid_name in sort!(collect(keys(s.traces)))
        doc_html!(buf, :traces, kid_name, s.traces[kid_name])
    end
    # doc_html!(buf, :traces, :scatter, s.traces[:scatter])

    println(buf, "
 ")
    println(buf, "
 ")

    # now do layout
    println(buf,
    """

Layout attributes

    """)
    doc_html!(buf, :layout, :layout, s.layout)

    println(buf, "
 ")
    println(buf, "
 ")
    println(buf, "")
    println(buf, "")
    buf
end

open(joinpath(dirname(@__FILE__), "schema.html"), "w") do f
  doc_html!(f, Schema())
end

end  # module


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

## core string functions ##

endof(s::AbstractString) = error("you must implement endof(", typeof(s), ")")
next(s::AbstractString, i::Int) = error("you must implement next(", typeof(s), ",Int)")
next(s::DirectIndexString, i::Int) = (s[i],i+1)
next(s::AbstractString, i::Integer) = next(s,Int(i))

string() = ""
string(s::AbstractString) = s

"""
    String(s::AbstractString)

Convert a string to a contiguous byte array representation encoded as UTF-8 bytes.
This representation is often appropriate for passing strings to C.
"""
String(s::AbstractString) = print_to_string(s)

convert(::Type{Vector{UInt8}}, s::AbstractString) = convert(Vector{UInt8}, String(s))
convert(::Type{Array{UInt8}}, s::AbstractString) = convert(Vector{UInt8}, s)
convert(::Type{String}, s::AbstractString) = String(s)
convert(::Type{Vector{Char}}, s::AbstractString) = collect(s)
convert(::Type{Symbol}, s::AbstractString) = Symbol(s)
convert(::Type{String}, s::Symbol) = unsafe_string(Cstring(s))

## generic supplied functions ##

start(s::AbstractString) = 1
done(s::AbstractString,i) = (i > endof(s))
getindex(s::AbstractString, i::Int) = next(s,i)[1]
getindex(s::AbstractString, i::Integer) = s[Int(i)]
getindex(s::AbstractString, i::Colon) = s
getindex{T<:Integer}(s::AbstractString, r::UnitRange{T}) = s[Int(first(r)):Int(last(r))]
# TODO: handle other ranges with stride ±1 specially?
getindex{T<:Integer}(s::AbstractString, v::AbstractVector{T}) =
    sprint(length(v), io->(for i in v; write(io,s[i]) end))
getindex(s::AbstractString, v::AbstractVector{Bool}) =
    throw(ArgumentError("logical indexing not supported for strings"))

Symbol(s::AbstractString) = Symbol(String(s))

"""
    sizeof(s::AbstractString)

The number of bytes in string `s`.

```jldoctest
julia> sizeof("❤")
3
```
"""
sizeof(s::AbstractString) = error("type $(typeof(s)) has no canonical binary representation")

eltype{T<:AbstractString}(::Type{T}) = Char

"""
```
*(s::AbstractString, t::AbstractString)
```

Concatenate strings. The `*` operator is an alias to this function.

```jldoctest
julia> "Hello " * "world"
"Hello world"
```
"""
(*)(s1::AbstractString, ss::AbstractString...) = string(s1, ss...)

one{T<:AbstractString}(::Union{T,Type{T}}) = convert(T, "")

length(s::DirectIndexString) = endof(s)

"""
    length(s::AbstractString)

The number of characters in string `s`.

```jldoctest
julia> length("jμΛIα")
5
```
"""
function length(s::AbstractString)
    i = start(s)
    if done(s,i)
        return 0
    end
    n = 1
    while true
        c, j = next(s,i)
        if done(s,j)
            return n
        end
        n += 1
        i = j
    end
end

## string comparison functions ##

function cmp(a::AbstractString, b::AbstractString)
    if a === b
        return 0
    end
    i = start(a)
    j = start(b)
    while !done(a,i)
        if done(b,j)
            return +1
        end
        c, i = next(a,i)
        d, j = next(b,j)
        if c != d
            return c < d ? -1 : +1
        end
    end
    done(b,j) ? 0 : -1
end

==(a::AbstractString, b::AbstractString) = cmp(a,b) == 0
isless(a::AbstractString, b::AbstractString) = cmp(a,b) < 0

# faster comparisons for symbols

cmp(a::Symbol, b::Symbol) = Int(sign(ccall(:strcmp, Int32, (Cstring, Cstring), a, b)))

isless(a::Symbol, b::Symbol) = cmp(a,b) < 0

## Generic validation functions ##

isvalid(s::DirectIndexString, i::Integer) = (start(s) <= i <= endof(s))

"""
    isvalid(str::AbstractString, i::Integer)

Tells whether index `i` is valid for the given string.

```jldoctest
julia> str = "αβγdef";

julia> isvalid(str, 1)
true

julia> str[1]
'α': Unicode U+03b1 (category Ll: Letter, lowercase)

julia> isvalid(str, 2)
false

julia> str[2]
ERROR: UnicodeError: invalid character index
[...]
```
"""
function isvalid(s::AbstractString, i::Integer)
    i < 1 && return false
    done(s,i) && return false
    try
        next(s,i)
        true
    catch
        false
    end
end

## Generic indexing functions ##

prevind(s::DirectIndexString, i::Integer) = Int(i)-1
prevind(s::AbstractArray    , i::Integer) = Int(i)-1
nextind(s::DirectIndexString, i::Integer) = Int(i)+1
nextind(s::AbstractArray    , i::Integer) = Int(i)+1

"""
    prevind(str::AbstractString, i::Integer)

Get the previous valid string index before `i`.
Returns a value less than `1` at the beginning of the string.

```jldoctest
julia> prevind("αβγdef", 3)
1

julia> prevind("αβγdef", 1)
0
```
"""
function prevind(s::AbstractString, i::Integer)
    e = endof(s)
    if i > e
        return e
    end
    j = Int(i)-1
    while j >= 1
        if isvalid(s,j)
            return j
        end
        j -= 1
    end
    return 0 # out of range
end

"""
    nextind(str::AbstractString, i::Integer)

Get the next valid string index after `i`.
Returns a value greater than `endof(str)` at or after the end of the string.

```jldoctest
julia> str = "αβγdef";

julia> nextind(str, 1)
3

julia> endof(str)
9

julia> nextind(str, 9)
10
```
"""
function nextind(s::AbstractString, i::Integer)
    e = endof(s)
    if i < 1
        return 1
    end
    if i > e
        return Int(i)+1
    end
    for j = Int(i)+1:e
        if isvalid(s,j)
            return j
        end
    end
    next(s,e)[2] # out of range
end

checkbounds(s::AbstractString, i::Integer) = start(s) <= i <= endof(s) || throw(BoundsError(s, i))
checkbounds{T<:Integer}(s::AbstractString, r::Range{T}) = isempty(r) || (minimum(r) >= start(s) && maximum(r) <= endof(s)) || throw(BoundsError(s, r))
# The following will end up using a deprecated checkbounds, when T is not Integer
checkbounds{T<:Real}(s::AbstractString, I::AbstractArray{T}) = all(i -> checkbounds(s, i), I)
checkbounds{T<:Integer}(s::AbstractString, I::AbstractArray{T}) = all(i -> checkbounds(s, i), I)

ind2chr(s::DirectIndexString, i::Integer) = begin checkbounds(s,i); i end
chr2ind(s::DirectIndexString, i::Integer) = begin checkbounds(s,i); i end

"""
    ind2chr(s::AbstractString, i::Integer)

Convert a byte index `i` to a character index with
respect to string `s`.

See also [`chr2ind`](@ref).

```jldoctest
julia> str = "αβγdef";

julia> ind2chr(str, 3)
2

julia> chr2ind(str, 2)
3
```
"""
function ind2chr(s::AbstractString, i::Integer)
    s[i] # throws error if invalid
    j = 1
    k = start(s)
    while true
        c, l = next(s,k)
        if i <= k
            return j
        end
        j += 1
        k = l
    end
end

"""
    chr2ind(s::AbstractString, i::Integer)

Convert a character index `i` to a byte index.

See also [`ind2chr`](@ref).

```jldoctest
julia> str = "αβγdef";

julia> chr2ind(str, 2)
3

julia> ind2chr(str, 3)
2
```
"""
function chr2ind(s::AbstractString, i::Integer)
    i < start(s) && throw(BoundsError(s, i))
    j = 1
    k = start(s)
    while true
        c, l = next(s,k)
        if i == j
            return k
        end
        j += 1
        k = l
    end
end

immutable EachStringIndex{T<:AbstractString}
    s::T
end
eachindex(s::AbstractString) = EachStringIndex(s)

length(e::EachStringIndex) = length(e.s)
start(e::EachStringIndex) = start(e.s)
next(e::EachStringIndex, state) = (state, nextind(e.s, state))
done(e::EachStringIndex, state) = done(e.s, state)
eltype(::Type{EachStringIndex}) = Int

## character column width function ##

"""
    strwidth(s::AbstractString)

Gives the number of columns needed to print a string.

```jldoctest
julia> strwidth("March")
5
```
"""
strwidth(s::AbstractString) = (w=0; for c in s; w += charwidth(c); end; w)

"""
    isascii(c::Union{Char,AbstractString}) -> Bool

Tests whether a character belongs to the ASCII character set, or whether this is true for
all elements of a string.
"""
isascii(c::Char) = c < Char(0x80)
isascii(s::AbstractString) = all(isascii, s)

## string promotion rules ##

promote_rule{S<:AbstractString,T<:AbstractString}(::Type{S}, ::Type{T}) = String

"""
    isxdigit(c::Char) -> Bool

Tests whether a character is a valid hexadecimal digit. Note that this does not
include `x` (as in the standard `0x` prefix).

```jldoctest
julia> isxdigit('a')
true

julia> isxdigit('x')
false
```
"""
isxdigit(c::Char) = '0'<=c<='9' || 'a'<=c<='f' || 'A'<=c<='F'

## uppercase, lowercase, and titlecase transformations ##

"""
    uppercase(s::AbstractString)

Returns `s` with all characters converted to uppercase.

```jldoctest
julia> uppercase("Julia")
"JULIA"
```
"""
uppercase(s::AbstractString) = map(uppercase, s)

"""
    lowercase(s::AbstractString)

Returns `s` with all characters converted to lowercase.

```jldoctest
julia> lowercase("STRINGS AND THINGS")
"strings and things"
```
"""
lowercase(s::AbstractString) = map(lowercase, s)

"""
    titlecase(s::AbstractString)

Capitalizes the first character of each word in `s`.

```jldoctest
julia> titlecase("the julia programming language")
"The Julia Programming Language"
```
"""
function titlecase(s::AbstractString)
    startword = true
    b = IOBuffer()
    for c in s
        if isspace(c)
            print(b, c)
            startword = true
        else
            print(b, startword ? titlecase(c) : c)
            startword = false
        end
    end
    return String(take!(b))
end

"""
    ucfirst(s::AbstractString)

Returns `string` with the first character converted to uppercase.

```jldoctest
julia> ucfirst("python")
"Python"
```
"""
function ucfirst(s::AbstractString)
    isempty(s) || isupper(s[1]) ? s : string(uppercase(s[1]),s[nextind(s,1):end])
end

"""
    lcfirst(s::AbstractString)

Returns `string` with the first character converted to lowercase.

```jldoctest
julia> lcfirst("Julia")
"julia"
```
"""
function lcfirst(s::AbstractString)
    isempty(s) || islower(s[1]) ? s : string(lowercase(s[1]),s[nextind(s,1):end])
end

## string map, filter, has ##

function map(f, s::AbstractString)
    out = IOBuffer(StringVector(endof(s)),true,true)
    truncate(out,0)
    for c in s
        c2 = f(c)
        if !isa(c2,Char)
            throw(ArgumentError("map(f,s::AbstractString) requires f to return Char; try map(f,collect(s)) or a comprehension instead"))
        end
        write(out, c2::Char)
    end
    String(take!(out))
end

function filter(f, s::AbstractString)
    out = IOBuffer(StringVector(endof(s)),true,true)
    truncate(out,0)
    for c in s
        if f(c)
            write(out, c)
        end
    end
    String(take!(out))
end


text:
println("these tests will print warnings: ")
@testset "Not installed" begin
    eval(Base, :(is_interactive = true)) # for interactive error handling

    add_format(format"NotInstalled", (), ".not_installed", [:NotInstalled])
    stdin_copy = STDIN
    stderr_copy = STDERR
    rs, wr = redirect_stdin()
    rserr, wrerr = redirect_stderr()
    ref = @async save("test.not_installed")
    println(wr, "y")
    if VERSION < v"0.4.0-dev"
        @test_throws ErrorException wait(ref) #("unknown package NotInstalled")
    else
        @test_throws CompositeException wait(ref) #("unknown package NotInstalled")
    end
    ref = @async save("test.not_installed")
    println(wr, "invalid") #test invalid input
    println(wr, "n") # don't install
    wait(ref)
    @test istaskdone(ref)

    close(rs);close(wr);close(rserr);close(wrerr)
    redirect_stdin(stdin_copy)
    redirect_stderr(stderr_copy)

    eval(Base, :(is_interactive = false)) # for interactive error handling

end

# Missing load/save functions
module BrokenIO
using FileIO
end
add_format(format"BROKEN", (), ".brok", [:BrokenIO])

@testset "Absent implementation" begin
    stderr_copy = STDERR
    rserr, wrerr = redirect_stderr()
    @test_throws FileIO.LoaderError load(Stream(format"BROKEN",STDIN))
    @test_throws FileIO.WriterError save(Stream(format"BROKEN",STDOUT))
    redirect_stderr(stderr_copy)
    close(rserr);close(wrerr)
end

module MultiError1
import FileIO: @format_str, File
load(file::File{format"MultiError"}) = error("1")
end
module MultiError2
import FileIO: @format_str, File, magic
load(file::File{format"MultiError"}) = error("2")
end
@testset "multiple errors" begin
    println("this test will print warnings: ")
    add_format(
        format"MultiError",
        (),
        ".multierr",
        [:MultiError1],
        [:MultiError2]
    )
    ret = @test_throws ErrorException load("test.multierr")
    #@test ret.value.msg == "1" # this is 0.5 only
end


text:
using VDPTag2
using Plots

pyplot()

pomdp = VDPTagPOMDP(mdp=VDPTagMDP(barriers=CardinalBarriers(0.2, 2.8)))

quiver(pomdp)


text:
function filter_all_subpaths!(root_dir, filter_fun::Function,
        container::Vector{AbstractString})
    subs = joinpath.(root_dir, readdir(root_dir));
    for sub in subs
        # filter
        if filter_fun(sub)
            push!(container, sub)
        end
        # recursive call
        if isdir(sub)
            filter_all_subpaths!(sub, filter_fun, container);
        end
    end
    return container;
end

function find_repos(working_dir)
    filter_fun = (path) -> isdir(path) && is_valid_res_repo_tree(path);
    return Vector{ResRepoTree}(ResRepoTree.(filter(filter_fun , joinpath.(working_dir, readdir(working_dir)))));
end

find_by_time_tag(working_dir, tt::Int) =
    find_by_time_tag(working_dir, TimeTag(tt));

function filter_over_all_repos(working_dir, filter_fun)
    repos = find_repos(working_dir);
    res = Vector{AbstractString}();
    for repo in repos
        filter_all_subpaths!(repo.res_dir, filter_fun, res);
    end
    return res;
end

function find_used_time_tags(working_dir)
    return sort!([repo.time_tag for repo in find_repos(working_dir)]);
end

function find_used_max_time_tag(working_dir)
    used_time_tags = find_used_time_tags(working_dir);
    if isempty(used_time_tags)
        return -1;
    else
        return last(used_time_tags)
    end
end

find_versions(working_dir, res_name) =
    sort_by_time_tag(filter((_repo) -> _repo.res_name == res_name,
        find_repos(working_dir)));


text:
x ^ 2

1:ngrid

[end]

x[end] = y[end] + x[end] + z[begin,end]


text:
module SMTPClient

using LibCURL

import Base.convert, Base.show

export send
export SendOptions, SendResponse

def_rto = 0.0

##############################
# Struct definitions
##############################

type SendOptions
    blocking::Bool
    isSSL::Bool
    username::AbstractString
    passwd::AbstractString

    SendOptions(; blocking=true,isSSL=false,  username="", passwd="", ) =
    new(blocking, isSSL, username, passwd)
end

type SendResponse
    body::IO
    code::Int
    total_time::Float64
    SendResponse() = new(IOBuffer(), 0, 0.0)
end

function show(io::IO, o::SendResponse)
    println(io, "Return Code   :", o.code)
    println(io, "Time :", o.total_time)
    println(io, "Response:", takebuf_string(o.body))
end

type ReadData
    typ::Symbol
    src::Any
    str::AbstractString
    offset::Csize_t
    sz::Csize_t

    ReadData() = new(:undefined, false, "", 0, 0)
end

type ConnContext
    curl::Ptr{CURL}
    url::AbstractString
    rd::ReadData
    resp::SendResponse
    options::SendOptions
    close_ostream::Bool
    bytes_recd::Integer

    ConnContext(options::SendOptions) = new(C_NULL, " ", ReadData(), SendResponse(), options, false, 0)
end

##############################
# Callbacks
##############################

function write_cb(buff::Ptr{UInt8}, sz::Csize_t, n::Csize_t, p_ctxt::Ptr{Void})
#    println("@write_cb")
    ctxt = unsafe_pointer_to_objref(p_ctxt)
    nbytes = sz * n
    write(ctxt.resp.body, buff, nbytes)
    ctxt.bytes_recd = ctxt.bytes_recd + nbytes

    nbytes::Csize_t
end

c_write_cb = cfunction(write_cb, Csize_t, (Ptr{UInt8}, Csize_t, Csize_t, Ptr{Void}))

function curl_read_cb(out::Ptr{Void}, s::Csize_t, n::Csize_t, p_ctxt::Ptr{Void})
#    println("@curl_read_cb")

    ctxt = unsafe_pointer_to_objref(p_ctxt)
    bavail::Csize_t = s * n
    breq::Csize_t = ctxt.rd.sz - ctxt.rd.offset
    b2copy = bavail > breq ? breq : bavail

    if (ctxt.rd.typ == :buffer)
        ccall(:memcpy, Ptr{Void}, (Ptr{Void}, Ptr{Void}, UInt),
                out, convert(Ptr{UInt8}, ctxt.rd.str) + ctxt.rd.offset, b2copy)
    elseif (ctxt.rd.typ == :io)
        b_read = read(ctxt.rd.src, UInt8, b2copy)
        ccall(:memcpy, Ptr{Void}, (Ptr{Void}, Ptr{Void}, UInt), out, b_read, b2copy)
    end
    ctxt.rd.offset = ctxt.rd.offset + b2copy

    r = convert(Csize_t, b2copy)
    r::Csize_t
end

c_curl_read_cb = cfunction(curl_read_cb, Csize_t, (Ptr{Void}, Csize_t, Csize_t, Ptr{Void}))

function curl_multi_timer_cb(curlm::Ptr{Void}, timeout_ms::Clong, p_muctxt::Ptr{Void})
    muctxt = unsafe_pointer_to_objref(p_muctxt)
    muctxt.timeout = timeout_ms / 1000.0

#    println("Requested timeout value : " * string(muctxt.timeout))

    ret = convert(Cint, 0)
    ret::Cint
end

c_curl_multi_timer_cb = cfunction(curl_multi_timer_cb, Cint, (Ptr{Void}, Clong, Ptr{Void}))

##############################
# Utility functions
##############################

macro ce_curl(f, args...)
    quote
        cc = CURLE_OK
        cc = $(esc(f))(ctxt.curl, $(args...))

        if(cc != CURLE_OK)
            error(string($f) * "() failed: " * bytestring(curl_easy_strerror(cc)))
        end
    end
end

macro ce_curlm(f, args...)
    quote
        cc = CURLM_OK
        cc = $(esc(f))(curlm, $(args...))

        if(cc != CURLM_OK)
            error(string($f) * "() failed: " * bytestring(curl_multi_strerror(cc)))
        end
    end
end

null_cb(curl) = return nothing

function set_opt_blocking(options::SendOptions)
        o2 = deepcopy(options)
        o2.blocking = true
        return o2
end

function setup_easy_handle(url, options::SendOptions)
    ctxt = ConnContext(options)

    curl = curl_easy_init()
    if (curl == C_NULL) throw("curl_easy_init() failed") end

    ctxt.curl = curl

    ctxt.url = url

    p_ctxt = pointer_from_objref(ctxt)

    @ce_curl curl_easy_setopt CURLOPT_URL url
    @ce_curl curl_easy_setopt CURLOPT_WRITEFUNCTION c_write_cb
    @ce_curl curl_easy_setopt CURLOPT_UPLOAD 1
    @ce_curl curl_easy_setopt CURLOPT_WRITEDATA p_ctxt

    if options.isSSL
        @ce_curl curl_easy_setopt CURLOPT_USE_SSL CURLUSESSL_ALL
    end

    if length(options.username) > 0
        @ce_curl curl_easy_setopt  CURLOPT_USERNAME options.username
        @ce_curl curl_easy_setopt  CURLOPT_PASSWORD options.passwd
    end

    ctxt
end

function cleanup_easy_context(ctxt::Union{ConnContext,Bool})
    if isa(ctxt, ConnContext)

        if (ctxt.curl != C_NULL)
            curl_easy_cleanup(ctxt.curl)
        end

        if ctxt.close_ostream
            close(ctxt.resp.body)
            ctxt.resp.body = nothing
            ctxt.close_ostream = false
        end
    end
end

function process_response(ctxt)
    http_code = Array(Int,1)
    @ce_curl curl_easy_getinfo CURLINFO_RESPONSE_CODE http_code

    total_time = Array(Float64,1)
    @ce_curl curl_easy_getinfo CURLINFO_TOTAL_TIME total_time

    ctxt.resp.code = http_code[1]
    ctxt.resp.total_time = total_time[1]
end

##############################
# Library initializations
##############################

init() = curl_global_init(CURL_GLOBAL_ALL)
cleanup() = curl_global_cleanup()

function send(url::AbstractString, to::Vector, from::AbstractString, body::IO, options::SendOptions=SendOptions())
    if (options.blocking)
        rd::ReadData = ReadData()

            rd.typ = :io
            rd.src = body
            seekend(body)
            rd.sz = position(body)
            seekstart(body)

        return _do_send(url,to, from, options, rd)
    else
        return remotecall(myid(), send, url, to, from, body, set_opt_blocking(options))
    end

end

function _do_send(url::AbstractString, to::Vector, from::AbstractString, options::SendOptions, rd::ReadData)
    ctxt = false
    slist::Ptr{Void} = C_NULL
    try
        ctxt = setup_easy_handle(url, options)
        ctxt.rd = rd
        # rd.typ is always IO for smtp

        p_ctxt = pointer_from_objref(ctxt)
        @ce_curl curl_easy_setopt CURLOPT_READDATA p_ctxt

        @ce_curl curl_easy_setopt CURLOPT_READFUNCTION c_curl_read_cb

        for tos in to
            slist = curl_slist_append(slist, tos)
        end

        @ce_curl curl_easy_setopt CURLOPT_MAIL_RCPT slist

        @ce_curl curl_easy_setopt CURLOPT_MAIL_FROM from

        # return exec_as_multi(ctxt)

        @ce_curl curl_easy_perform
        process_response(ctxt)
        return ctxt.resp
    finally
        if (slist != C_NULL)
            curl_slist_free_all(slist)
        end
        cleanup_easy_context(ctxt)
    end
end

function exec_as_multi(ctxt)
    curl = ctxt.curl
    curlm = curl_multi_init()

    if (curlm == C_NULL) error("Unable to initialize curl_multi_init()") end

    try
        if isa(ctxt.options.callback, Function) ctxt.options.callback(curl) end

        @ce_curlm curl_multi_add_handle curl

        n_active = Array(Int,1)
        n_active[1] = 1

        no_to = 30 * 24 * 3600.0
        request_timeout = 0.001 + (ctxt.options.request_timeout == 0.0 ? no_to : ctxt.options.request_timeout)

        started_at = time()
        time_left = request_timeout

    # poll_fd is unreliable when multiple parallel fds are active, hence using curl_multi_perform

# START curl_multi_perform  mode

        cmc = curl_multi_perform(curlm, n_active);
        while (n_active[1] > 0) &&  (time_left > 0)
            nb1 = ctxt.bytes_recd
            cmc = curl_multi_perform(curlm, n_active);
            if(cmc != CURLM_OK) error("curl_multi_perform() failed: " * bytestring(curl_multi_strerror(cmc))) end

            nb2 = ctxt.bytes_recd

            if (nb2 > nb1)
                yield() # Just yield to other tasks
            else
                sleep(0.005) # Just to prevent unnecessary CPU spinning
            end

            time_left = request_timeout - (time() - started_at)
        end

# END OF curl_multi_perform

        if (n_active[1] == 0)
            msgs_in_queue = Array(Cint,1)
            p_msg::Ptr{CURLMsg2} = curl_multi_info_read(curlm, msgs_in_queue)

            while (p_msg != C_NULL)
#                println("Messages left in Q : " * string(msgs_in_queue[1]))
                msg = unsafe_load(p_msg)

                if (msg.msg == CURLMSG_DONE)
                    ec = convert(Int, msg.data)
                    if (ec != CURLE_OK)
#                        println("Result of transfer: " * string(msg.data))
                        throw("Error executing request : " * bytestring(curl_easy_strerror(ec)))
                    else
                        process_response(ctxt)
                    end
                end

                p_msg = curl_multi_info_read(curlm, msgs_in_queue)
            end
        else
            error("request timed out")
        end

    finally
        curl_multi_remove_handle(curlm, curl)
        curl_multi_cleanup(curlm)
    end

    ctxt.resp
end

end # module


text:
include("..\\libs\\FunctionsRobot.jl")

function moveToWall!(r::Robot)
    c = 0
    while !isborder(r, West)
        putmarker!(r)
        move!(r, West)
        c+=1
    end
    putmarker!(r)
    return c
end

function moveAndPuting!(r)
    c = moveToWall!(r)
    while c > 0
        if isborder(r, Nord)
            break
        end
        move!(r, Nord)
        for i in 1:c-1
            move!(r, Ost)
        end
        c = moveToWall!(r)
    end
end

function lestniza!(r::Robot)
    sides = moveAndReturnDirections!(r)
    moveAndPuting!(r)
    moveToStartplace!(r)
    moveToBeginplace!(r, sides)
end