text:
# Polyharmonic functions used below.
# function polyharmonic_phi1(r) return r; end
# function polyharmonic_phi2(r) return r>1e-16 ? r^2 * log(r) : r; end
# function polyharmonic_phi3(r) return r^3; end
# function polyharmonic_phi4(r) return r>1e-16 ? r^4 * log(r) : r^3; end

# Uses polyharmonic interpolation to interpolate at x using values u at centers c.
# sizes are:
# x is d x nx
# c is d x nc
# u is nc
# (optional) ord determines the type of basis function. default is 3. 1,2,3,4 are available
# Note: This now redirects to a dimension specific function
function polyharmonic_interp(x, c, u; ord=3, limiter=0, neighbor_vals=[])
    dim = size(x,1);
    if dim==1
        return polyharmonic_interp_1d(x, c, u, ord=ord, limiter=limiter, neighbor_vals=neighbor_vals);
    elseif dim==2
        return polyharmonic_interp_2d(x, c, u, ord=ord, limiter=limiter, neighbor_vals=neighbor_vals);
    else
        return polyharmonic_interp_3d(x, c, u, ord=ord, limiter=limiter, neighbor_vals=neighbor_vals);
    end
end

# Uses polyharmonic interpolation to interpolate at x using values u at centers c.
# See above for input details.
function polyharmonic_interp_1d(x, c, u; ord=3, limiter=0, neighbor_vals=[])
    dim = 1;
    nx = size(x,2); # number of places to evaluate

    # If no centers are provided, this is meaningless, but return zero.
    if length(c)==0
        return zeros(nx);
    end

    nc = size(c,2); # number of centers

    # Constant if only one center is given
    if nc == 1
        if nx==1
            return u;
        else
            return u[1] .* ones(nx);
        end

    elseif nc == 2 # Linear if two centers are given
        du = u[2] - u[1];
        dx1 = c[2] - c[1];
        dx2 = x .- c[1];
        if abs(dx1) > 1e-20
            sol = u[1] .+ (du/dx1 .* dx2);
        else # overlapping points
            sol = u[1] .* ones(nx);
        end

    else # Use polyharmonics if more than two centers are given
        # build A matrix
        A = zeros(nc + dim+1, nc + dim+1);
        # Handle each order separately for efficiency
        if ord == 1
            for i=1:nc
                for j=(i+1):nc
                    A[i,j] = norm(c[:,j] - c[:,i]); # phi(norm(c[:,j] - c[:,i]))
                    A[j,i] = A[i,j];
                end
            end
        elseif ord == 2
            for i=1:nc
                for j=(i+1):nc
                    # A[i,j] = norm(c[:,j] - c[:,i]); # phi(norm(c[:,j] - c[:,i]))
                    tmp = sum((c[:,j] - c[:,i]).^2);
                    if tmp > 1e-30
                        A[i,j] = tmp * log(tmp) * 0.5;
                    else
                        A[i,j] = sqrt(tmp);
                    end
                    A[j,i] = A[i,j];
                end
            end
        elseif ord == 3
            for i=1:nc
                for j=(i+1):nc
                    A[i,j] = (norm(c[:,j] - c[:,i]))^3; # phi(norm(c[:,j] - c[:,i]))
                    A[j,i] = A[i,j];
                end
            end
        else # max order is 4
            for i=1:nc
                for j=(i+1):nc
                    # A[i,j] = norm(c[:,j] - c[:,i]); # phi(norm(c[:,j] - c[:,i]))
                    tmp = sum((c[:,j] - c[:,i]).^2);
                    if tmp > 1e-30
                        A[i,j] = tmp*tmp * log(tmp) * 0.5;
                    else
                        A[i,j] = tmp*sqrt(tmp);
                    end
                    A[j,i] = A[i,j];
                end
            end
        end

        # build B matrix (within A matrix)
        A[1:nc, nc+1] .= 1;
        A[nc+1, 1:nc] .= 1;
        A[1:nc, (nc+2):(nc+dim+1)] = c';
        A[(nc+2):(nc+dim+1), 1:nc] = c;

        # build RHS vector
        b = vcat(u, zeros(dim+1));

        wv = A\b; # coefficients for the polyharmonic

        sol = zeros(nx);
        for i=1:nx
            # evaluate at x[i]
            sol[i] = wv[nc+1];
            for k=1:dim
                sol[i] += wv[nc+k+1] * x[k,i];
            end
            if ord == 1
                for j=1:nc
                    sol[i] += wv[j]*norm(c[:,j] - x[:,i]); # wv[j]*phi(norm(c[:,j] - x[:,i]));
                end
            elseif ord == 2
                for j=1:nc
                    # sol[i] += wv[j]*phi(norm(c[:,j] - x[:,i]));
                    tmp = sum((c[:,j] - x[:,i]).^2);
                    if tmp > 1e-30
                        sol[i] += wv[j] * tmp * log(tmp) * 0.5;
                    else
                        sol[i] += wv[j] * sqrt(tmp);
                    end
                end
            elseif ord == 3
                for j=1:nc
                    sol[i] += wv[j]*norm(c[:,j] - x[:,i])^3;
                end
            else # 4+
                for j=1:nc
                    # sol[i] += wv[j]*phi(norm(c[:,j] - x[:,i]));
                    tmp = sum((c[:,j] - x[:,i]).^2);
                    if tmp > 1e-30
                        sol[i] += wv[j] * tmp * tmp * log(tmp) * 0.5;
                    else
                        sol[i] += wv[j] * tmp*sqrt(tmp);
                    end
                end
            end
        end # building sol
    end
    # At this point sol contains the extrapolated value

    # If a limiting is desired, limit the extrapolated values to the range of the neighboring cell values
    if limiter > 0 && length(neighbor_vals) == 2
        maxv = max(neighbor_vals[1], neighbor_vals[2]);
        minv = min(neighbor_vals[1], neighbor_vals[2]);

        for i=1:nx
            #sol[i] = max(sol[i], minv);
            #sol[i] = min(sol[i], maxv);
            if sol[i] < minv
                sol[i] = limiter * minv + (1-limiter) * sol[i];
            elseif sol[i] > maxv
                sol[i] = limiter * maxv + (1-limiter) * sol[i];
            end
        end
    end

    return sol;
end

# Uses polyharmonic interpolation to interpolate at x using values u at centers c.
# See above polyharmonic_interp() for input details
function polyharmonic_interp_2d(x, c, u; ord=3, limiter=0, neighbor_vals=[])
    dim = 2;
    nx = size(x,2); # number of places to evaluate

    # If no centers are provided, this is meaningless, but return zero.
    if length(c)==0
        return zeros(nx);
    end

    nc = size(c,2); # number of centers

    # Constant if only one center is given
    if nc == 1
        if nx==1
            return u;
        else
            return u[1] .* ones(nx);
        end

    elseif nc == 2 # NOT A UNIQUE INTERPOLATION, so assume the gradient is parallel to the vector between the two points
        sol = linear_interp(x, c, u, 2);

    else # Use polyharmonics if more than two centers are given and at least three create a nondegenerate triangle
        # First check for nondegeneracy
        degenerate = true;
        v1 = c[:,2] - c[:,1];
        for i=3:nc
            v2 = c[:,i] - c[:,1];
            if (v1[1]*v2[2] - v1[2]*v2[1]) > 1e-20 # rotate v2 by pi/2 and dot
                degenerate = false;
            end
        end

        if degenerate
            sol = linear_interp(x, c[:,1:2], u[1:2], 2);

        else
            # build A matrix
            A = zeros(nc + dim+1, nc + dim+1);
            # Handle each order separately for efficiency
            if ord == 1
                for i=1:nc
                    for j=(i+1):nc
                        A[i,j] = norm(c[:,j] - c[:,i]); # phi(norm(c[:,j] - c[:,i]))
                        A[j,i] = A[i,j];
                    end
                end
            elseif ord == 2
                for i=1:nc
                    for j=(i+1):nc
                        # A[i,j] = norm(c[:,j] - c[:,i]); # phi(norm(c[:,j] - c[:,i]))
                        tmp = sum((c[:,j] - c[:,i]).^2);
                        if tmp > 1e-30
                            A[i,j] = tmp * log(tmp) * 0.5;
                        else
                            A[i,j] = sqrt(tmp);
                        end
                        A[j,i] = A[i,j];
                    end
                end
            elseif ord == 3
                for i=1:nc
                    for j=(i+1):nc
                        A[i,j] = (norm(c[:,j] - c[:,i]))^3; # phi(norm(c[:,j] - c[:,i]))
                        A[j,i] = A[i,j];
                    end
                end
            else # max order is 4
                for i=1:nc
                    for j=(i+1):nc
                        # A[i,j] = norm(c[:,j] - c[:,i]); # phi(norm(c[:,j] - c[:,i]))
                        tmp = sum((c[:,j] - c[:,i]).^2);
                        if tmp > 1e-30
                            A[i,j] = tmp*tmp * log(tmp) * 0.5;
                        else
                            A[i,j] = tmp*sqrt(tmp);
                        end
                        A[j,i] = A[i,j];
                    end
                end
            end

            # build B matrix (within A matrix)
            A[1:nc, nc+1] .= 1;
            A[nc+1, 1:nc] .= 1;
            A[1:nc, (nc+2):(nc+dim+1)] = c';
            A[(nc+2):(nc+dim+1), 1:nc] = c;

            # build RHS vector
            b = vcat(u, zeros(dim+1));

            wv = A\b; # coefficients for the polyharmonic

            sol = zeros(nx);
            for i=1:nx
                # evaluate at x[i]
                sol[i] = wv[nc+1];
                for k=1:dim
                    sol[i] += wv[nc+k+1] * x[k,i];
                end
                if ord == 1
                    for j=1:nc
                        sol[i] += wv[j]*norm(c[:,j] - x[:,i]); # wv[j]*phi(norm(c[:,j] - x[:,i]));
                    end
                elseif ord == 2
                    for j=1:nc
                        # sol[i] += wv[j]*phi(norm(c[:,j] - x[:,i]));
                        tmp = sum((c[:,j] - x[:,i]).^2);
                        if tmp > 1e-30
                            sol[i] += wv[j] * tmp * log(tmp) * 0.5;
                        else
                            sol[i] += wv[j] * sqrt(tmp);
                        end
                    end
                elseif ord == 3
                    for j=1:nc
                        sol[i] += wv[j]*norm(c[:,j] - x[:,i])^3;
                    end
                else # 4+
                    for j=1:nc
                        # sol[i] += wv[j]*phi(norm(c[:,j] - x[:,i]));
                        tmp = sum((c[:,j] - x[:,i]).^2);
                        if tmp > 1e-30
                            sol[i] += wv[j] * tmp * tmp * log(tmp) * 0.5;
                        else
                            sol[i] += wv[j] * tmp*sqrt(tmp);
                        end
                    end
                end
            end # building sol
        end

    end
    # At this point sol contains the extrapolated value

    # If a limiting is desired, limit the extrapolated values to the range of the neighboring cell values
    if limiter > 0 && length(neighbor_vals) == 2
        maxv = max(neighbor_vals[1], neighbor_vals[2]);
        minv = min(neighbor_vals[1], neighbor_vals[2]);

        for i=1:nx
            #sol[i] = max(sol[i], minv);
            #sol[i] = min(sol[i], maxv);
            if sol[i] < minv
                sol[i] = limiter * minv + (1-limiter) * sol[i];
            elseif sol[i] > maxv
                sol[i] = limiter * maxv + (1-limiter) * sol[i];
            end
        end
    end

    return sol;
end

# Uses polyharmonic interpolation to interpolate at x using values u at centers c.
# See above polyharmonic_interp() for input details
function polyharmonic_interp_3d(x, c, u; ord=3, limiter=0, neighbor_vals=[])
    dim = 2;
    nx = size(x,2); # number of places to evaluate

    # If no centers are provided, this is meaningless, but return zero.
    if length(c)==0
        return zeros(nx);
    end

    nc = size(c,2); # number of centers

    # Constant if only one center is given
    if nc == 1
        if nx==1
            return u;
        else
            return u[1] .* ones(nx);
        end

    else
        #TODO
        sol = u[1] .* ones(nx);
    end

    return sol;
end

# Linear interpolation used when polyharmonic won't work.
# It is assumed that there are two centers
function linear_interp(x, c, u, dim)
    nx = size(x,2); # number of places to evaluate
    if dim==1
        du = u[2] - u[1];
        dx1 = c[2] - c[1];
        dx2 = x .- c[1];
        if abs(dx1) > 1e-20
            sol = u[1] .+ (du/dx1 .* dx2);
        else # overlapping points
            sol = u[1] .* ones(nx);
        end

    elseif dim==2 # NOT A UNIQUE INTERPOLATION, so assume the gradient is parallel to the vector between the two points
        du = u[2] - u[1];
        dx1 = [c[1,2] - c[1,1]; c[2,2] - c[2,1]];
        if nx == 1
            dx2 = [x[1,1] - c[1,1]; x[2,1] - c[2,1]];
        else
            dx2 = similar(x);
            for xi=1:nx
                dx2[:,xi] = [x[1,xi] - c[1,1]; x[2,xi] - c[2,1]];
            end
        end
        if abs(dx1[1])+abs(dx1[2]) > 1e-20
            # u + grad(u) . dx2
            if nx == 1
                sol = [u[1] + (du/(dx1[1]^2 + dx1[2]^2) * (dx1[1]*dx2[1] + dx1[2]*dx2[2]))];
            else
                sol = zeros(nx);
                for xi=1:nx
                    sol[xi] = u[1] .+ (du/(dx1[1]^2 + dx1[2]^2) .* (dx1[1]*dx2[1,xi] + dx1[2]*dx2[2,xi]));
                end
            end
        else # overlapping points
            sol = u[1] .* ones(nx);
        end

    else # 3D
        #TODO
        sol = u[1] .* ones(nx);
    end

    return sol;
end

max_stars_repo_name: paralab/Finch

max_issues_repo_name: paralab/Finch

max_forks_repo_name: paralab/Finch


text:
@testset "CacheImage" begin
    @test (CacheImage <: Augmentor.AffineOperation) == false
    @test (CacheImage <: Augmentor.ImageOperation) == true
    @test typeof(@inferred(CacheImage())) <: CacheImage
    @test str_show(CacheImage()) == "Augmentor.CacheImage()"
    @test str_showconst(CacheImage()) == "CacheImage()"
    @test str_showcompact(CacheImage()) == "Cache into temporary buffer"

    @test @inferred(Augmentor.applyeager(CacheImage(),square)) === square
    @test @inferred(Augmentor.applyeager(CacheImage(),(square,))) === (square,)
    @test @inferred(Augmentor.applyeager(CacheImage(),(square,square2))) === (square,square2)

    # SubArray
    v = view(square, :, :)
    img = @inferred Augmentor.applyeager(CacheImage(), v)
    @test typeof(img) <: Array
    @test eltype(img) == eltype(v)
    @test img !== square
    @test img == square
    # Identidy ranges
    v = view(square, IdentityRange(1:3), IdentityRange(1:3))
    img = @inferred Augmentor.applyeager(CacheImage(), v)
    @test typeof(img) <: OffsetArray
    @test eltype(img) == eltype(v)
    @test img == OffsetArray(square, 0, 0)
    # Affine
    v = Augmentor.prepareaffine(square)
    img = @inferred Augmentor.applyeager(CacheImage(), v)
    @test typeof(img) <: OffsetArray
    @test eltype(img) == eltype(v)
    @test img == OffsetArray(square, 0, 0)
    # Array and SubArray
    v = view(square, :, :)
    tmp,img = @inferred Augmentor.applyeager(CacheImage(), (square,v))
    @test typeof(img) <: Array
    @test eltype(img) == eltype(v)
    @test tmp === square
    @test img !== square
    @test img == square
    # OffsetArray
    o = OffsetArray(square, (-1,2))
    @test @inferred(Augmentor.applyeager(CacheImage(),o)) === o
    @test @inferred(Augmentor.applyeager(CacheImage(),(o,square))) === (o,square)

    @test Augmentor.supports_eager(CacheImage) === true
    @test Augmentor.supports_lazy(CacheImage) === false
    @test Augmentor.supports_view(CacheImage) === false
    @test Augmentor.supports_stepview(CacheImage) === false
    @test Augmentor.supports_permute(CacheImage) === false
    @test Augmentor.supports_affine(CacheImage) === false
    @test Augmentor.supports_affineview(CacheImage) === false

    @test_throws MethodError Augmentor.applylazy(CacheImage(), v)
    @test_throws MethodError Augmentor.applylazy(CacheImage(), (v,v))
    @test_throws MethodError Augmentor.applyview(CacheImage(), v)
    @test_throws MethodError Augmentor.applystepview(CacheImage(), v)
    @test_throws MethodError Augmentor.applypermute(CacheImage(), v)
    @test_throws MethodError Augmentor.applyaffine(CacheImage(), v)
    @test_throws MethodError Augmentor.applyaffineview(CacheImage(), v)
end

# --------------------------------------------------------------------

@testset "CacheImageInto" begin
    @test_throws UndefVarError CacheImageInto
    @test (Augmentor.CacheImageInto <: Augmentor.AffineOperation) == false
    @test (Augmentor.CacheImageInto <: Augmentor.ImageOperation) == true
    @test_throws MethodError Augmentor.CacheImageInto()

    @testset "single image" begin
        buf = copy(rect)
        @test typeof(@inferred(CacheImage(buf))) <: Augmentor.CacheImageInto
        op = @inferred CacheImage(buf)
        @test Augmentor.CacheImageInto(buf) === op
        @test str_show(op) == "Augmentor.CacheImageInto(::Array{Gray{N0f8},2})"
        @test str_showconst(op) == "CacheImage(Array{Gray{N0f8}}(2, 3))"
        op2 = @inferred CacheImage(Array{Gray{N0f8}}(2, 3))
        @test typeof(op) == typeof(op2)
        @test typeof(op.buffer) == typeof(op2.buffer)
        @test size(op.buffer) == size(op2.buffer)
        @test str_showcompact(op) == "Cache into preallocated 2×3 Array{Gray{N0f8},2}"

        v = Augmentor.applylazy(Resize(2,3), camera)
        res = @inferred Augmentor.applyeager(op, v)
        @test res == v
        @test typeof(res) <: OffsetArray
        @test parent(res) === op.buffer

        res = @inferred Augmentor.applyeager(op, rect)
        @test res == rect
        @test res === op.buffer

        res = @inferred Augmentor.applylazy(op, v)
        @test res == v
        @test typeof(res) <: OffsetArray
        @test parent(res) === op.buffer

        res = @inferred Augmentor.applylazy(op, rect)
        @test res == rect
        @test res === op.buffer

        @test_throws BoundsError Augmentor.applyeager(op, camera)
        @test_throws MethodError Augmentor.applyview(CacheImage(buf), v)
        @test_throws MethodError Augmentor.applystepview(CacheImage(buf), v)
        @test_throws MethodError Augmentor.applypermute(CacheImage(buf), v)
        @test_throws MethodError Augmentor.applyaffine(CacheImage(buf), v)
        @test_throws MethodError Augmentor.applyaffineview(CacheImage(buf), v)
    end

    @testset "multiple images" begin
        buf1 = copy(square)
        buf2 = copy(rgb_rect)
        @test typeof(@inferred(CacheImage(buf1,buf2))) <: Augmentor.CacheImageInto
        op = @inferred CacheImage(buf1,buf2)
        @test op === @inferred CacheImage((buf1,buf2))
        @test Augmentor.CacheImageInto((buf1,buf2)) === op
        @test str_show(op) == "Augmentor.CacheImageInto((::Array{Gray{N0f8},2}, ::Array{RGB{N0f8},2}))"
        @test str_showconst(op) == "CacheImage(Array{Gray{N0f8}}(3, 3), Array{RGB{N0f8}}(2, 3))"
        op2 = @inferred CacheImage(Array{Gray{N0f8}}(3, 3), Array{RGB{N0f8}}(2, 3))
        @test typeof(op) == typeof(op2)
        @test typeof(op.buffer) == typeof(op2.buffer)
        @test size.(op.buffer) === size.(op2.buffer)
        @test str_showcompact(op) == "Cache into preallocated (3×3 Array{Gray{N0f8},2}, 2×3 Array{RGB{N0f8},2})"

        @test buf1 == square
        @test buf2 == rgb_rect
        v1 = Augmentor.applylazy(Resize(3,3), camera)
        v2 = Augmentor.applylazy(Resize(2,3), RGB.(camera))
        res = @inferred Augmentor.applyeager(op, (v1,v2))
        @test buf1 != square
        @test buf2 != rgb_rect
        @test res == (v1, v2)
        @test typeof(res) <: NTuple{2,OffsetArray}
        @test parent.(res) === (op.buffer[1], op.buffer[2])

        @test_throws BoundsError Augmentor.applyeager(op, (camera,buf1))
        @test_throws MethodError Augmentor.applylazy(op, v1)
        @test_throws BoundsError Augmentor.applylazy(op, (buf2,buf1))
        @test_throws BoundsError Augmentor.applylazy(op, (buf1,))
        # ?
        @test_throws DimensionMismatch Augmentor.applylazy(op, (v1,v1))
    end

    @test Augmentor.supports_eager(Augmentor.CacheImageInto) === true
    @test Augmentor.supports_lazy(Augmentor.CacheImageInto) === true
    @test Augmentor.supports_view(Augmentor.CacheImageInto) === false
    @test Augmentor.supports_stepview(Augmentor.CacheImageInto) === false
    @test Augmentor.supports_permute(Augmentor.CacheImageInto) === false
    @test Augmentor.supports_affine(Augmentor.CacheImageInto) === false
    @test Augmentor.supports_affineview(Augmentor.CacheImageInto) === false
end

ext: jl


text:
"""
    Resize <: Augmentor.ImageOperation

Description
--------------

Rescales the image to a fixed pre-specified pixel size.

This operation does not take any measures to preserve aspect
ratio of the source image. Instead, the original image will
simply be resized to the given dimensions. This is useful when
one needs a set of images to all be of the exact same size.

Usage
--------------

    Resize(; height=64, width=64)

    Resize(size)

    Resize(size...)

Arguments
--------------

- **`size`** : `NTuple` or `Vararg` of `Int` that denote the
    output size in pixel for each dimension.

See also
--------------

[`CropSize`](@ref), [`augment`](@ref)

Examples
--------------

```julia
using Augmentor
img = testpattern()

augment(img, Resize(30, 40))
```
"""
struct Resize{N} <: ImageOperation
    size::NTuple{N,Int}

    function Resize{N}(size::NTuple{N,Int}) where N
        all(s->s>0, size) || throw(ArgumentError("Specified sizes must be strictly greater than 0. Actual: $size"))
        new{N}(size)
    end
end
Resize(::Tuple{}) = throw(MethodError(Resize, ((),)))
Resize(size::NTuple{N,Int}) where {N} = Resize{N}(size)
Resize(size::Int...) = Resize(size)
Resize(; width=64, height=64) = Resize((height,width))

@inline supports_affineview(::Type{<:Resize}) = true

function toaffinemap(op::Resize{2}, img::AbstractMatrix)
    # emulate behaviour of ImageTransformations.imresize!
    Rin  = CartesianRange(indices(img))
    sf = map(/, op.size, (last(Rin)-first(Rin)+1).I)
    offset = map((io,ir,s)->io - 0.5 - s*(ir-0.5), first(Rin).I, (1, 1), map(inv,sf))
    ttrans = AffineMap(@SMatrix([1. 0.; 0. 1.]), SVector(offset))
    tscale = recenter(@SMatrix([sf[1] 0.; 0. sf[2]]), @SVector([1., 1.]))
    tscale ∘ ttrans
end

function applyeager(op::Resize, img::AbstractArray, param)
    plain_array(imresize(img, op.size))
end

function applylazy(op::Resize, img::AbstractArray, param)
    applyaffineview(op, prepareaffine(img), param)
end

function padrange(range::AbstractUnitRange, pad)
    first(range)-pad:last(range)+pad
end

function applyaffineview(op::Resize{N}, img::AbstractArray{T,N}, param) where {T,N}
    Rin, Rout = CartesianRange(indices(img)), CartesianRange(op.size)
    sf = map(/, (last(Rout)-first(Rout)+1).I, (last(Rin)-first(Rin)+1).I)
    # We have to extrapolate if the image is upscaled,
    # otherwise the original border will only cause a single pixel
    tinv = toaffinemap(op, img, param)
    inds = ImageTransformations.autorange(img, tinv)
    pad_inds = map((s,r)-> s>=1 ? padrange(r,ceil(Int,s/2)) : r, sf, inds)
    wv = invwarpedview(img, tinv, pad_inds)
    # expanding will cause an additional pixel that has to be skipped
    indirect_view(wv, map(s->2:s+1, op.size))
end

function showconstruction(io::IO, op::Resize)
    print(io, typeof(op).name.name, '(', join(map(string, op.size),", "), ')')
end

function Base.show(io::IO, op::Resize{N}) where {N}
    if get(io, :compact, false)
        if N == 2
            print(io, "Resize to $(op.size[1])×$(op.size[2])")
        else
            print(io, "Resize to $(op.size)")
        end
    else
        print(io, typeof(op), "($(op.size))")
    end
end


text:
Pkg.add("PyPlot")
Pkg.build("PyPlot")


text:
# generalised OU
include("common.jl")
import SDELab2
################### problem parameters
T=1 # time step, final time
x0=1. # initial data
alp=-1/2 # drift parameter
bet= 4/5 # drift parameter
sig=sqrt(1/2) # diffusion param
quad0=SDELab2.delta_measure(x0)
###########################################
function drift(u,alp,bet)
  x=u[1]; mf=u[2]
  return alp*x+bet*mf
end
##########################################
function diff(u,sig)
  return sig
end
#######################################
function mf_r(x)
  return [x]
end
##########################################
fcn=SDELab2.set_fcn_mf(u->drift(u,alp,bet),
          u->diff(u,sig),
          x->mf_r(x))
# test function
phi=x->[x;x.^2]
###################################### exact
exact_mean=x0*exp((alp+bet));
exact_var=(sig^2)/(2*alp)*exp(2*alp)-sig^2/(2*alp);
exact_mom2=exact_var+exact_mean^2;
exact=[exact_mean,exact_mom2];
println("ode soln ",exact)
##################################### method parameters
dt=1e-3
solver="MFEM0"
params=SDELab2.set_opt(dt,solver)
params["cutoff"]=10 # radius R for Alg 3.1
params["const"]=1e3 # for no_Gauss_points
params["length_scale"]=1 # for no_Gauss_points
params["skiphalf"]=true
params["periodic"]=false
params["no_partition"]=1
dtpts=[0.2,0.1,0.05,0.025,0.01,0.005]#,0.0025,0.001]
errors=zeros(length(dtpts),2)
errors_no=deepcopy(errors)
errors_nd=deepcopy(errors)
cpu_times=zeros(length(dtpts))
cpu_times_no=deepcopy(cpu_times)
cpu_times_nd=deepcopy(cpu_times)

params["Solver"]="MFEM0"
dtpts_no=dtpts/4
for i=1:(length(dtpts))
    params["MaxStepSize"]=dtpts_no[i]
    qoi,cpu_time,stats,q=get_soln(fcn,quad0,params,phi,false)
    errors_no[i,:],cpu_times_no[i]=helper(qoi, exact, cpu_time, "normal")
end
########################
params["Solver"]="MFEM0"
for i=1:(length(dtpts))
    params["MaxStepSize"]=dtpts[i]
    qoi,cpu_time,stats,q=get_soln(fcn,quad0,params,phi,true)
    errors[i,:],cpu_times[i]=helper(qoi, exact, cpu_time, "normal")
end
################
params["Solver"]="MFEM2"
for i=1:(length(dtpts))
    params["MaxStepSize"]=dtpts[i]
    qoi,cpu_time,stats,q=get_soln(fcn,quad0,params,phi,false)
    errors_nd[i,:],cpu_times_nd[i]=helper(qoi, exact, cpu_time, "normal")
end
#########################################

figure(1);
# first functional
i=1
time_step_err(dtpts,errors[:,i],
  dtpts_no,errors_no[:,i],
  dtpts,errors_nd[:,i],[1e0,1e0],utf8("mom1.pdf"))
err_cpu(errors[:,i],cpu_times[:],
  errors_no[:,i],cpu_times_no[:],
  errors_nd[:,i],cpu_times_nd[:],[1e0,1e0], utf8("mom1_.pdf"))
# second functional
i=2
time_step_err(dtpts,errors[:,i],
      dtpts_no,errors_no[:,i],
      dtpts,errors_nd[:,i],[1e0,1],utf8("mom2.pdf"))
err_cpu(errors[:,i], cpu_times[:],
      errors_no[:,i],cpu_times_no[:],
      errors_nd[:,i],cpu_times_nd[:],[1e0,1], utf8("mom2_.pdf"))
#########################################
println("finished preparing plots")


text:
module SkipLists

using Base
import ..AbstractDataStructure

export SkipList, ins!, del!, at, search, len, ordered_data_structure_p

struct _EOL end
struct _HEAD end

Base.:(==)(::_EOL, ::Any) = false
Base.:(<)(::_EOL, ::Any) = false
Base.:(<)(::Any, ::_EOL) = true
Base.:(>)(::_EOL, ::Any) = true
Base.:(>)(::Any, ::_EOL) = false

Base.:(==)(::_HEAD, ::Any) = true
Base.:(<)(::_HEAD, ::Any) = true
Base.:(<)(::Any, ::_HEAD) = false
Base.:(>)(::_HEAD, ::Any) = false
Base.:(>)(::Any, ::_HEAD) = true

EOL = _EOL()
HEAD = _HEAD()

mutable struct LaneNode{V}
    data::Union{V, _HEAD, _EOL}
    node_width::Int  # Number of nodes at the base lane until the next node at this lane
    next::Union{LaneNode, Nothing}
    lower_lane_node::Union{LaneNode, Nothing}
    higher_lane_node::Union{LaneNode, Nothing}
    function LaneNode{V}(
            data::Union{V, _HEAD, _EOL},
            node_width::Int,
            next::Union{LaneNode{V}, Nothing},
            lower_lane_node::Union{LaneNode{V}, Nothing},
            higher_lane_node::Union{LaneNode{V}, Nothing},
    ) where {V}
        @assert((data === EOL && next === nothing && node_width == 0)
                || (data !== EOL && node_width > 0))
        new{V}(data, node_width, next, lower_lane_node, higher_lane_node)
    end
end

lower_lane_repr(l::LaneNode)::String = l.lower_lane_node===nothing ? "⤈" : "↓"
higher_lane_repr(l::LaneNode)::String = l.higher_lane_node===nothing ? "⤉" : "↑"
width_repr(l::LaneNode)::String = "w:$(l.node_width)"

# Lane iteration over all data nodes
Base.iterate(l::LaneNode) = ((l.data, width_repr(l), lower_lane_repr(l), higher_lane_repr(l)), l.next)
Base.iterate(::LaneNode, state::LaneNode) = ((
    state.data, width_repr(state), lower_lane_repr(state), higher_lane_repr(state)), state.next)
Base.iterate(::LaneNode, ::Nothing) = nothing

mutable struct SkipList{V} <: AbstractDataStructure.DataStructure{Int, V}
    p::Real  # The lane promotion probability for new nodes
    lanes::Vector{LaneNode{V}}  # where lanes[1] is the full lane
    length::Int
    function SkipList{V}() where {V}
        end_node = LaneNode{V}(EOL, 0, nothing, nothing, nothing)
        start_node = LaneNode{V}(HEAD, 1, end_node, nothing, nothing)
        lanes::Vector{LaneNode{V}} = [start_node]
        new{V}(0.5, lanes, 0)
    end
end

function SkipList{V}(vec::Array{Any, 1}) where {V}
    if length(vec) > 0
        vec::Vector{V} = convert(Array{V}, vec)
        return SkipList{V}(vec)
    end
    SkipList{V}()
end

function SkipList{V}(vec::Vector{V}) where {V}
    l = SkipList{V}()
    adjust_lane_count!(l, length(vec))
    for e in vec
        insert_nodes, _ = search_last_lt(l, e)
        create_new_node!(l.p, insert_nodes, e)
    end
    l.length = length(vec)
    l
end

function Base.setproperty!(obj::SkipList, attr::Symbol, rhs)
    if attr == :p
        @assert(0 < rhs < 1)
    end
    setfield!(obj, attr, rhs)
end

AbstractDataStructure.ordered_data_structure_p(::SkipList) = true

mutable struct SkipListPtr{V}
    skip_list::SkipList{V}
    lane::Int
    node::Union{LaneNode{V}, Nothing}
    bottom_idx::Int
end

function hopn(l::SkipList, n::Int)::Vector{LaneNode}
    @assert(n ≥ 0)
    result_vec = Vector{LaneNode}(undef, length(l.lanes))
    copyto!(result_vec, CartesianIndices(l.lanes), l.lanes, CartesianIndices(l.lanes))
    hopn(SkipListPtr(l, length(l.lanes), l.lanes[end], 0), result_vec, n)
end

function hopn(ptr::SkipListPtr, result_vec::Vector{LaneNode}, n::Int)::Vector{LaneNode}
    # Iterate at highest lane, accummulating node_width while remaining under n.
    # Then go to the next lane and repeat until index n is found.
    # Return the corresponding node at lowest lane and all predecesors at higher lanes.
    while ptr.bottom_idx + ptr.node.node_width ≤ n && ptr.node.data !== EOL
        ptr.bottom_idx += ptr.node.node_width
        ptr.node = ptr.node.next
        if ptr.node === nothing
            throw(KeyError(n))
        end
    end
    if ptr.bottom_idx == n
        while ptr.node !== nothing
            result_vec[ptr.lane] = ptr.node
            ptr.lane -= 1
            ptr.node = ptr.node.lower_lane_node
        end
        return result_vec
    end
    result_vec[ptr.lane] = ptr.node
    if ptr.lane == 1
        return result_vec
    end
    ptr.lane -= 1
    ptr.node = ptr.node.lower_lane_node
    hopn(ptr, result_vec, n)
end

"""
Find all lanes's last nodes with a value that is strictly lower than val, and the bottom-lane
index. If the search value is lower than any data in the list, it will be the vector of all
HEADs, and an index of 0.
"""
function search_last_lt(l::SkipList{V}, val::V)::Tuple{Vector{LaneNode{V}}, Int} where {V}
    if l.lanes[1].next === missing
        throw("Corrupted SkipList, missing a HEAD or EOL node.")
    end
    result_vec = Vector{LaneNode{V}}(undef, length(l.lanes))
    search_last_lt(SkipListPtr(l, length(l.lanes), last(l.lanes), 0), result_vec, val)
end

function search_last_lt(
        ptr::SkipListPtr{V},
        result_vec::Vector{LaneNode{V}},
        val::V
)::Tuple{Vector{LaneNode{V}}, Int} where {V}
    previous = undef
    while ptr.node.data < val
        previous = ptr.node
        ptr.node = ptr.node.next
        ptr.bottom_idx += previous.node_width
    end
    # We have overshot: backtrack
    result_vec[ptr.lane] = previous
    if ptr.lane > 1
        # and go one lane lower if possible
        ptr.lane -= 1
        ptr.node = previous.lower_lane_node
        ptr.bottom_idx -= previous.node_width
    else  # ptr.lane == 1
        # We are at the bottom lane, we have our result
        @assert(ptr.node.lower_lane_node === nothing && !(undef in result_vec))
        return result_vec, ptr.bottom_idx
    end
    search_last_lt(ptr, result_vec, val)
end

Base.iterate(l::SkipList) = iterate(l, l.lanes[1].next)  # Skip the HEAD
Base.iterate(::SkipList, state::LaneNode) = state.data === EOL ? nothing : (state.data, state.next)
Base.iterate(::SkipList, ::Nothing) = nothing

function add_lanes!(l::SkipList{V}, how_many::Int) where {V}
    head_top_node = last(l.lanes)
    eol_top_node = head_top_node
    while eol_top_node.data !== EOL; eol_top_node = eol_top_node.next end
    top_lane_idx = length(l.lanes)
    new_lanes_vector = Vector{LaneNode{V}}(undef, top_lane_idx + how_many)
    copyto!(new_lanes_vector, CartesianIndices(l.lanes), l.lanes, CartesianIndices(l.lanes))
    lower_lane_head = head_top_node
    lower_lane_eol = eol_top_node
    for i in top_lane_idx+1 : top_lane_idx+how_many
        eol = LaneNode{V}(EOL, 0, nothing, lower_lane_eol, nothing)
        lower_lane_eol.higher_lane_node = eol
        head = LaneNode{V}(HEAD, l.length+1, eol, lower_lane_head, nothing)
        lower_lane_head.higher_lane_node = head
        lower_lane_eol = eol
        lower_lane_head = head
        new_lanes_vector[i] = head
    end
    # Run stochastic promotion on all nodes sticking out of the previous top lane
    for node in head_top_node
        # TODO(agravier)
    end
    l.lanes = new_lanes_vector
end

function delete_lanes!(l::SkipList{V}, how_many::Int) where {V}
    cur_len = length(l.lanes)
    node = l.lanes[cur_len-how_many]
    while node !== nothing
        node.higher_lane_node = nothing
        node = node.next
    end
    l.lanes = l.lanes[1:cur_len-how_many]
end

function target_lane_count(lane_promotion_probability::Real, elem_count::Int)::Int
    try
        Int(ceil(log(1/lane_promotion_probability, elem_count))) + 1
    catch e
        if isa(e, InexactError)
            return 1  # we've got 0 element
        end
        throw(e)
    end
end

"""
elem_count: new number of elements
"""
function adjust_lane_count!(l::SkipList, elem_count::Int)
    lane_c = length(l.lanes)
    target_lane_c = target_lane_count(l.p, elem_count)
    δ = target_lane_c - lane_c
    if δ > 0
        add_lanes!(l, δ)
    elseif δ < 0
        delete_lanes!(l, abs(δ))
    end
end

function nhops(a::LaneNode, b::LaneNode)::Union{Int, Nothing}
    i = 0
    while a != b
        i += a.node_width
        a = a.next
        if a === nothing
            return nothing
        end
    end
    i
end

"""
`after` is a slice of the lanes after which the new node should be introduced.
This DOES NOT increase the `length` attribute of the skiplist.
This DOES NOT create new lanes when necessary.
"""
function create_new_node!(promotion_p::Real, after::Vector{LaneNode{V}}, val::V) where {V}
    lower_node = nothing
    promote = true
    for i in 1:length(after)
        if promote
            # Insert node: update all links and counts
            prev = after[i]
            next = prev.next
            lower_prev_node = prev.lower_lane_node
            lower_next_node = next.lower_lane_node
            lower_prev_width_to_anchor = lower_node === nothing ? 1 : nhops(lower_prev_node, lower_node)
            lower_next_width_from_anchor = lower_node === nothing ? 1 : nhops(lower_node, lower_next_node)
            new_node = LaneNode{V}(val, lower_next_width_from_anchor, next, lower_node, nothing)
            prev.next = new_node
            prev.node_width = lower_prev_width_to_anchor
            if lower_node !== nothing
                lower_node.higher_lane_node = new_node
            end
            lower_node = new_node
            promote = rand(Float64) < promotion_p
        else
            # No more node promotion: update upper lanes's previous node width
            after[i].node_width += 1
        end
    end
end

"""
ins!(ds::DataStructure, val::Value) -> Ref
"""
function AbstractDataStructure.ins!(l::SkipList{V}, val::V)::Int where {V}
    adjust_lane_count!(l, l.length+1)
    insert_nodes, bottom_prev_idx = search_last_lt(l, val)
    create_new_node!(l.p, insert_nodes, val)
    l.length += 1
    bottom_prev_idx
end

"""
del!(ds::DataStructure, pos::Ref) -> Value
raises KeyError
"""
function AbstractDataStructure.del!(l::SkipList{V}, pos::Int)::V where {V}
    @assert(pos > 0)
    node_vec = hopn(l, pos-1)
    if node_vec[1].next === nothing || node_vec[1].next.data === EOL
        throw(KeyError(pos))
    end
    adjust_lane_count!(l, l.length-1)
    val = node_vec[1].next.data
    height_at_pos = 1
    node = node_vec[1].next
    while node.higher_lane_node !== nothing
        height_at_pos += 1
        node = node.higher_lane_node
    end
    i = height_at_pos
    while i > 0
        node_vec[i].node_width += node_vec[i].next.node_width - 1
        node_vec[i].next = node_vec[i].next.next
        i -= 1
    end
    for node in node_vec[height_at_pos+1:end]
        node.node_width -= 1
    end
    l.length -= 1
    val
end

"""
at(ds::DataStructure, pos::Ref) -> Value
raises KeyError
"""
function AbstractDataStructure.at(l::SkipList{V}, pos::Int)::V where {V}
    @assert(pos > 0)
    node_vec = hopn(l, pos)
    if node_vec[1].data === EOL; throw(KeyError(pos)) end
    node_vec[1].data
end

"""
search(ds::DataStructure, val::Value) -> Vector{Ref}
"""
function AbstractDataStructure.search(ds::SkipList{V}, val::V)::Vector{Int} where {V}
    last_lt_node_vec, bottom_idx = search_last_lt(ds, val)
    node::LaneNode{V} = last_lt_node_vec[1].next
    result = Vector{Int}()
    while node.data == val
        push!(result, bottom_idx)
        bottom_idx += 1
        node = node.next
    end
    result
end

"""
len(ds::DataStructure) -> Int
"""
function AbstractDataStructure.len(ds::SkipList)::Int
    return ds.length
end

_iobuf = IOBuffer()
printstyled(IOContext(_iobuf, :color => true), "[ ", color=:green)
_LIST_OPEN_COLOR = String(take!(_iobuf))
printstyled(IOContext(_iobuf, :color => true), " ]", color=:green)
_LIST_CLOSE_COLOR = String(take!(_iobuf))
printstyled(IOContext(_iobuf, :color => true), " → ", color=:light_green)
_LIST_SEP_COLOR = String(take!(_iobuf))
_iobuf = nothing
MAX_LENGTH_DISPLAY_LIST = 10

function Base.show(io::IO, x::SkipList)
    if get(io, :color, false)
        open = _LIST_OPEN_COLOR
        close = _LIST_CLOSE_COLOR
        sep = _LIST_SEP_COLOR
    else
        open = "[ "
        close = " ]"
        sep = " → "
    end
    lane_c = length(x.lanes)
    el_word = "element" * (x.length>1 ? "s" : "")
    l_word = "lane" * (lane_c>1 ? "s" : "")
    x_disp = nothing
    if x.length ≤ MAX_LENGTH_DISPLAY_LIST
        x_disp = x
    else
        x_disp = Vector{Any}(undef, MAX_LENGTH_DISPLAY_LIST)
        for (i, e) in enumerate(x)
            x_disp[i] = e
            if i > MAX_LENGTH_DISPLAY_LIST - 1
                break
            end
        end
        x_disp[MAX_LENGTH_DISPLAY_LIST] = "…"
    end
    print(io, "SkipList with $(x.length) $el_word and $(lane_c) $l_word: " * open)
    join(io, x_disp, sep)
    print(io, close)
end

function Base.show(io::IO, ::_HEAD)
    print(io, "HEAD")
end

function Base.show(io::IO, ::_EOL)
    print(io, "EOL")
end

_iobuf = IOBuffer()
printstyled(IOContext(_iobuf, :color => true), "⟦ ", color=:red)
_LANE_NODE_OPEN_COLOR = String(take!(_iobuf))
printstyled(IOContext(_iobuf, :color => true), " ⟧", color=:red)
_LANE_NODE_CLOSE_COLOR = String(take!(_iobuf))
printstyled(IOContext(_iobuf, :color => true), " → ", color=:light_red)
_LANE_NODE_SEP_COLOR = String(take!(_iobuf))
_iobuf = nothing

function Base.show(io::IO, x::LaneNode)
    if get(io, :color, false)
        open = _LANE_NODE_OPEN_COLOR
        close = _LANE_NODE_CLOSE_COLOR
        sep = _LANE_NODE_SEP_COLOR
    else
        open = "⟦ "
        close = " ⟧"
        sep = " → "
    end
    print(io, open)
    join(io, x, sep)
    print(io, close)
end

function Base.show(io::IO, x::Vector{LaneNode})
    for i in length(x):-1:1
       println(io, "$i: $(x[i])")
    end
end

end  # module SkipLists


text:
function surround2D(b, i, j)
    h, w = size(b)
    [b[x,y] for x in i-1:i+1, y in j-1:j+1 if (0 < x <= h && 0 < y <= w)]
end

surroundhas1or2(b, i, j) = 0 < sum(map(x->Char(x)=='H', surround2D(b, i, j))) <= 2 ? 'H' : '.'

function boardstep!(currentboard, nextboard)
    x, y = size(currentboard)
    for j in 1:y, i in 1:x
        ch = Char(currentboard[i, j])
        if ch == ' '
            continue
        else
            nextboard[i, j] = (ch == 'H') ? 't' : (ch == 't' ? '.' :
                   surroundhas1or2(currentboard, i, j))
        end
    end
end

const b1 = "         " *
           "  tH     " *
           " .  .... " *
           "  ..     " *
           "         "
const mat = reshape(map(x->UInt8(x[1]), split(b1, "")), (9, 5))'
const mat2 = copy(mat)

function printboard(mat)
    for i in 1:size(mat)[1]
        println("\t", join([Char(c) for c in mat[i,:]], ""))
    end
end

println("Starting Wireworld board:")
printboard(mat)
for step in 1:8
    boardstep!(mat, mat2)
    println(" Step $step:")
    printboard(mat2)
    mat .= mat2
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using Weave

```
Converts `exercise.jl` into `exercise_skeleton.ipynb`.
```
function _generate_skeleton(file_name::AbstractString)
    lines = readlines(file_name, keep=true)
    skeleton_file = replace(file_name, ".jl" => "_skeleton.jl")

    open(skeleton_file, "w+") do f
        for line in lines
            startswith(strip(line), "#%") ? write(f, "\n") : write(f, line)
        end
    end

    convert_doc(skeleton_file, replace(skeleton_file, ".jl" => ".ipynb"))
    rm(skeleton_file)
end

```
Converts a file `exercise_complete.ipynb` into `exercise.jl`. Use this to generate
.jl scripts which will have the solution lines marked with '#%'.

Assumes that inside each code cell, the parts that should be removed from the skeleton are
separated by lines corresponding to the parameter _keyword_.
```
function _generate_script_from_solution(file_name::AbstractString; keyword="#%")
    doc = Weave.WeaveDoc(file_name)
    for chunk in filter(chunk -> typeof(chunk) == Weave.CodeChunk, doc.chunks)
        flag = false
        new_content = ""
        for line in split(chunk.content,"\n")
            if flag
                if strip(line) == keyword
                    flag = false
                else
                    new_content *= "#% " * line * "\n"
                end
            else
                if strip(line) == keyword
                    flag = true
                    new_content *= "# TODO: add your code here\n"
                else
                    new_content *= line * "\n"
                end
            end
        end
        if flag
            @warn("Unmatched " * keyword * " in code cell number $(chunk.number), cell is passed through as is!")
        else
            chunk.content = new_content
        end
    end

    converted = Weave._convert_doc(doc, "script")

    outfile = replace(file_name, ".ipynb" => ".jl")
    outfile = replace(outfile, "_complete" => "")
    open(outfile, "w") do f
        write(f, converted)
    end
    return outfile
end

```
Generates both skeleton and complete notebooks.
```
function generate_skeleton(file_name::AbstractString)
    jl_file = _generate_script_from_solution(file_name)
    _generate_skeleton(jl_file)
    rm(jl_file)
end


text:
#Adapted from https://github.com/github/developer.github.com/blob/master/Rakefile#L21-L48

cd(@__DIR__)

last_commit=readchomp(`git --no-pager log -1 --pretty=format:"%h:%s"`)

ENV["GIT_DIR"]=abspath(chomp(readall(`git rev-parse --git-dir`)))

old_sha = chomp(readall(`git rev-parse refs/remotes/origin/gh-pages`))

run(`jekyll build`)

cd("_site") do

	gif="/tmp/dev.gh.i"
	ENV["GIT_INDEX_FILE"]=gif
	ENV["GIT_WORK_TREE"]=pwd()
	run(`git add -A`)
	tsha=chomp(readall(`git write-tree`))
	mesg="Deploy docs for master@$last_commit"

	if length(old_sha) == 40
	  	csha = chomp(readall(`git commit-tree $tsha -p $old_sha -m $(mesg)`))
	else
		csha = chomp(readall(`git commit-tree $tsha -m $(mesg)`))
	end

	 print("Created commit $csha")

	 run(`git --no-pager show $csha --stat`)

	 run(`git update-ref refs/heads/gh-pages $csha `)

	 run(`git push origin gh-pages `)

end


text: a = 10


text:
using Distributions
using SpecialFunctions
using JLD2
using FileIO
using Distributed
using LaTeXStrings
using Glob
using QuadGK

#This has the exponent (alpha d)^phi

function norm_gen_gaus(ϕ,α)
    return (2 /α * gamma(1+ 1/ϕ))
end

function integral_KB_del_ph(δ,ϕ,α)
    integral, err = quadgk(x -> exp(-α^(ϕ) *(abs(x))^ϕ)* (-α^(ϕ)* (abs((x)^ϕ) -abs((x-δ)^ϕ) )), - Inf, Inf, rtol=1e-3)
    return integral/norm_gen_gaus(ϕ,α)
end

function takememory_phi(β,v,α,ϕ)
    prob =probmemory_phi(β,v,α,ϕ)
    wahrscheinlikeit = rand()
    result=false
    if wahrscheinlikeit < prob
        result=true
    end
    result
end

function probmemory_phi(β,v,α,ϕ)
    prob = 1- exp(-β *energiememory_phi(v,α,ϕ) )
    prob
end

function energiememory_phi(v,α,ϕ)
    energie =exp(-(α^(ϕ)*(abs(v))^ϕ))/norm_gen_gaus(ϕ,α)
end

function freeenergy_real_phi(β,v,α,ϕ)
    energy = log( exp(β *energiememory_phi(v,α,ϕ)) -1)/β
    energy
end

function dissipation_phi(β,v,α,ϕ)
    diss = integral_KB_del_ph(v,ϕ,α)/β
    diss
end

function netUtility_phi(β,v,α,ϕ)
    unet=freeenergy_real_phi(β,v,α,ϕ)- dissipation_phi(β,v,α,ϕ)
    unet
end

function newv_drift(δ,n)
    vnew = sqrt(n)δ+randn()*δ*0.05*sqrt(n)
    #vnew = sqrt(n)δ+randn()*δ*0.1
    vnew
end

function onestep_dritf_real_phi(β,δ,oldv,α,n,ϕ)

    if (takememory_phi(β,oldv,α,ϕ))
        free=freeenergy_real_phi(β,oldv,α,ϕ)
        diss=dissipation_phi(β,oldv,α,ϕ)
        unet=free -diss
        nextv=newv_drift(δ,n)
        n+=1
        naiveuse=0
    else
        unet = 0
        free=0
        diss=0
        nextv=newv_drift(δ,1)
        naiveuse=1
        n=2
    end

    return (unet,nextv,naiveuse,free,diss,n)
end

function livesycle_drift_real_phi(β,δ,α,meanencounters,reps,ϕ)
    meanunet=0 #initializing for all the rounds for fist cost
    meanfree=0
    meandiss=0
    naiveuses = 0
    for round in 1:reps

        localnetu= 0
        netu = 0
        localfree= 0
        free = 0
        localdiss= 0
        diss = 0

        v=newv_drift(δ,1) #start the cirus at the memory wich is allways at 0
        naiveuses += 1/reps
        n=2
        for encounter in 2:meanencounters  #start at 1 because initialization is allways at 0 in the first round without any prior memory
             naiveuse =0    # Used naive for the first virus

             (netu,v,naiveuse,free,diss,n)=onestep_dritf_real_phi(β,δ,v,α,n,ϕ)
            localnetu+=netu
            localfree+=free
            localdiss+=diss
            naiveuses += naiveuse/reps
        end

        meanunet+= localnetu/reps
        meanfree+= localfree/reps
        meandiss+= localdiss/reps
    end

    return (meanunet,naiveuses,meanfree,meandiss)
end

function findoptalphaoflive_drift_real_phi(β,δ,meanencounters,reps,αsteps,αmax,ϕ)

    start= 0.01
    α = range(start,stop=αmax,length=αsteps)

    data=livesycle_drift_real_phi.(β,δ,α,meanencounters,reps,ϕ)
    usesmaive=zeros(αsteps)
    memuses=zeros(αsteps)
    utility=zeros(αsteps)
    free=zeros(αsteps)
    diss=zeros(αsteps)
    for ii in 1:αsteps
        utility[ii] = data[ii][1]
        usesmaive[ii] = data[ii][2]
        free[ii] = data[ii][3]
        diss[ii] = data[ii][4]
    end
    memuses= meanencounters./usesmaive

    amax= argmax(utility)
    optalpha = α[amax]
    Nusesatmax=memuses[amax]
    argmaxn = argmax(memuses)
    maxNpossible=memuses[argmaxn]
    utiliymax = utility[amax]
    freemax=free[amax]
    dissmax=diss[amax]

    return Nusesatmax,utiliymax,optalpha,maxNpossible,freemax,dissmax
end

function scandeltaliveGeqBparallel_drift_real_phi(β,meanencounters,reps,αsteps,αmax,deltasteps,ϕ)

    δ = range(0.005,stop=1.5/αmax,length= deltasteps)
    memstepsatmax =zeros(deltasteps)
    utilityopt =zeros(deltasteps)
    alphaopt =zeros(deltasteps)
    posnmax =zeros(deltasteps)
    freemax =zeros(deltasteps)
    dissmax =zeros(deltasteps)
    Threads.@threads for i =1:deltasteps
        memstepsatmax[i],utilityopt[i],alphaopt[i],posnmax[i],freemax[i],dissmax[i]=findoptalphaoflive_drift_real_phi(β,δ[i],meanencounters,reps,αsteps,αmax,ϕ)
    end
    return memstepsatmax,utilityopt,alphaopt,posnmax,δ,freemax,dissmax

end

function scanbetaeqGdeltaliveparallel_drift_real_phi(meanencounters,reps,αsteps,αmax,deltasteps,betasteps,ϕ)

    β = range(0.1,length=betasteps,stop=10.0)

    data = scandeltaliveGeqBparallel_drift_real_phi.(β,meanencounters,reps,αsteps,αmax,deltasteps,ϕ)

    memstepsatmax= zeros(betasteps,deltasteps)
    utilityopt= zeros(betasteps,deltasteps)
    alphaopt= zeros(betasteps,deltasteps)
    posnmax= zeros(betasteps,deltasteps)
    freemax =zeros(betasteps,deltasteps)
    dissmax =zeros(betasteps,deltasteps)
    for ii in 1:betasteps
        memstepsatmax[ii,:]=data[ii][1]
        utilityopt[ii,:]=data[ii][2]
        alphaopt[ii,:] = data[ii][3]
        posnmax[ii,:]=data[ii][4]
        freemax[ii,:]=data[ii][6]
        dissmax[ii,:]=data[ii][7]

    end
    δ = data[1][5]

    introduction="This gives the data meory_used_at_the_max  utility_atmax  otimal_alpha posible_maximum_memoryuses delta and beta and free_at_the_max and diss_at_the_max the data is Gamma first and delta second index"
 return memstepsatmax,utilityopt,alphaopt,posnmax,δ,β,freemax,dissmax,introduction

end

####Data Production

function simulation_DM_cluster_drift_real_oneB_phi(meanencounters,reps,asteps,bsteps,dsteps,amax,thisb,ϕ)

    βas = range(0.1,length=bsteps,stop=10.0)
    β=βas[thisb]
    parameters = (meanencounters,reps,asteps,bsteps,dsteps,amax,β,ϕ)

    filename = "Data/Gen_Bell/Singel_beta/meanencounters_$(meanencounters)/reps_$(reps)/Phi_$(ϕ)/DATA_Asteps_$(asteps)-Bsteps_$(bsteps)-Dsteps_$(dsteps)-amax_$(amax)_real-thisb_$(thisb).jld2"
    println("I want to do  $(filename)")
    if isfile(filename)
            println("Data $(filename) allready exists")
    else
        if !isdir("Data")
            mkdir("Data")
        end
        if !isdir("Data/Gen_Bell")
            mkdir("Data/Gen_Bell")
        end
        if !isdir("Data/Gen_Bell/Singel_beta")
            mkdir("Data/Gen_Bell/Singel_beta")
        end

        if !isdir("Data/Gen_Bell/Singel_beta/meanencounters_$(meanencounters)")
            mkdir("Data/Gen_Bell/Singel_beta/meanencounters_$(meanencounters)")
        end
        if !isdir("Data/Gen_Bell/Singel_beta/meanencounters_$(meanencounters)/reps_$(reps)")
            mkdir("Data/Gen_Bell/Singel_beta/meanencounters_$(meanencounters)/reps_$(reps)")
        end
        if !isdir("Data/Gen_Bell/Singel_beta/meanencounters_$(meanencounters)/reps_$(reps)/Phi_$(ϕ)")
            mkdir("Data/Gen_Bell/Singel_beta/meanencounters_$(meanencounters)/reps_$(reps)/Phi_$(ϕ)")
        end

        memstepsatmax,utilityopt,alphaopt,posnmax,δ,freemax,dissmax = scandeltaliveGeqBparallel_drift_real_phi(β,meanencounters,reps,asteps,amax,dsteps,ϕ)

    #filename = "/Users/oskar/Dropbox/PHD/Simulations/b-cell/Simulation/src/julia/data/test.jld"
        save(filename,"meory_used_at_the_min",  memstepsatmax,"max_utility",utilityopt,"otimal_alpha",alphaopt,"posible_maximum_memoryuses",posnmax,"parameters",parameters,"beta",β,"delta",δ,"freemax",freemax,"dissmax",dissmax,"thisb",thisb)
        println("I'm saving the file $(filename)")
    end
end

function addup_simulation_DM_cluster_drift_real_oneB_phi(meanencounters,reps,asteps,bsteps,dsteps,amax,ϕ)

    filename = "Data/Gen_Bell/full/meanencounters_$(meanencounters)/reps_$(reps)/Phi_$(ϕ)/DATA_Asteps_$(asteps)-Bsteps_$(bsteps)-Dsteps_$(dsteps)-amax_$(amax)_real.jld2"
    deltas = range(0.0001,length=dsteps,stop=0.4)

    parameters = (meanencounters,reps,asteps,bsteps,dsteps,amax,ϕ)

    if !isdir("Data/Gen_Bell/full")
        mkdir("Data/Gen_Bell/full")
    end
        if !isdir("Data/Gen_Bell/full/meanencounters_$(meanencounters)")
            mkdir("Data/Gen_Bell/full/meanencounters_$(meanencounters)")
        end
        if !isdir("Data/Gen_Bell/full/meanencounters_$(meanencounters)/reps_$(reps)")
            mkdir("Data/Gen_Bell/full/meanencounters_$(meanencounters)/reps_$(reps)")
        end
        if !isdir("Data/Gen_Bell/full/meanencounters_$(meanencounters)/reps_$(reps)/Phi_$(ϕ)")
            mkdir("Data/Gen_Bell/full/meanencounters_$(meanencounters)/reps_$(reps)/Phi_$(ϕ)")
        end

        βas = range(0.1,length=bsteps,stop=10.0)

    direct= "Data/Gen_Bell/Singel_beta/meanencounters_$(meanencounters)/reps_$(reps)/Phi_$(ϕ)"
    filenames = "DATA_Asteps_$(asteps)-Bsteps_$(bsteps)-Dsteps_$(dsteps)-amax_$(amax)_real-thisb_*.jld2"

    files2 = get_files(direct,filenames)
    l=length(files2)
    if l<bsteps
        println("!!!!  Achtung !!!!!!!!!!")
        println("es gibt nicht alle files ")
        println("!!!!  Achtung !!!!!!!!!!")
        println(" ")
        println("I have $(l)  of $(bsteps)")
        println(" ")
        println("Please do again ")
        println(" meanencounters_$(meanencounters)/reps_$(reps)/Phi_$(ϕ)/DATA_Asteps_$(asteps)-Bsteps_$(bsteps)-Dsteps_$(dsteps)-amax_$(amax)")
        #return -1
    end

    memstepsatmax= zeros(bsteps,dsteps)
    utilityopt= zeros(bsteps,dsteps)
    alphaopt= zeros(bsteps,dsteps)
    posnmax= zeros(bsteps,dsteps)
    freemax =zeros(bsteps,dsteps)
    dissmax =zeros(bsteps,dsteps)
    for ii in 1:length(files2)
        thisb = load(files2[ii],"thisb")
        memstepsatmax[thisb,:]=load(files2[ii],"meory_used_at_the_min")
        utilityopt[thisb,:]=load(files2[ii],"max_utility")
        alphaopt[thisb,:] = load(files2[ii],"otimal_alpha")
        posnmax[thisb,:]=load(files2[ii],"posible_maximum_memoryuses")
        freemax[thisb,:]=load(files2[ii],"freemax")
        dissmax[thisb,:]=load(files2[ii],"dissmax")
    end

    δ = load(files2[1],"delta")

    #filename = "/Users/oskar/Dropbox/PHD/Simulations/b-cell/Simulation/src/julia/data/test.jld"
        save(filename,"meory_used_at_the_min",  memstepsatmax,"max_utility",utilityopt,"otimal_alpha",alphaopt,"posible_maximum_memoryuses",posnmax,"parameters",parameters,"beta",βas,"delta",δ,"freemax",freemax,"dissmax",dissmax,"phi",ϕ)
            println("I'm saving the file $(filename)")

end

function get_files(direct,filenames)
    filenames=glob(filenames,direct)
    return filenames
end

max_stars_repo_name: StatPhysBio/ImmuneMemoryDM

max_issues_repo_name: StatPhysBio/ImmuneMemoryDM

max_forks_repo_name: StatPhysBio/ImmuneMemoryDM


text:
using PartialLeastSquaresRegressor
using Test
using MLJBase
using Statistics
using Random

# load some locally defined models needed for testing:
include("./_standardizer.jl")
import .Stand

include("./utils_test.jl")
include("./pls1_test.jl")
include("./pls2_test.jl")
include("./kpls_test.jl")


text:
module Command

using JLD2, FileIO
using LinearAlgebra, SpecialFunctions, NMF
using Distributions, Statistics, StatsBase
using Plots

include("../src/scrna.jl")
include("../src/mle.jl")
include("../src/util.jl")
include("../src/pointcloud.jl")
include("../src/SeqSpace.jl")

using .scRNA, .SeqSpace

# ------------------------------------------------------------------------
# variable inputs

mutable struct Parameters
    name      :: String
    threshold :: NamedTuple{(:gene,:cell), Tuple{Float64,Float64}}
    subdir    :: Union{Nothing, AbstractString}
    filter    :: Function
    plot      :: Bool
    δ         :: Int
end

const default = Parameters(
    "",
    (gene = 5e-3, cell = 2e-1),
    nothing,
    (X) -> X,
    true,
    0
)

Parameters(name;
    threshold=default.threshold,
    subdir=default.subdir,
    filter=default.filter,
    plot=default.plot,
    δ=default.δ
) = Parameters(name, threshold, subdir, filter, plot, δ)

include("utils.jl")
using .CommandUtility

# ------------------------------------------------------------------------
# plotting code

module Plot

using Statistics, StatsBase
using Plots, ColorSchemes

rank(x) = invperm(sortperm(x))

cdf(x;  kwargs...) = plot(sort(x),  range(0,1,length=length(x)); kwargs...)
cdf!(x; kwargs...) = plot!(sort(x), range(0,1,length=length(x)); kwargs...)

function marginals(matrix; ϵ=1e-6)
    p₁ = cdf(vec(mean(matrix,dims=1)).+ϵ, xscale=:log10, label="", linewidth=2)
    xaxis!("mean count/cell")
    yaxis!("CDF")

    p₂ = cdf(vec(mean(matrix,dims=2)).+ϵ, xscale=:log10, label="", linewidth=2)
    xaxis!("mean count/gene")
    yaxis!("CDF")

    return plot(p₁, p₂)
end

function qq(x, model, param)
    p = scatter(rank(x)./length(x), model.cumulative(param), linewidth=2, label="")
    plot!(0:1, 0:1, color=:red, linewidth=2, linestyle=:dashdot, label="ideal")
    xaxis!("empirical quantile")
    yaxis!("model quantile")

    return p
end

function mlefits(matrix, param)
    mₓ = vec(mean(matrix, dims=2))
    Mₓ = vec(maximum(matrix, dims=2))

    p₁ = scatter(mₓ, param.α,
        alpha=0.1,
        xscale=:log10,
        marker_z=log10.(Mₓ),
        label=false
    )
    xaxis!("expression / cell")
    yaxis!("estimated α")

    p₂ = scatter(mₓ, param.β,
        alpha=0.1,
        xscale=:log10,
        marker_z=log10.(Mₓ),
        label=false
    )
    xaxis!("expression / cell")
    yaxis!("estimated β")

    p₃ = scatter(mₓ, param.γ,
        alpha=0.1,
        xscale=:log10,
        yscale=:log10,
        marker_z=log10.(Mₓ),
        label=false
    )
    xaxis!("expression / cell")
    yaxis!("estimated γ")

    return p₁, p₂, p₃
end

function mpmaxeigval(λ, x₀, rank)
    p = plot(λ, 1:length(λ),
             xscale=:log10,
             yscale=:log10,
             linewidth=2,
             label="Empirical distribution"
    )

    vline!([x₀], linewidth=2, linestyle=:dashdot, label="MP maximum eigenvalue")
    title!("scRNAseq data (k=$rank)")
    xaxis!("singular value")
    yaxis!("CDF")

    p
end

function ipr(ψ)
    IPR = vec(sum(ψ.^2,dims=2).^2 ./ sum(ψ.^4,dims=2))
    p = plot(1:length(IPR), IPR,
             xscale=:log10,
             xlabel="principal component",
             ylabel="participation ratio",
             yscale=:log10,
             label="",
    )
    return p
end

end

# ------------------------------------------------------------------------
# the bulk of the code

function data(dir::AbstractString, subdir, filter::Function)
    N = if subdir === nothing
        scRNA.load(dir)
    else
        reduce(∪, scRNA.load("$dir/$d") for d ∈ readdir(dir) if occursin(subdir,d))
    end

    return filter(N)
end

function cutoff(counts, threshold)
    counts = scRNA.filtergene(counts) do gene, _
        mean(gene) >= threshold.gene && length(unique(gene)) > 3
    end

    counts = scRNA.filtercell(counts) do cell, _
        mean(cell) >= threshold.cell
    end

    return counts
end

function normalize(dir::AbstractString, param::Parameters, figs::AbstractString)
    # load in raw data. impose data specific filter
    alert("-->loading raw data")
    counts = data(dir, param.subdir, param.filter)
    param.plot && savefig(Plot.marginals(counts), "$figs/raw_count_marginals.png")

    # impose (conservative) count-based cutoffs
    alert("-->threshold raw data")
    counts = cutoff(counts, param.threshold)
    param.plot && savefig(Plot.marginals(counts), "$figs/threshold_count_marginals.png")

    # estimate marginalized overdispersion distribution from highly expressed genes
    alert("-->estimating overdispersion prior")
    p₀ = MLE.fit_glm(:negative_binomial, counts;
        Γ=(β̄=1, δβ¯²=10, Γᵧ=nothing),
        run=(x) -> mean(x) > 1
    )
    param.plot && savefig(
        Plot.cdf(p₀.γ, xlabel="estimated γ", ylabel="CDF", xscale=:log10, label="", linewidth=2),
        "$figs/overdispersion_distribution_estimate.png"
    )

    logγ  = log.(p₀.γ)
    model = MLE.generalized_normal(logγ)
    prior = MLE.fit(model)
    param.plot && savefig(Plot.qq(logγ, model, prior), "$figs/overdispersion_distribution_fit.png")

    # use overdispersion distribution as prior into the full fit
    alert("-->fitting overdispersion per gene")
    p₀ = MLE.fit_glm(:negative_binomial, counts; Γ=(β̄=1, δβ¯²=10, Γᵧ=prior))
    param.plot && let
        p = Plot.mlefits(counts, p₀)
        savefig(p[1], "$figs/nb_α_fit.png")
        savefig(p[2], "$figs/nb_β_fit.png")
        savefig(p[3], "$figs/nb_γ_fit.png")
    end

    # normalize variance of count matrix. estimate rank
    alert("-->normalizing count variance & estimating rank")
    N, σ², u², v² = let
        σ² = counts.*(counts.+p₀.γ) ./ (1 .+ p₀.γ)
        u², v², _ = Utility.sinkhorn(σ²)

        (Diagonal(.√u²) * counts * Diagonal(.√v²)), (Diagonal(u²) * σ² * Diagonal(v²)), u², v²
    end

    Λ = svd(N)
    λ = Λ.S
    R = sum(λ .> (sqrt(size(N,1))+sqrt(size(N,2)))) + param.δ
    param.plot && let
        savefig(Plot.mpmaxeigval(λ, √(size(N,1)) + √(size(N,2)), R), "$figs/rank_estimation.png")
        savefig(Plot.ipr(Λ.Vt), "$figs/participation_ratio.png")
    end

    # reduce rank to signal directions
    N, c = let
        r = nnmf(N, R; alg=:cd)
        m = r.W*r.H
        m, cor(m[:], N[:])
    end

    # renormalize
    alert("-->normalizing reduced rank matrix")
    N, u, v = let
        u, v, _ = Utility.sinkhorn(N)
        (Diagonal(u) * N * Diagonal(v)), u, v
    end

    return (
        raw       = counts.data,
        gene      = counts.gene,
        cell      = counts.cell,
        rank      = R,
        corr      = c,
        data      = N,
        varnorm   = (row=u², col=v²),
        scalenorm = (row=u,  col=v),
        mleparam  = p₀,
    )
end

# ------------------------------------------------------------------------
# main point of entry

if abspath(PROGRAM_FILE) == @__FILE__
    params = (
        p = Arg("inpath", loadparams, default),
        o = Arg("outpath", (path)->path, nothing)
    )

    args = argparse(ARGS, params)
    length(args) > 1 && error("too many input paths")

    input = args[1]
    !isdir(input) && error("directory $input not found")

    outdir = dirname(params.o.value)
    !isdir(outdir) && mkdir(outdir)

    figdir = "$(outdir)/figs/norm"
    !isdir(figdir) && mkpath(figdir)

    ENV["GKSwstype"] = "nul"
    allparams = isa(params.p.value, AbstractArray) ? params.p.value : [ params.p.value ]

    jldopen("$outdir/norms.jld2", "w") do data
        for (i,param) in enumerate(allparams)
            name = length(param.name) > 0 ? param.name : "set$i"
            alert("-->processing group $(name)")

            figsubdir = "$figdir/$(name)"
            !isdir(figsubdir) && mkpath(figsubdir)
            result = normalize(input, param, figsubdir)

            data["$name/normparam"] = param
            for (key,val) in pairs(result)
                data["$name/$key"] = val
            end
        end
    end
end

end

max_stars_repo_name: nnoll/seqspace

max_issues_repo_name: nnoll/seqspace

max_forks_repo_name: nnoll/seqspace


text:
module FisherRaoDistance

using KernelDensityEstimate
using MultivariateStats
using StatsBase

export
    get_low_dim_points,
    calc_vec_dist,
    fisherraodistance,
    fisherraotest,
    transform, #from MultivariateStats
    kde!, #from package KernelDensityEstimate
    evaluateDualTree,
    sample

# source files
include("dfr_functions.jl")

end #module


text:
import Base: show

mutable struct Indicator
    fun::Function
    paramset::ParameterSet
    data::TS
    function Indicator(fun::Function, paramset::ParameterSet, data::TS)
        return new(fun, paramset, data)
    end
end

function Indicator(fun::Function, paramset::ParameterSet)
    return Indicator(fun, paramset, TS())
end

function calculate(indicator::Indicator, input::TS; arg_values::Vector)::TS
    return indicator.fun(input; generate_dict(indicator.paramset, arg_values)...)
end

# TODO: add information about the calculation function
function show(io::IO, indicator::Indicator)
    show(io, indicator.paramset)
end

max_stars_repo_name: lithiferous/Strathreads

max_issues_repo_name: lithiferous/Strathreads

max_forks_repo_name: lithiferous/Strathreads


text:
#!/bin/env julia

function first_part(freqs)
    total = 0

    for f in split(freqs)
        num = parse.(Int, f)
        total += num
    end

    println(total)
end

function second_part(freqs)
    total = 0
    seen = Set([0])
    twice = -1

    while twice <= -1
        for f in split(freqs)
            num = parse.(Int, f)
            total += num

            if total in seen
                twice = total
                break
            end

            push!(seen, total)
        end
    end

    println(twice)
end

file = open("../input/input.txt")
freqs = read(file, String)

first_part(freqs)
second_part(freqs)

close(file)

max_stars_repo_name: GrayJack/Advent-Of-Code

max_issues_repo_name: GrayJack/Advent-Of-Code

max_forks_repo_name: GrayJack/Advent-Of-Code


text:
using RandomMatrixDistributions
using Random: seed!
using Test

@testset "RandomMatrixDistributions" begin
    seed!(1)
    include("MatrixSamplers.jl")
    include("Densities.jl")
end


text:
ia> [RpcApi.SimpleOrganism(0x000000000000c617,:(function (o::Creature.Organism,)
     local var_7::Int64 = -2881162876892958953
     local var_1::Int8 = 121
     Creature.eatUp(o,Int(var_1))
     for i = 1:var_1
         local var_8::Int16 = 23832
         local var_6::Int8 = 66
         Creature.eatRight(o,Int(var_6))
         Creature.eatRight(o,Int(var_6))
         Creature.eatRight(o,Int(var_6))
         Creature.eatRight(o,Int(var_6))
         Creature.eatRight(o,Int(var_6))
         Creature.eatRight(o,Int(var_6))
         Creature.eatDown(o,Int(var_6))
         Creature.eatRight(o,Int(var_6))
         Creature.eatRight(o,Int(var_6))
         Creature.eatUp(o,Int(var_6))
         var_8 = if haskey(o.mem,var_8)
                 o.mem[var_8]
             else
                 var_8
             end
         Creature.stepRight(o)
         Creature.stepRight(o)
     end
     Creature.eatUp(o,Int(var_1))
     Creature.stepLeft(o)
     Creature.stepLeft(o)
     for i = 1:var_1
         Creature.stepDown(o)
         Creature.stepLeft(o)
         Creature.stepRight(o)
     end
     Creature.stepUp(o)
     var_7 = var_7 / var_7
     Creature.stepLeft(o)
     Creature.stepLeft(o)
     Creature.stepUp(o)
     return true
 end),[100,300,95,0,1,1,1],1,5000,1,12758650,0x01070550,Dict{Int16,Int16}(),[446,333],31),RpcApi.SimpleOrganism(0x000000000000c498,:(function (o::Creature.Organism,)
     local var_7::Int64 = -2881162876892958953
     local var_1::Int8 = 121
     Creature.eatUp(o,Int(var_1))
     for i = 1:var_1
         local var_8::Int16 = 23832
         local var_6::Int8 = 66
         Creature.eatRight(o,Int(var_6))
         Creature.eatRight(o,Int(var_6))
         Creature.eatRight(o,Int(var_6))
         Creature.eatRight(o,Int(var_6))
         Creature.eatRight(o,Int(var_6))
         Creature.eatRight(o,Int(var_6))
         Creature.eatDown(o,Int(var_6))
         Creature.eatRight(o,Int(var_6))
         Creature.eatRight(o,Int(var_6))
         Creature.eatUp(o,Int(var_6))
         var_8 = if haskey(o.mem,var_8)
                 o.mem[var_8]
             else
                 var_8
             end
         Creature.stepRight(o)
         Creature.stepRight(o)
     end
     Creature.eatUp(o,Int(var_1))
     Creature.stepLeft(o)
     Creature.stepLeft(o)
     for i = 1:var_1
         Creature.stepDown(o)
         Creature.stepLeft(o)
         Creature.stepRight(o)
     end
     Creature.stepUp(o)
     var_7 = var_7 / var_7
     Creature.stepLeft(o)
     Creature.stepLeft(o)
     Creature.stepUp(o)
     return true
 end),[100,300,95,0,1,1,1],1,5000,1,12703705,0x01070550,Dict{Int16,Int16}(),[562,334],31),RpcApi.SimpleOrganism(0x000000000000c3c7,:(function (o::Creature.Organism,)
     local var_7::Int64 = -2881162876892958953
     local var_1::Int8 = 121
     Creature.eatUp(o,Int(var_1))
     for i = 1:var_1
         local var_8::Int16 = 23832
         local var_6::Int8 = 66
         Creature.eatRight(o,Int(var_6))
         Creature.eatRight(o,Int(var_6))
         Creature.eatRight(o,Int(var_6))
         Creature.eatRight(o,Int(var_6))
         Creature.eatRight(o,Int(var_6))
         Creature.eatRight(o,Int(var_6))
         Creature.eatDown(o,Int(var_6))
         Creature.stepRight(o)
         Creature.eatRight(o,Int(var_6))
         Creature.eatUp(o,Int(var_6))
         var_8 = if haskey(o.mem,var_8)
                 o.mem[var_8]
             else
                 var_8
             end
         Creature.stepRight(o)
         Creature.stepRight(o)
     end
     Creature.eatUp(o,Int(var_1))
     Creature.stepLeft(o)
     Creature.stepLeft(o)
     for i = 1:var_1
         Creature.stepDown(o)
         Creature.stepLeft(o)
         Creature.stepRight(o)
     end
     Creature.stepUp(o)
     var_7 = var_7 / var_7
     Creature.stepLeft(o)
     Creature.stepLeft(o)
     Creature.stepUp(o)
     return true
 end),[100,300,95,0,1,1,1],1,5000,1,12687734,0x0106ef78,Dict{Int16,Int16}(),[385,333],31),RpcApi.SimpleOrganism(0x000000000000c5f1,:(function (o::Creature.Organism,)
     local var_7::Int64 = -2881162876892958953
     local var_1::Int8 = 121
     for i = 1:var_1
         local var_8::Int16 = 23832
         local var_6::Int8 = 66
         Creature.eatRight(o,Int(var_6))
         Creature.eatRight(o,Int(var_6))
         Creature.eatRight(o,Int(var_6))
         Creature.eatLeft(o,Int(var_6))
         Creature.eatRight(o,Int(var_6))
         Creature.stepUp(o)
         Creature.eatDown(o,Int(var_6))
         Creature.stepRight(o)
         Creature.eatRight(o,Int(var_6))
         Creature.eatUp(o,Int(var_6))
         var_8 = if haskey(o.mem,var_8)
                 o.mem[var_8]
             else
                 var_8
             end
         Creature.eatDown(o,Int(var_6))
         Creature.stepRight(o)
     end
     Creature.eatUp(o,Int(var_1))
     Creature.stepLeft(o)
     Creature.stepLeft(o)
     for i = 1:var_1
         Creature.stepDown(o)
         Creature.stepLeft(o)
         Creature.stepRight(o)
     end
     Creature.stepUp(o)
     var_7 = var_7 / var_7
     Creature.stepLeft(o)
     Creature.stepLeft(o)
     var_1 = var_1 * var_1
     return true
 end),[100,300,95,0,1,1,1],1,5000,1,11431879,0x01070550,Dict{Int16,Int16}(),[637,121],30),RpcApi.SimpleOrganism(0x000000000000c3c6,:(function (o::Creature.Organism,)
     local var_7::Int64 = -2881162876892958953
     local var_1::Int8 = 121
     Creature.eatUp(o,Int(var_1))
     Creature.eatLeft(o,Int(var_1))
     for i = 1:var_1
         local var_8::Int16 = 23832
         local var_6::Int8 = 66
         Creature.eatRight(o,Int(var_6))
         Creature.eatRight(o,Int(var_6))
         Creature.eatLeft(o,Int(var_6))
         Creature.eatRight(o,Int(var_6))
         Creature.stepUp(o)
         Creature.stepRight(o)
         Creature.stepUp(o)
         Creature.stepLeft(o)
         Creature.stepRight(o)
         Creature.eatDown(o,Int(var_6))
         Creature.eatDown(o,Int(var_6))
         Creature.stepRight(o)
     end
     Creature.stepLeft(o)
     Creature.stepLeft(o)
     Creature.stepDown(o)
     for i = 1:var_1
         Creature.stepDown(o)
         Creature.stepLeft(o)
     end
     Creature.stepLeft(o)
     Creature.stepLeft(o)
     Creature.stepRight(o)
     return true
 end),[100,300,95,0,1,1,1],1,5000,1,11411174,0x01069818,Dict{Int16,Int16}(),[392,238],28),RpcApi.SimpleOrganism(0x000000000000c5f3,:(function (o::Creature.Organism,)
     local var_7::Int64 = -2881162876892958953
     local var_1::Int8 = 121
     Creature.eatUp(o,Int(var_1))
     for i = 1:var_1
         local var_8::Int16 = 23832
         local var_6::Int8 = 66
         Creature.eatRight(o,Int(var_6))
         Creature.eatRight(o,Int(var_6))
         Creature.eatRight(o,Int(var_6))
         Creature.eatRight(o,Int(var_6))
         Creature.eatRight(o,Int(var_6))
         Creature.eatRight(o,Int(var_6))
         Creature.eatDown(o,Int(var_6))
         Creature.eatRight(o,Int(var_6))
         Creature.eatRight(o,Int(var_6))
         Creature.eatUp(o,Int(var_6))
         var_8 = if haskey(o.mem,var_8)
                 o.mem[var_8]
             else
                 var_8
             end
         Creature.stepRight(o)
         Creature.stepRight(o)
     end
     Creature.eatUp(o,Int(var_1))
     Creature.stepLeft(o)
     Creature.stepLeft(o)
     for i = 1:var_1
         Creature.stepDown(o)
         Creature.stepLeft(o)
         Creature.stepRight(o)
     end
     Creature.stepUp(o)
     var_7 = var_7 / var_7
     Creature.stepLeft(o)
     Creature.stepLeft(o)
     Creature.stepUp(o)
     return true
 end),[100,300,95,0,1,1,1],1,5000,1,10993446,0x01070550,Dict{Int16,Int16}(),[440,333],31),RpcApi.SimpleOrganism(0x000000000000c5cc,:(function (o::Creature.Organism,)
     local var_7::Int64 = -2881162876892958953
     local var_1::Int8 = 121
     Creature.eatRight(o,Int(var_1))
     Creature.eatUp(o,Int(var_1))
     for i = 1:var_1
         local var_8::Int16 = 23832
         local var_6::Int8 = 66
         Creature.eatLeft(o,Int(var_6))
         Creature.eatRight(o,Int(var_6))
         Creature.eatRight(o,Int(var_6))
         Creature.eatLeft(o,Int(var_6))
         Creature.eatRight(o,Int(var_6))
         Creature.stepUp(o)
         Creature.eatLeft(o,Int(var_6))
         Creature.eatRight(o,Int(var_6))
         Creature.eatLeft(o,Int(var_6))
         Creature.stepRight(o)
     end
     Creature.stepLeft(o)
     Creature.stepLeft(o)
     for i = 1:var_1
         Creature.stepLeft(o)
         Creature.stepUp(o)
         Creature.stepLeft(o)
     end
     for i = 1:var_1
     end
     Creature.stepLeft(o)
     Creature.stepLeft(o)
     Creature.stepUp(o)
     return true
 end),[100,300,95,0,1,1,1],1,5000,1,10218775,0x01069818,Dict{Int16,Int16}(),[4,1],27),RpcApi.SimpleOrganism(0x000000000000c516,:(function (o::Creature.Organism,)
     local var_7::Int64 = -2881162876892958953
     local var_1::Int8 = 121
     Creature.eatUp(o,Int(var_1))
     for i = 1:var_1
         local var_8::Int16 = 23832
         local var_6::Int8 = 66
         Creature.eatRight(o,Int(var_6))
         Creature.eatRight(o,Int(var_6))
         Creature.eatRight(o,Int(var_6))
         Creature.eatRight(o,Int(var_6))
         Creature.eatRight(o,Int(var_6))
         Creature.eatRight(o,Int(var_6))
         Creature.eatDown(o,Int(var_6))
         Creature.stepRight(o)
         Creature.eatRight(o,Int(var_6))
         Creature.eatUp(o,Int(var_6))
         var_8 = if haskey(o.mem,var_8)
                 o.mem[var_8]
             else
                 var_8
             end
         Creature.stepRight(o)
         Creature.stepRight(o)
     end
     Creature.eatUp(o,Int(var_1))
     Creature.stepLeft(o)
     Creature.stepLeft(o)
     for i = 1:var_1
         Creature.stepDown(o)
         Creature.stepLeft(o)
         Creature.stepRight(o)
     end
     Creature.stepUp(o)
     var_7 = var_7 / var_7
     Creature.stepLeft(o)
     Creature.stepLeft(o)
     Creature.stepUp(o)
     return true
 end),[100,300,95,0,1,1,1],1,5000,1,9228585,0x01070550,Dict{Int16,Int16}(),[582,333],31),RpcApi.SimpleOrganism(0x000000000000c54e,:(function (o::Creature.Organism,)
     function func_10(var_9::Int64=-1465987312700749350)
         return var_9
     end
     local var_7::Int64 = -2881162876892958953
     local var_1::Int8 = 121
     Creature.eatUp(o,Int(var_1))
     for i = 1:var_1
         local var_8::Int16 = 23832
         local var_6::Int8 = 66
         Creature.eatRight(o,Int(var_6))
         Creature.eatRight(o,Int(var_6))
         Creature.eatRight(o,Int(var_6))
         Creature.eatRight(o,Int(var_6))
         Creature.eatRight(o,Int(var_6))
         Creature.eatRight(o,Int(var_6))
         Creature.eatDown(o,Int(var_6))
         Creature.stepRight(o)
         Creature.eatRight(o,Int(var_6))
         Creature.eatUp(o,Int(var_6))
         var_8 = if haskey(o.mem,var_8)
                 o.mem[var_8]
             else
                 var_8
             end
         Creature.stepRight(o)
         Creature.stepRight(o)
     end
     Creature.eatUp(o,Int(var_1))
     Creature.stepLeft(o)
     Creature.stepLeft(o)
     for i = 1:var_1
         Creature.stepDown(o)
         Creature.stepLeft(o)
         Creature.stepRight(o)
     end
     Creature.stepUp(o)
     var_7 = var_7 / var_7
     Creature.stepLeft(o)
     Creature.stepLeft(o)
     Creature.stepUp(o)
     return true
 end),[100,300,95,0,1,1,1],1,5000,1,8601752,0x01070550,Dict{Int16,Int16}(),[439,333],32),RpcApi.SimpleOrganism(0x000000000000c47e,:(function (o::Creature.Organism,)
     function func_10(var_9::Int64=-1465987312700749350)
         return var_9
     end
     local var_7::Int64 = -2881162876892958953
     local var_1::Int8 = 121
     Creature.eatUp(o,Int(var_1))
     for i = 1:var_1
         local var_8::Int16 = 23832
         local var_6::Int8 = 66
         Creature.eatRight(o,Int(var_6))
         Creature.eatRight(o,Int(var_6))
         Creature.eatRight(o,Int(var_6))
         Creature.eatRight(o,Int(var_6))
         Creature.eatRight(o,Int(var_6))
         Creature.eatRight(o,Int(var_6))
         Creature.eatDown(o,Int(var_6))
         Creature.stepRight(o)
         Creature.eatRight(o,Int(var_6))
         Creature.eatUp(o,Int(var_6))
         var_8 = if haskey(o.mem,var_8)
                 o.mem[var_8]
             else
                 var_8
             end
         Creature.stepRight(o)
         Creature.stepRight(o)
     end
     Creature.eatUp(o,Int(var_1))
     Creature.stepLeft(o)
     Creature.stepLeft(o)
     for i = 1:var_1
         Creature.stepDown(o)
         Creature.stepLeft(o)
         Creature.stepRight(o)
     end
     Creature.stepUp(o)
     var_7 = var_7 / var_7
     Creature.stepLeft(o)
     Creature.stepLeft(o)
     Creature.stepUp(o)
     return true
 end),[100,300,95,0,1,1,1],1,5000,1,8416280,0x01070550,Dict{Int16,Int16}(),[390,333],32)]

max_stars_repo_name: tmptrash/jevo

max_issues_repo_name: tmptrash/jevo

max_forks_repo_name: tmptrash/jevo


text:
prune(strategy, m) = mappruneable(strategy, m; exclude = x -> isleaf(x) && x isa AbstractArray)
prune(strategies::Union{<:Tuple, <:AbstractVector}, m::Chain) = Chain(prune.(strategies, m)...)

function iterativeprune(finetune!, stages, m; maxtries = 5)
    for stage in stages
        # prune to current stage
        m̄ = prune(stage, m)
        @show count(iszero, m̄[end].weight) / length(m̄[end].weight)

        # finetune pruned model as needed
        attempts = 0
        while (attempts < maxtries) && !finetune!(m̄)
            attempts += 1
        end
        m = m̄
    end

    return m
end


text:
module WaveletTypes
export  DiscreteWavelet,
        ContinuousWavelet,
        FilterWavelet,
        LSWavelet,
        WaveletType,
        #
        WaveletBoundary,
        PerBoundary,
        Periodic,
        #ZPBoundary,
        #CPBoundary,
        #
        OrthoFilter,
        #BiOrthoFilter,
        LSstep,
        GLS,
        #
        scale,
        wavelet,
        waveletfilter,
        waveletls
        #daubechies

using Compat

# TYPE HIERARCHY

abstract DiscreteWavelet{T}
abstract ContinuousWavelet{T}
# discrete transforms via filtering
abstract FilterWavelet{T} <: DiscreteWavelet{T}
# discrete transforms via lifting
abstract LSWavelet{T} <: DiscreteWavelet{T}
# all wavelet types
typealias WaveletType Union(DiscreteWavelet, ContinuousWavelet)

# BOUNDARY TYPES

abstract WaveletBoundary
# periodic (default)
immutable PerBoundary <: WaveletBoundary end
# zero padding
#immutable ZPBoundary <: WaveletBoundary end
# constant padding
#immutable CPBoundary <: WaveletBoundary end
# and so on...

const DEF_BOUNDARY = PerBoundary()

const Periodic = PerBoundary()

# ------------------------------------------------------------

# Compute filters from the Daubechies class
# N is the number of zeros at -1
function daubechies(N::Int)
    @assert N > 0
    # Create polynomial
    C = Array(Int, N)
    @inbounds for n = 0:N-1
        C[N-n] = binomial(N-1+n, n)
    end

    # Find roots in y domain (truncated binomial series; (1 - y)^{-N})
    Y = roots(C)

    # Find roots in z domain:
    # z + z^{-1} = 2 - 4*y
    # where y is a root from above
    Z = zeros(Complex128, 2*N-2)
    @inbounds for i = 1:N-1
        Yi = Y[i]
        d = 2*sqrt( Yi*Yi - Yi )
        y2 = 1 - 2*Yi
        Z[i] = y2 + d
        Z[i+N-1] = y2 -d
    end

    # Retain roots inside unit circle
    nr = 0  # count roots
    @inbounds for i = 1:length(Z)
        if abs(Z[i]) <= 1 + eps()
            nr += 1
        end
    end

    # Find coefficients of the polynomial
    # (1 + z)^N * \prod_i (z - z_i)
    R = Array(Complex128, N+nr)
    @inbounds for i = 1:N
        R[i] = -1
    end
    k = N
    @inbounds for i = 1:length(Z)
        if abs(Z[i]) <= 1 + eps()
            k += 1
            R[k] = Z[i]
        end
    end
    HH = vieta( R )

    # Normalize coefficients
    scale!(HH, 1/norm(HH))
    return real(HH)
end

# Compute roots of polynomial
# Input is a coefficient vector with highest powers first
function roots(C::AbstractVector)
    A = compan(C)
    return eigvals(A)
end

# Create companion matrix for a polynomial
# Input is a coefficient vector with highest powers first
function compan(C::AbstractVector)
    n = length(C)
    A = zeros(n-1, n-1)

    if n > 1
        @inbounds A[1,:] = -C[2:end] ./ C[1]
        @inbounds A[2:n:end] = 1
    end
    return A
end

# Vieta-like formula for computing polynomial coefficients from roots
# See
# http://www.mathworks.se/help/matlab/ref/poly.html
function vieta(R::AbstractVector)
    n = length( R )
    C = zeros(Complex128, n+1)
    C[1] = 1
    Ci::Complex128 = 0
    Cig::Complex128 = 0

    @inbounds for k = 1:n
        Ci = C[1]
        for i = 1:k
            Cig = C[i+1]
            C[i+1] = Cig - R[k] * Ci
            Ci = Cig
        end
    end
    return C
end

# IMPLEMENTATIONS OF FilterWavelet

immutable OrthoFilter{T<:WaveletBoundary} <: FilterWavelet{T}
    qmf::Vector{Float64}        # quadrature mirror filter
    name::ASCIIString           # filter short name
    OrthoFilter(qmf, name) = new(qmf, name)
end
function OrthoFilter{T<:WaveletBoundary}(name::String, ::T=DEF_BOUNDARY)
    if name[1:2] == "db"
        d = name[3:end]
        if isdigit(d)
            qmf = daubechies(int(d))
        else
            error("filter not found")
        end
    else
        qmf = get(FILTERS, name, nothing)
        qmf == nothing && error("filter not found")
    end
    # make sure it is normalized in l2-norm
    return OrthoFilter{T}(qmf./norm(qmf), name)
end
Base.length(f::OrthoFilter) = length(f.qmf)

#immutable BiOrthoFilter{T<:WaveletBoundary} <: FilterWavelet{T}
#    qmf1::Vector{Float64}       # quadrature mirror filter 1
#    qmf2::Vector{Float64}       # quadrature mirror filter 2
#    name::ASCIIString           # filter short name
#    BiOrthoFilter(qmf1, qmf2, name) = new(qmf1, qmf2, name)
#end

# IMPLEMENTATIONS OF LSWavelet

immutable LSstep
    stept::Char             # step type: 'p' for predict, 'u' for update
    coef::Vector{Float64}   # lifting coefficients
    shift::Int              # + left shift, - right shift
    function LSstep(stept,coef,shift)
        @assert stept=='p' || stept=='u'
        return new(stept,coef,shift)
    end
end
# general lifting scheme
immutable GLS{T<:WaveletBoundary} <: LSWavelet{T}
    step::Vector{LSstep}    # steps to be taken
    norm1::Float64          # normalization of scaling coefs.
    norm2::Float64          # normalization of detail coefs.
    name::ASCIIString       # name of scheme
    GLS(step, norm1, norm2, name) = new(step, norm1, norm2, name)
end
function GLS{T<:WaveletBoundary}(name::String, ::T=DEF_BOUNDARY)
    schemedef = get(SCHEMES, name, nothing)
    schemedef == nothing && error("scheme not found")
    return GLS{T}(schemedef..., name)
end

# IMPLEMENTATIONS OF ContinuousWavelet

# ...

# CLASSES

# convert class and number to name
function getname(class::String, n::Union(Integer,String))
    namebase = get(FILTERC2N, class, nothing)
    namebase == nothing && error("wavelet type not found")
    return string(namebase,n)
end

for f in (:OrthoFilter, :GLS)
@eval begin
	# convert class and number to wavelet name
	($f)(class::String, n::Union(Integer,String), args...) = ($f)(getname(class, n), args...)
end
end

# UTIL

# scale filter by scalar
function scale{T<:WaveletBoundary}(f::OrthoFilter{T}, a::Number)
    return OrthoFilter{T}(f.qmf.*a, f.name)
end

# TRANSFORM TYPE CONSTRUCTORS

function wavelet(name::String, boundary::WaveletBoundary=DEF_BOUNDARY)
    return waveletfilter(name, boundary)
end
function waveletfilter(name::String, boundary::WaveletBoundary=DEF_BOUNDARY)
    return OrthoFilter(name, boundary)
end
#function waveletfilterbo(name::String; boundary::WaveletBoundary=DEF_BOUNDARY)
#    return BiOrthoFilter(name, boundary)
#end
function waveletls(name::String, boundary::WaveletBoundary=DEF_BOUNDARY)
    return GLS(name, boundary)
end

for f in (:wavelet, :waveletfilter, :waveletls)
@eval begin
	# convert to type name
	($f)(class::String, n::Union(Integer,String), arg...) = ($f)(getname(class,n), arg...)
end
end

# class => namebase
const FILTERC2N=@compat Dict{ASCIIString,ASCIIString}(
"Haar" => "haar",
"Coiflet" => "coif",
"Daubechies" => "db",
"Symmlet" => "sym",
"Symlet" => "sym",
"Battle" => "batt",
"Beylkin" => "beyl",
"Vaidyanathan" => "vaid",
"CDF" => "cdf"
)

# scaling filters h (low pass)
# the number at end of a filter name is the
# number of vanishing moments of the mother wavelet function
# sources:
# http://statweb.stanford.edu/~wavelab/ (Orthogonal/MakeONFilter.m)
# http://www.mathworks.com/matlabcentral/fileexchange/5502-filter-coefficients-to-popular-wavelets
### https://github.com/nigma/pywt/blob/master/src/wavelets_coeffs.template.h
# name => qmf
const FILTERS=@compat Dict{ASCIIString,Vector{Float64}}(
# Haar filter, same as db1
"haar" =>
[0.7071067811865475,0.7071067811865475]
,
# Daubechies filters
"db1" =>
[0.7071067811865475,0.7071067811865475]
,
"db2" =>
[0.4829629131445341,0.8365163037378077,0.2241438680420134,-0.1294095225512603]
,
"db3" =>
[0.3326705529500827,0.8068915093110928,0.4598775021184915,-0.1350110200102546,-0.0854412738820267,0.0352262918857096]
,
"db4" =>
[0.2303778133074431,0.7148465705484058,0.6308807679358788,-0.0279837694166834,-0.1870348117179132,0.0308413818353661,0.0328830116666778,-0.0105974017850021]
,
"db5" =>
[0.1601023979741930,0.6038292697971898,0.7243085284377729,0.1384281459013204,-0.2422948870663824,-0.0322448695846381,0.0775714938400459,-0.0062414902127983,-0.0125807519990820,0.0033357252854738]
,
"db6" =>
[0.1115407433501094,0.4946238903984530,0.7511339080210954,0.3152503517091980,-0.2262646939654399,-0.1297668675672624,0.0975016055873224,0.0275228655303053,-0.0315820393174862,0.0005538422011614,0.0047772575109455,-0.0010773010853085]
,
"db7" =>
[0.0778520540850081,0.3965393194819136,0.7291320908462368,0.4697822874052154,-0.1439060039285293,-0.2240361849938538,0.0713092192668312,0.0806126091510820,-0.0380299369350125,-0.0165745416306664,0.0125509985560993,0.0004295779729214,-0.0018016407040474,0.0003537137999745]
,
"db8" =>
[0.0544158422431049,0.3128715909143031,0.6756307362972904,0.5853546836541907,-0.0158291052563816,-0.2840155429615702,0.0004724845739124,0.1287474266204837,-0.0173693010018083,-0.0440882539307952,0.0139810279173995,0.0087460940474061,-0.0048703529934518,-0.0003917403733770,0.0006754494064506,-0.0001174767841248]
,
"db9" =>
[0.0380779473638791,0.2438346746125939,0.6048231236901156,0.6572880780512955,0.1331973858249927,-0.2932737832791761,-0.0968407832229524,0.1485407493381306,0.0307256814793395,-0.0676328290613302,0.0002509471148340,0.0223616621236805,-0.0047232047577520,-0.0042815036824636,0.0018476468830564,0.0002303857635232,-0.0002519631889427,0.0000393473203163]
,
"db10" =>
[0.0266700579005546,0.1881768000776863,0.5272011889317202,0.6884590394536250,0.2811723436606485,-0.2498464243272283,-0.1959462743773399,0.1273693403357890,0.0930573646035802,-0.0713941471663697,-0.0294575368218480,0.0332126740593703,0.0036065535669880,-0.0107331754833036,0.0013953517470692,0.0019924052951930,-0.0006858566949566,-0.0001164668551285,0.0000935886703202,-0.0000132642028945]
,
# Coiflet filters
"coif2" =>
[-0.072732619513,0.337897662458,0.852572020212,0.384864846864,-0.072732619513,-0.015655728135]
,
"coif4" =>
[0.0163873364635998,-0.0414649367819558,-0.0673725547222826,0.3861100668229939,0.8127236354493977,0.4170051844236707,-0.0764885990786692,-0.0594344186467388,0.0236801719464464,0.0056114348194211,-0.0018232088707116,-0.0007205494453679]
,
"coif6" =>
[-0.0037935128644910,0.0077825964273254,0.0234526961418362,-0.0657719112818552,-0.0611233900026726,0.4051769024096150,0.7937772226256169,0.4284834763776168,-0.0717998216193117,-0.0823019271068856,0.0345550275730615,0.0158805448636158,-0.0090079761366615,-0.0025745176887502,0.0011175187708906,0.0004662169601129,-0.0000709833031381,-0.0000345997728362]
,
"coif8" =>
[0.0008923136685824,-0.0016294920126020,-0.0073461663276432,0.0160689439647787,0.0266823001560570,-0.0812666996808907,-0.0560773133167630,0.4153084070304910,0.7822389309206135,0.4343860564915321,-0.0666274742634348,-0.0962204420340021,0.0393344271233433,0.0250822618448678,-0.0152117315279485,-0.0056582866866115,0.0037514361572790,0.0012665619292991,-0.0005890207562444,-0.0002599745524878,0.0000623390344610,0.0000312298758654,-0.0000032596802369,-0.0000017849850031]
,
"coif10" =>
[-0.0002120808398259,0.0003585896879330,0.0021782363583355,-0.0041593587818186,-0.0101311175209033,0.0234081567882734,0.0281680289738655,-0.0919200105692549,-0.0520431631816557,0.4215662067346898,0.7742896037334738,0.4379916262173834,-0.0620359639693546,-0.1055742087143175,0.0412892087544753,0.0326835742705106,-0.0197617789446276,-0.0091642311634348,0.0067641854487565,0.0024333732129107,-0.0016628637021860,-0.0006381313431115,0.0003022595818445,0.0001405411497166,-0.0000413404322768,-0.0000213150268122,0.0000037346551755,0.0000020637618516,-0.0000001674428858,-0.0000000951765727]
,
# Symmlet filter
"sym4" =>
[0.0455703458960000,-0.0178247014420000,-0.1403176241790000,0.4212345342040000,1.1366582434079999,0.7037390686560000,-0.0419109651250000,-0.1071489014180000]
,
"sym5" =>
[0.0276321529580000,-0.0298424998690000,-0.2479513626130000,0.0234789231360000,0.8965816483800000,1.0230529668940000,0.2819906968540000,-0.0553441861170000,0.0417468644220000,0.0386547959550000]
,
"sym6" =>
[-0.0110318675090000,0.0024999220930000,0.0632505626600000,-0.0297837512990000,-0.1027249698620000,0.4779043713330000,1.1138927839260000,0.6944579729580000,-0.0683231215870000,-0.1668632154120000,0.0049366123720000,0.0217847003270000]
,
"sym7" =>
[0.0145213947620000,0.0056713426860000,-0.1524638718960000,-0.1980567068070000,0.4081839397250000,1.0857827098140000,0.7581626019640000,0.0246656594890000,-0.0700782912220000,0.0960147679360000,0.0431554525820000,-0.0178704316510000,-0.0014812259150000,0.0037926585340000]
,
"sym8" =>
[-0.0047834585120000,-0.0007666908960000,0.0448236230420000,0.0107586117510000,-0.2026486552860000,-0.0866536154060000,0.6807453471900000,1.0991066305370001,0.5153986703740000,-0.0734625087610000,-0.0384935212630000,0.0694904659110000,0.0053863887540000,-0.0211456865280000,-0.0004283943000000,0.0026727933930000]
,
"sym9" =>
[0.0019811937360000,0.0008765025390000,-0.0187693968360000,-0.0163033512260000,0.0427444336020000,0.0008251409290000,-0.0771721610970000,0.3376589236020000,1.0152597908320000,0.8730484073490000,0.0498828309590000,-0.2708937835030000,-0.0257864459300000,0.0877912515540000,0.0125288962420000,-0.0145155785530000,-0.0006691415090000,0.0015124873090000]
,
"sym10" =>
[-0.0006495898960000,0.0000806612040000,0.0064957283750000,-0.0011375353140000,-0.0287862319260000,0.0081528167990000,0.0707035675500000,-0.0452407722180000,-0.0502565400920000,0.5428130112130000,1.0882515305000000,0.6670713381540000,-0.1002402150310000,-0.2255589722340000,0.0164188694260000,0.0649509245790000,-0.0020723639230000,-0.0122206426300000,0.0001352450200000,0.0010891704470000]
,
# Battle-Lemarie filter
"batt2" =>
[-0.0000867523000000,-0.0001586010000000,0.0003617810000000,0.0006529220000000,-0.0015570100000000,-0.0027458800000000,0.0070644200000000,0.0120030000000000,-0.0367309000000000,-0.0488618000000000,0.2809310000000000,0.5781630000000000,0.2809310000000000,-0.0488618000000000,-0.0367309000000000,0.0120030000000000,0.0070644200000000,-0.0027458800000000,-0.0015570100000000,0.0006529220000000,0.0003617810000000,-0.0001586010000000,-0.0000867523000000]
,
"batt4" =>
[0.0001033070000000,-0.0001642640000000,-0.0002018180000000,0.0003267490000000,0.0003959460000000,-0.0006556200000000,-0.0007804680000000,0.0013308600000000,0.0015462400000000,-0.0027452900000000,-0.0030786300000000,0.0057993200000000,0.0061414300000000,-0.0127154000000000,-0.0121455000000000,0.0297468000000000,0.0226846000000000,-0.0778079000000000,-0.0354980000000000,0.3068300000000000,0.5417360000000000,0.3068300000000000,-0.0354980000000000,-0.0778079000000000,0.0226846000000000,0.0297468000000000,-0.0121455000000000,-0.0127154000000000,0.0061414300000000,0.0057993200000000,-0.0030786300000000,-0.0027452900000000,0.0015462400000000,0.0013308600000000,-0.0007804680000000,-0.0006556200000000,0.0003959460000000,0.0003267490000000,-0.0002018180000000,-0.0001642640000000,0.0001033070000000]
,
"batt6" =>
[0.0001011130000000,0.0001107090000000,-0.0001591680000000,-0.0001726850000000,0.0002514190000000,0.0002698420000000,-0.0003987590000000,-0.0004224850000000,0.0006355630000000,0.0006628360000000,-0.0010191200000000,-0.0010420700000000,0.0016465900000000,0.0016413200000000,-0.0026864600000000,-0.0025881600000000,0.0044400200000000,0.0040788200000000,-0.0074684800000000,-0.0063988600000000,0.0128754000000000,0.0099063500000000,-0.0229951000000000,-0.0148537000000000,0.0433544000000000,0.0208414000000000,-0.0914068000000000,-0.0261771000000000,0.3128690000000000,0.5283740000000000,0.3128690000000000,-0.0261771000000000,-0.0914068000000000,0.0208414000000000,0.0433544000000000,-0.0148537000000000,-0.0229951000000000,0.0099063500000000,0.0128754000000000,-0.0063988600000000,-0.0074684800000000,0.0040788200000000,0.0044400200000000,-0.0025881600000000,-0.0026864600000000,0.0016413200000000,0.0016465900000000,-0.0010420700000000,-0.0010191200000000,0.0006628360000000,0.0006355630000000,-0.0004224850000000,-0.0003987590000000,0.0002698420000000,0.0002514190000000,-0.0001726850000000,-0.0001591680000000,0.0001107090000000,0.0001011130000000]
,
# Beylkin filter
"beyl" =>
[0.0993057653740000,0.4242153608130000,0.6998252140570000,0.4497182511490000,-0.1109275983480000,-0.2644972314460000,0.0269003088040000,0.1555387318770000,-0.0175207462670000,-0.0885436306230000,0.0196798660440000,0.0429163872740000,-0.0174604086960000,-0.0143658079690000,0.0100404118450000,0.0014842347820000,-0.0027360316260000,0.0006404853290000]
,
# Vaidyanathan filter
"vaid" =>
[-0.0000629061180000,0.0003436319050000,-0.0004539566200000,-0.0009448971360000,0.0028438345470000,0.0007081375040000,-0.0088391034090000,0.0031538470560000,0.0196872150100000,-0.0148534480050000,-0.0354703986070000,0.0387426192930000,0.0558925236910000,-0.0777097509020000,-0.0839288843660000,0.1319716614170000,0.1350842271290000,-0.1944504717660000,-0.2634948024880000,0.2016121617750000,0.6356010598720000,0.5727977932110000,0.2501841295050000,0.0457993341110000]

)

# biortho filters
# name => (qmf1,qmf2)
const BIFILTERS=@compat Dict{ASCIIString,NTuple{2,Vector{Float64}}}(
# test
"test" =>
([0.7071067811865475,0.7071067811865475],
 [0.7071067811865475,0.7071067811865475])

)

# in matlab (d,p) -> (p,u)
const SCHEMES=@compat Dict{ASCIIString,NTuple{3}}(
# cdf 5/3 -> bior 2.2, cdf 9/7 -> bior 4.4
# Cohen-Daubechies-Feauveau [Do Quan & Yo-Sung Ho. Optimized median lifting scheme for lossy image compression.]
"cdf9/7" => ([ LSstep('u',1.5861343420604  *[1.0,1.0],0),
			LSstep('p',0.05298011857291494 *[1.0,1.0],1),
			LSstep('u',-0.882911075531393  *[1.0,1.0],0),
            LSstep('p',-0.44350685204384654*[1.0,1.0],1)],
            1.1496043988603355,
            0.8698644516247099)
,

# Haar, same as db1
"haar" => ([ LSstep('p',[-1.0],0),
            LSstep('u',[0.5],0)],
            0.7071067811865475,
            1.4142135623730951)
,
# Daubechies
"db1" => ([ LSstep('p',[-1.0],0),
            LSstep('u',[0.5],0)],
            0.7071067811865475,
            1.4142135623730951)
,
"db2" => ([ LSstep('p',[-1.7320508075688772],0),
            LSstep('u',[-0.0669872981077807,0.4330127018922193],1),
            LSstep('p',[1.0],-1)],
            0.5176380902050414,
            1.9318516525781364)

)

end


text:
function tss(y::AbstractVector, hasintercept::Bool, weights::AbstractWeights)
    if hasintercept
        m = mean(y, weights)
        sum(@inbounds (y[i] - m)^2 * weights[i] for i in eachindex(y))
    else
        sum(@inbounds y[i]^2 * weights[i] for i in eachindex(y))
    end
end


text:
# Left-looking Cholesky
@inline function ExaTronKernels.dicf(n::Int,L::ROCDeviceArray{Float64,2})
    tx = workitemIdx().x

    @inbounds for j=1:n
        # Apply the pending updates.
        if tx >= j && tx <= n
            for k=1:j-1
                L[tx,j] -= L[tx,k] * L[j,k]
            end
        end
        AMDGPU.sync_workgroup()

        if (L[j,j] <= 0)
            AMDGPU.sync_workgroup()
            return -1
        end

        Ljj = sqrt(L[j,j])
        if tx >= j && tx <= n
            L[tx,j] /= Ljj
        end
        AMDGPU.sync_workgroup()
    end

    if tx <= n
        @inbounds for j=1:n
            if tx > j
                L[j,tx] = L[tx,j]
            end
        end
    end
    AMDGPU.sync_workgroup()

    return 0
end


text:
using HeatEqCalc
sol = HeatEqCalc.calc_sol()
HeatEqCalc.save_sol(sol)

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public