text:
module SurveyDataWeighting

using Base: Tuple
using Base.Threads

export do_chi_square_reweighting, do_reweighting
export DistanceFunctionType, chi_square, d_and_s_type_a, d_and_s_type_b, constrained_chi_square, d_and_s_constrained
export ITERATIONS_EXCEEDED

"""
   Implements the micro data weighting procedures from:
   * Creedy 2003 http://www.treasury.govt.nz/publications/research-policy/wp/2003/03-17/twp03-17.pdf
   * Creedy 2003  http://www.business.curtin.edu.au/files/creedy2.pdf
   * Jean-Claude Deville and Carl-Erik Sarndal http://www.jstor.org/stable/2290268
"""

"""
Possible distance types; see Creedy.
"""
@enum DistanceFunctionType chi_square d_and_s_type_a d_and_s_type_b constrained_chi_square d_and_s_constrained

const ITERATIONS_EXCEEDED = -1

""""
Make a weights vector which weights the matrix `data`
so when summed the col totals sum to `target_populations`
See the Creedy Paper for `function_type`
If using one of the constrained types,
the output weights should be no more than ru*the initial weight,
no less than rl
Returns a Dict with :=>weights and some extra info on convergence.
data : KxJ matrix where k is num observations and J is num constraints;
see:
Microdata Adjustment by the Minimum Information Loss Principle Joachim Merz; FFB Discussion Paper No. 10 July 1994
for a good discussion on how to lay out the dataset

intial_weights, new_weights : K length vector
target_populations - J length vector;

tolx, tolf, max_iterations : see Solve_Non_Linear_Equation_System in the parent
ru/rl max/min acceptable values of ratio of final_weight/initial_weight (for constrained distance functions)

note: chi-square is just there for checking purposes; use `do_chi_square_reweighting` if that's all you need.

"""
function do_reweighting(
    ;
    data               :: AbstractArray{ <:Real, 2 },
    initial_weights    :: AbstractArray{ <:Real, 1 }, # a column
    target_populations :: AbstractArray{ <:Real, 1 }, # a row
    functiontype       :: DistanceFunctionType,
    upper_multiple     :: Real = 0.0,
    lower_multiple     :: Real = 0.0,
    tolx               :: Real = 0.000001,
    tolf               :: Real = 0.000001 ) :: NamedTuple

    nrows = size( data )[1]
    ncols = size( data )[2]
    @assert ncols == size( target_populations )[1]
    @assert nrows == size( initial_weights )[1]
    a = target_populations - (initial_weights'*data)'
    lamdas = zeros( Float64, ncols )
    ru = upper_multiple # shorthand
    rl = lower_multiple
    ##
    ## horror cpde: document problem with calling func with
    ## variable input in J1.0. Hence this version
    ## where lamdas is a semi-global so we don't modify it
    ## before we call `thefunc`.
    ##
    function local_solve_non_linear_equation_system(
        thefunc,
        numtrials :: Integer = 50 ) :: Dict{Symbol,Any}

        deltas = zeros( ncols )
        error = 0
        iterations = 0
        for k in 1:numtrials
            iterations += 1
            errf = 0.0
            errx = 0.0
            outx :: Dict{Symbol, Any} = thefunc()
            gradient = outx[:gradient]
            hessian = outx[:hessian]
            for i in 1:ncols
                 errf += abs( gradient[i])
            end
            if errf <= tolf
                 break
            end
            deltas = hessian \ gradient
            lamdas += deltas
            for i in 1:ncols
                 errx += abs( deltas[i])
            end
            if errx <= tolx
                 break
            end
        end # trials
        if iterations == numtrials
             error = ITERATIONS_EXCEEDED
        end
        return Dict( :iterations=>iterations, :error=>error )
    end

    function compute_lamdas_and_hessian()  :: Dict{ Symbol, Any }
        gradient = zeros( Float64, ncols, 1 )
        hessian = zeros( Float64, ncols, ncols )
        z = zeros( Float64, ncols, 1 )
        for row in 1:nrows
            rv = data[row,:]
            u = (rv' * lamdas)[1]
            # println("u=$u lamdas=$lamdas")
            d_g_m1 = 0.0
            g_m1 = 0.0
            if functiontype == chi_square
                d_g_m1 = 1.0;
                g_m1 = 1.0 + u;
            elseif functiontype == constrained_chi_square
                if( u < ( rl - 1.0 ))
                   g_m1 = rl
                   d_g_m1 = 0.0
                elseif( u > ( ru - 1.0 ))
                   g_m1 = ru
                   d_g_m1 = 0.0
                else
                   g_m1 = 1.0 + u
                   d_g_m1 = 1.0
                end
            elseif functiontype == d_and_s_type_a
               g_m1 = ( 1.0 -  u/2.0 ) ^ ( -2 )
               d_g_m1 = ( 1.0 - u/2.0 ) ^ ( -3 )
            elseif functiontype == d_and_s_type_b
               g_m1 = ( 1.0- u ) ^ (-1 )
               d_g_m1 = ( 1.0 - u ) ^ ( -2 )
            elseif functiontype == d_and_s_constrained
               alpha = ( ru - rl ) / (( 1.0 - rl )*( ru - 1.0 ))
               g_m1 = rl*(ru-1.0)+ru*(1.0-rl)*exp( alpha*u )/((ru-1.0)+(1.0-rl)*(exp( alpha*u )))
               d_g_m1 = g_m1 * ( ru - g_m1 ) *
                 ((( 1.0 - rl )*alpha*exp( alpha*u )) /
                  (( ru - 1.0 ) + (( 1.0 - rl ) * exp( alpha*u ))))
           end # function cases
           for col in 1:ncols
               z[col] += initial_weights[row]*data[row,col]*(g_m1-1.0)
               ## the hessian
               for c2 in 1:ncols
                   zz :: Float64 = initial_weights[row]*data[row,col]*data[row,c2]
                   hessian[col,c2] += zz*d_g_m1
               end
           end
        end # obs loop
        gradient = a - z
        # println( "gradient $gradient")
        # println( "hessian $hessian")
        d = Dict(:gradient=>gradient,:hessian=>hessian )
        return d
    end # nested function

    function make_start_stops( nrows::Int, num_threads::Int )::Tuple
        start = zeros(Int, num_threads)
        stop = zeros(Int, num_threads)
        chunk_size = Int(trunc(nrows/num_threads))-1;
        p = 1;
        for i in 1:num_threads
            start[i] = p
            p = p + chunk_size;
            if i < num_threads;
                stop[i] = p;
            else;
                stop[i] = nrows; # possibly different number on last thread
            end;
            p = p + 1;
        end;
        return (start, stop)
    end

    function compute_lamdas_and_hessian_threaded()  :: Dict{ Symbol, Any }
        num_threads = nthreads()
        start,stop = make_start_stops( nrows, num_threads )
        t_hessian = Array{Matrix{Float64}}(undef, num_threads )
        t_z = Array{Matrix{Float64}}(undef, num_threads ) # FIXME really a vector

        @threads for thread in 1:num_threads
            # println( "starting thread $thread with obs $(start[thread]) - $(stop[thread])")
            gradient = zeros( Float64, ncols, 1 )
            hessian = zeros( Float64, ncols, ncols )
            z = zeros( Float64, ncols, 1 )
            for row in start[thread]:stop[thread]
                rv = data[row,:]
                u = (rv' * lamdas)[1]
                # println("u=$u lamdas=$lamdas")
                d_g_m1 = 0.0
                g_m1 = 0.0
                if functiontype == chi_square
                    d_g_m1 = 1.0;
                    g_m1 = 1.0 + u;
                elseif functiontype == constrained_chi_square
                    if( u < ( rl - 1.0 ))
                        g_m1 = rl
                        d_g_m1 = 0.0
                    elseif( u > ( ru - 1.0 ))
                        g_m1 = ru
                        d_g_m1 = 0.0
                    else
                        g_m1 = 1.0 + u
                        d_g_m1 = 1.0
                    end
                elseif functiontype == d_and_s_type_a
                        g_m1 = ( 1.0 -  u/2.0 ) ^ ( -2 )
                        d_g_m1 = ( 1.0 - u/2.0 ) ^ ( -3 )
                    elseif functiontype == d_and_s_type_b
                        g_m1 = ( 1.0- u ) ^ (-1 )
                        d_g_m1 = ( 1.0 - u ) ^ ( -2 )
                    elseif functiontype == d_and_s_constrained
                        alpha = ( ru - rl ) / (( 1.0 - rl )*( ru - 1.0 ))
                        g_m1 = rl*(ru-1.0)+ru*(1.0-rl)*exp( alpha*u )/((ru-1.0)+(1.0-rl)*(exp( alpha*u )))
                        d_g_m1 = g_m1 * ( ru - g_m1 ) *
                            ((( 1.0 - rl )*alpha*exp( alpha*u )) /
                            (( ru - 1.0 ) + (( 1.0 - rl ) * exp( alpha*u ))))
                    end # function cases
                for col in 1:ncols
                    z[col] += initial_weights[row]*data[row,col]*(g_m1-1.0)
                    ## the hessian
                    for c2 in 1:ncols
                        zz :: Float64 = initial_weights[row]*data[row,col]*data[row,c2]
                        hessian[col,c2] += zz*d_g_m1
                    end
                end # ncols
            end # obs loop
            t_z[thread] = z
            t_hessian[thread] = hessian
        end # threads

        g_hessian = zeros( Float64, ncols, ncols )
        g_z = zeros( Float64, ncols, 1 )
        for i in 1:num_threads
            g_hessian += t_hessian[i]
            g_z += t_z[i]
        end
        gradient = a - g_z
        # println( "gradient $gradient")
        # println( "hessian $g_hessian")
        # println( "z $g_z")

        d = Dict(:gradient=>gradient,:hessian=>g_hessian )
        return d
    end # nested function

    rc = local_solve_non_linear_equation_system(
        compute_lamdas_and_hessian_threaded )

    new_weights = copy(initial_weights)
    @assert rc[:error] == 0
    for r in 1:nrows
        row = data[r,:]
        u = (row'*lamdas)[1]
        g_m1 = 0.0
        if functiontype == chi_square
            g_m1 = 1.0 + u;
        elseif functiontype == constrained_chi_square
            if( u < ( rl - 1.0 ))
               g_m1 = rl
            elseif( u > ( ru - 1.0 ))
               g_m1 = ru
            else
               g_m1 = 1.0 + u
            end
        elseif functiontype == d_and_s_type_a
           g_m1 = ( 1.0 -  u/2.0 ) ^ ( -2 )
        elseif functiontype == d_and_s_type_b
           g_m1 = ( 1.0- u ) ^ (-1 )
        elseif functiontype == d_and_s_constrained
           alpha = ( ru - rl ) / (( 1.0 - rl )*( ru - 1.0 ))
           g_m1 = rl*(ru-1.0)+ru*(1.0-rl)*exp( alpha*u )/((ru-1.0)+(1.0-rl)*(exp( alpha*u )))
       end # function cases
        #
        # Creedy wp 03/17 table 3
        #
        new_weights[r] = initial_weights[r]*g_m1
    end
    if functiontype in [constrained_chi_square, d_and_s_constrained ]
     # check the constrainted methods keep things inside ll and ul
        for r in 1:nrows
            @assert new_weights[r] <= initial_weights[r]*upper_multiple
            @assert new_weights[r] >= initial_weights[r]*lower_multiple
        end
    end #
    weighted_totals = (new_weights' * data)'
    @assert weighted_totals ≈ target_populations

    return (lamdas=lamdas, rc = rc, weights=new_weights ) # , :converged => converge )
end # do reweighting

"
This is a route-1 approach to Chi-square reweighting.
The iterative main method should produce identical results
when method=chi_square. This is kept here mainly for testing.
Note the weights can be negative.
See the Creedy Papers.
"
function do_chi_square_reweighting(
    data               :: AbstractArray{ <:Real, 2 },
    initial_weights    :: AbstractArray{ <:Real, 1 }, # a row
    target_populations :: AbstractArray{ <:Real, 1 } ) :: Array{ <:Real }

    nrows = size( data )[1]
    ncols = size( data )[2]

    row = zeros( ncols )
    populations = zeros( ncols, 1 )
    lamdas  = zeros( ncols, 1 )
    weights = zeros( nrows, 1 )
    m = zeros( ncols, ncols )
    for r in 1:nrows
        row = data[r,:]
        m += initial_weights[r]*(row*row')
        for c in 1:ncols
            populations[c] += (row[c]*initial_weights[r])'
        end
    end
    lamdas = (m^-1)*(target_populations-populations)
    for r in 1:nrows
        row = data[r,:]
        weights[r] = initial_weights[r]*(1.0 + (row'*lamdas)[1])
    end
    return weights;
end

end # package

max_stars_repo_name: grahamstark/SurveyDataWeighting

max_issues_repo_name: grahamstark/SurveyDataWeighting

max_forks_repo_name: grahamstark/SurveyDataWeighting


text:
module DispatchTests
using  ExpressionPatterns.Dispatch
using  Base.Test
# works?

@macromethod  w1(x) x
@metafunction w1(x) x
@macromethod  w2(x) = x
@metafunction w2(x) = x

@test @w1(1) == @w2(1) == w1(1) == w2(1)

# parameter dispatch

@macromethod param(x)   1
@macromethod param(x,y) 2

@test @param(a)   == 1
@test @param(a,b) == 2

# expression dispatch

@macromethod disp( x+y ) 1
@macromethod disp( x-y ) 2
@macromethod disp((x,y)) 3
@macromethod disp([x+y]) 4

@test @disp( x+y ) == 1
@test @disp( x-y ) == 2
@test @disp((x,y)) == 3
@test @disp([x+y]) == 4

# inteligent dispatch: matches the most specific, regardless of definition order

@macromethod intdisp1(x)   1
@macromethod intdisp1(x+y) 2

@macromethod intdisp2(x+y) 2
@macromethod intdisp2(x)   1

@test @intdisp1(x)   == 1
@test @intdisp1(x+y) == 2

@test @intdisp2(x)   == 1
@test @intdisp2(x+y) == 2

@metafunction intdisp3(:T{x, Number})  Number
@metafunction intdisp3(:T{x, Integer}) Integer

@metafunction intdisp4(:T{x, Integer}) Integer
@metafunction intdisp4(:T{x, Number})  Number

@test intdisp3(10.0) == Number
@test intdisp3(1000) == Integer

@test intdisp4(10.0) == Number
@test intdisp4(1000) == Integer

# dispatch on slurps

@macromethod sdisp1(x)     '1'
@macromethod sdisp1(x,y)   '2'
@macromethod sdisp1(*{xs}) 'n'

@test @sdisp1(x)     == '1'
@test @sdisp1(x,y)   == '2'
@test @sdisp1(x,y,z) == 'n'

@macromethod sdisp2(*{xs}) 'n'
@macromethod sdisp2(x)     '1'

@test @sdisp2(x)   == '1'
@test @sdisp2(x,y) == 'n'

@macromethod sdisp3(f(*{as}))        (f,as...)
@macromethod sdisp3(f(*{as}; *{ks})) (f,as...,ks...)

@test @sdisp3(f(x,y))     == (:f, :x, :y)
@test @sdisp3(f(x,y;k=1)) == (:f, :x, :y, Expr(:kw, :k, 1))

# generic dispatch

@metadispatch metad_mfs() = 0
@metadispatch function metad_mfl() 0 end
@metadispatch macro metad_mm() 0 end

@metadestruct metad_fs() = 0
@metadestruct function metad_fl() 0 end

@metadestruct let (x,) = :(0,)
  @metadestruct let a = 0, b = 0
    @test 0 ==
      (@metadestruct ()->0)()    ==
      metad_fs()  == metad_fl()  ==
      metad_mfs() == metad_mfl() ==
      (@metad_mm) == x == a == b
  end
end

# metamethod overwriting

@metafunction o()[x] 1
@metafunction o()[x] 2

@test o() == 2

# errors

@test_throws MetaMethodError o(1)
@test_throws ErrorException try @eval @metafunction o(x)[x] 0; catch err throw(err.error) end

end

ext: jl


text:
struct Poisson <: Distribution{Int} end

"""
    poisson(lambda::Real)

Sample an `Int` from the Poisson distribution with rate `lambda`.
"""
const poisson = Poisson()

function logpdf(::Poisson, x::Int, lambda::Real)
    x < 0 ? -Inf : x * log(lambda) - lambda - loggamma(x+1)
end

function logpdf_grad(::Poisson, x::Int, lambda::Real)
    (nothing, x/lambda - 1)
end

function random(::Poisson, lambda::Real)
    rand(Distributions.Poisson(lambda))
end

(::Poisson)(lambda) = random(Poisson(), lambda)
is_discrete(::Poisson) = true

has_output_grad(::Poisson) = false
has_argument_grads(::Poisson) = (false,)

export poisson


text:
#ДАНО: Робот - в произвольной клетке ограниченного прямоугольного поля
#РЕЗУЛЬТАТ: Робот - в исходном положении, и все клетки поля промакированы
function zacras(r::Robot)
    v = 0
    q = 0
    while isborder(r,HorizonSide(0)) == false
        move!(r,HorizonSide(0))
        v = v+1
    end
    while isborder(r,HorizonSide(1)) == false
        move!(r,HorizonSide(1))
        q = q+1
    end
    while isborder(r,HorizonSide(2)) == false
        if isborder(r,HorizonSide(1)) == true
            idyvpravo(r)
        end
        if isborder(r,HorizonSide(3)) == true
            idyvlevo(r)
        end
    end
    if isborder(r,HorizonSide(1))
        idyvpravo(r)
    else
        idyvlevo(r)
    end
    ugol(r)
    for i in 1:v
        move!(r, HorizonSide(2))
    end
    for i in 1:q
        move!(r, HorizonSide(3))
    end
end
function idyvpravo(r::Robot)
    putmarker!(r)
    while isborder(r,HorizonSide(3))==false
        move!(r,HorizonSide(3))
        putmarker!(r)
    end
    if isborder(r,HorizonSide(2)) == false
        move!(r, HorizonSide(2))
    end
end
function idyvlevo(r::Robot)
    putmarker!(r)
    while isborder(r,HorizonSide(1) )==false
        move!(r,HorizonSide(1))
        putmarker!(r)
    end
    if isborder(r,HorizonSide(2)) == false
        move!(r, HorizonSide(2))
    end
end
function ugol(r)
    while isborder(r, HorizonSide(0)) == false
        move!(r,HorizonSide(0))
    end
    while isborder(r, HorizonSide(1)) == false
        move!(r,HorizonSide(1))
    end
end

max_stars_repo_name: listova-nastya/mirea-progs

max_issues_repo_name: listova-nastya/mirea-progs

max_forks_repo_name: listova-nastya/mirea-progs


text:
function print_header_1()
    println("---------------------------------------------------------------------------------------")
    println("=======================================================================================")
    println("                  ProxSDP : Proximal Semidefinite Programming Solver                   ")
    println("                         (c) Mario Souto and Joaquim D. Garcia, 2020                   ")
    println("                                                              v1.6.0                   ")
    println("---------------------------------------------------------------------------------------")
end

function print_parameters(opt::Options, conic_sets::ConicSets)
    println("    Solver parameters:")
    tol_str = "       tol_gap = $(opt.tol_gap) tol_feasibility = $(opt.tol_feasibility)\n"
    tol_str *= "       tol_primal = $(opt.tol_primal) tol_dual = $(opt.tol_dual)"
    if length(conic_sets.socone) >= 1
        tol_str *= " tol_soc = $(opt.tol_soc)"
    end
    if length(conic_sets.sdpcone) >= 1
        tol_str *= " tol_psd = $(opt.tol_psd)"
    end
    println(tol_str)

    println("       max_iter = $(opt.max_iter_local) time_limit = $(opt.time_limit)s")

    return nothing
end

eq_plural(val::Integer) = ifelse(val != 1, "ies", "y")
eqs(val::Integer) = val > 0 ? "$(val) linear equalit$(eq_plural(val))" : ""
ineqs(val::Integer) = val > 0 ? "$(val) linear inequalit$(eq_plural(val))" : ""
function print_constraints(aff::AffineSets)
    println("    Constraints:")
    println("       $(eqs(aff.p)) and $(eqs(aff.m))")
    return nothing
end

cone_plural(val::Integer) = ifelse(val != 1, "s", "")
function print_prob_data(conic_sets::ConicSets)
    soc_dict = Dict()
    for soc in conic_sets.socone
        if soc.len in keys(soc_dict)
            soc_dict[soc.len] += 1
        else
            soc_dict[soc.len] = 1
        end
    end
    psd_dict = Dict()
    for psd in conic_sets.sdpcone
        if psd.sq_side in keys(psd_dict)
            psd_dict[psd.sq_side] += 1
        else
            psd_dict[psd.sq_side] = 1
        end
    end
    println("    Cones:")
    if length(conic_sets.socone) > 0
        for (k, v) in soc_dict
            println("       $v second order cone$(cone_plural(v)) of size $k")
        end
    end
    if length(conic_sets.sdpcone) > 0
        for (k, v) in psd_dict
            println("       $v psd cone$(cone_plural(v)) of size $k")
        end
    end

    return nothing
end

function print_header_2(opt, beg = true)
    bar  = "---------------------------------------------------------------------------------------"
    name = "    Initializing Primal-Dual Hybrid Gradient method"
    cols = "|  iter  | prim obj | rel. gap |  feasb.  | prim res | dual res | tg. rank |  time(s) |"

    if opt.extended_log || opt.extended_log2
        bar  *= "-----------"
        cols *= " dual obj |"
    end
    if opt.extended_log2
        bar  *= "-----------"
        cols *= " d feasb. |"
    end

    if beg
        println(bar)
        println(name)
        println(bar)
    end
    println(cols)
    if beg
        println(bar)
    end

    return nothing
end

function print_progress(residuals::Residuals, p::Params, opt, val = -1.0)
    primal_res = residuals.primal_residual[p.iter]
    dual_res = residuals.dual_residual[p.iter]
    s_k = @sprintf("%d", p.iter)
    s_k *= " |"
    s_s = @sprintf("%.5f", residuals.dual_gap[p.iter])
    s_s *= " |"
    s_o = @sprintf("%.3f", residuals.prim_obj[p.iter])
    s_o *= " |"
    s_f = @sprintf("%.5f", residuals.feasibility[p.iter])
    s_f *= " |"
    s_p = @sprintf("%.5f", primal_res)
    s_p *= " |"
    s_d = @sprintf("%.5f", dual_res)
    s_d *= " |"
    s_target_rank = @sprintf("%.0f", sum(p.target_rank))
    s_target_rank *= " |"
    s_time = @sprintf("%.4f", time() - p.time0)
    s_time *= " |"
    a = "|"
    a *= " "^max(0, 9 - length(s_k))
    a *= s_k
    a *= " "^max(0, 11 - length(s_o))
    a *= s_o
    a *= " "^max(0, 11 - length(s_s))
    a *= s_s
    a *= " "^max(0, 11 - length(s_f))
    a *= s_f
    a *= " "^max(0, 11 - length(s_p))
    a *= s_p
    a *= " "^max(0, 11 - length(s_d))
    a *= s_d
    a *= " "^max(0, 11 - length(s_target_rank))
    a *= s_target_rank
    a *= " "^max(0, 11 - length(s_time))
    a *= s_time

    if opt.extended_log || opt.extended_log2
        str = @sprintf("%.3f", residuals.dual_obj[p.iter]) * " |"
        str = " "^max(0, 11 - length(str)) * str
        a *= str
    end
    if opt.extended_log2
        str = @sprintf("%.5f", val) * " |"
        str = " "^max(0, 11 - length(str)) * str
        a *= str
    end
    if opt.log_repeat_header
        print_header_2(opt, false)
    end

    println(a)

    return nothing
end

function print_result(stop_reason::Int, time_::Float64, residuals::Residuals, max_rank::Int,  p::Params)
    println("---------------------------------------------------------------------------------------")
    println("    Solver status:")
    println("       "*p.stop_reason_string)
    println("       Time elapsed     = $(round(time_; digits = 2)) seconds")
    println("       Primal objective = $(round(residuals.prim_obj[p.iter]; digits = 5))")
    println("       Dual objective   = $(round(residuals.dual_obj[p.iter]; digits = 5))")
    println("       Duality gap      = $(round(100*residuals.dual_gap[p.iter]; digits = 2)) %")
    println("---------------------------------------------------------------------------------------")
    println("    Primal feasibility:")
    println("       ||A(X) - b|| / (1 + ||b||) = $(round(residuals.equa_feasibility; digits = 6))    [linear equalities] ")
    println("       ||max(G(X) - h, 0)|| / (1 + ||h||) = $(round(residuals.ineq_feasibility; digits = 6))    [linear inequalities]")
    println("    Rank of p.s.d. variable is $max_rank.")
    println("=======================================================================================")

    return nothing
end


text:
include(joinpath(@__DIR__, "testcommon.jl"))

# Wrap the functions defined in C++
module ParametricTypes

using CxxWrap
@wrapmodule(CxxWrap.CxxWrapCore.libparametric)

function __init__()
    @initcxx
end

end

import .ParametricTypes.TemplateType, .ParametricTypes.NonTypeParam

p1 = TemplateType{ParametricTypes.P1, ParametricTypes.P2}()
p2 = TemplateType{ParametricTypes.P2, ParametricTypes.P1}()

println("Dumping object p1:")
dump(p1)

@testset "$(basename(@__FILE__)[1:end-3])" begin

@test ParametricTypes.get_first(p1) == 1
@test ParametricTypes.get_second(p2) == 1
@test typeof(ParametricTypes.get_first(p1)) == Int32
@test typeof(ParametricTypes.get_second(p2)) == Int32

@test ParametricTypes.get_first(CxxRef(p1)) == 1
@test ParametricTypes.get_first(CxxRef(p1)[]) == 1
@test length(typeof(CxxRef(p1)[]).parameters) == 2

@test ParametricTypes.get_first(p2) == 10.
@test ParametricTypes.get_second(p1) == 10.
@test typeof(ParametricTypes.get_first(p2)) == Float64
@test typeof(ParametricTypes.get_second(p1)) == Float64

@test ParametricTypes.TemplateDefaultType{ParametricTypes.P1}() != nothing

nontype1 = ParametricTypes.NonTypeParam{Int32, Int32(1)}()
@test ParametricTypes.get_nontype(nontype1) == 1

nontype2 = ParametricTypes.NonTypeParam{UInt32, UInt32(2)}()
@test ParametricTypes.get_nontype(nontype2) == UInt32(2)

nontype3 = ParametricTypes.NonTypeParam{Int32, Int32(1)}(3)
@test ParametricTypes.get_nontype(nontype3) == 3

nontype4 = ParametricTypes.NonTypeParam{CxxLong, CxxLong(64)}()
@test ParametricTypes.get_nontype(nontype4) == CxxLong(64)

concr = ParametricTypes.ConcreteTemplate{Float64}()
@test isa(concr, ParametricTypes.AbstractTemplate{Float64})
@test isa(concr, ParametricTypes.AbstractTemplate)
@test isa(concr, ParametricTypes.ConcreteTemplate)
abst = ParametricTypes.to_base(CxxPtr(concr))
@test isa(abst[], ParametricTypes.AbstractTemplate{Float64})
@test isa(abst[], ParametricTypes.AbstractTemplate)

f3 = ParametricTypes.Foo3{Int32, CxxWrap.CxxBool, Float32}()
@test length(methods(ParametricTypes.foo3_method)) == 6
f2 = ParametricTypes.Foo2{Float64}()
@test length(methods(ParametricTypes.foo2_method)) == 2

@test length(methods(ParametricTypes.foo3_free_method)) == 6
ParametricTypes.foo3_free_method(f3)

@test supertype(ParametricTypes.Foo3{Float64,ParametricTypes.P1,Float32}) == ParametricTypes.AbstractTemplate{Float64}

darr = [1.0, 2.0, 3.0]
carr = Complex{Float32}[1+2im, 3+4im]
vec1 = ParametricTypes.CppVector{Float64}(pointer(darr), 3)
vec2 = ParametricTypes.CppVector2{Float64, Float32}()
vec3 = ParametricTypes.CppVector{Complex{Float32}}(pointer(carr), 2)
@test isa(vec1, AbstractVector{Float64})
@test isa(vec2, AbstractVector{Float64})
@test isa(vec3, AbstractVector{Complex{Float32}})
@test ParametricTypes.get(vec1,0)[] == 1.0
@test ParametricTypes.get(vec1,1)[] == 2.0
@test ParametricTypes.get(vec1,2)[] == 3.0
@test ParametricTypes.get(vec3,0)[] == 1+2im
@test ParametricTypes.get(vec3,1)[] == 3+4im

end


text:
using QualitativeAudioFeatures
using Test

@testset "QualitativeAudioFeatures.jl" begin
    # Write your own tests here.
end


text:
#This program is to test if two Nodes exist
#with only one objective functions.
using Plasmo
using Ipopt
using PlasmoTest

G = OptiGraph()
O = Ipopt.Optimizer

#Add nodes to a ModelGraph
n1 = @optinode(G)
n2 = @optinode(G)

#Node 1 Model
@variable(n1,0 <= x1 <= 2)
@variable(n1,2 <= y1 <= 4)
@constraint(n1,x1 + y1 >= 5)

#Node 2 Model
@variable(n2,x2)
@NLnodeconstraint(n2,ref,exp(x2) >= 2)
@variable(n2,y2 >= 0)
@constraint(n2,x2 + y2 >= 3)

#Add linkconstraint
@linkconstraint(G,n1[:x1]==n2[:y2])

#Graph objective
@objective(G, Min, n1[:x1] + n2[:x2])

#Use Ipopt to solve
optimize!(G,O)

println("x1 = ",digit5(value(n1,n1[:x1])))
println("y1 = ",digit5(value(n1,n1[:y1])))
println("x2 = ",digit5(value(n2,n2[:x2])))
println("y2 = ",digit5(value(n2,n2[:y2])))
println("Objective = ", digit5(objective_value(G)))

max_stars_repo_name: BoyangSheng/Skill-studing

max_issues_repo_name: BoyangSheng/Skill-studying

max_forks_repo_name: BoyangSheng/Skill-studying


text:
# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

"""
    VectorizeBridge{T,S} <: Bridges.Variable.AbstractBridge

`VectorizeBridge` implements the following reformulations:

 * ``x \\ge a`` into ``[y] \\in \\mathbb{R}_+`` with the substitution
   rule ``x = a + y``
 * ``x \\le a`` into ``[y] \\in \\mathbb{R}_-`` with the substitution
   rule ``x = a + y``
 * ``x == a`` into ``[y] \\in \\{0\\}`` with the substitution
   rule ``x = a + y``

where `T` is the coefficient type of `a + y`.

## Source node

`VectorizeBridge` supports:

 * [`MOI.VariableIndex`](@ref) in [`MOI.GreaterThan{T}`](@ref)
 * [`MOI.VariableIndex`](@ref) in [`MOI.LessThan{T}`](@ref)
 * [`MOI.VariableIndex`](@ref) in [`MOI.EqualTo{T}`](@ref)

## Target nodes

`VectorizeBridge` creates:

 * One variable node: [`MOI.VectorOfVariables`](@ref) in `S`, where `S` is one
   of [`MOI.Nonnegatives`](@ref), [`MOI.Nonpositives`](@ref),
   [`MOI.Zeros`](@ref) depending on the type of ``S``.
"""
mutable struct VectorizeBridge{T,S} <: AbstractBridge
    variable::MOI.VariableIndex
    vector_constraint::MOI.ConstraintIndex{MOI.VectorOfVariables,S}
    set_constant::T # constant in scalar set
end

const Vectorize{T,OT<:MOI.ModelLike} =
    SingleBridgeOptimizer{VectorizeBridge{T},OT}

function bridge_constrained_variable(
    ::Type{VectorizeBridge{T,S}},
    model::MOI.ModelLike,
    set::MOI.Utilities.ScalarLinearSet{T},
) where {T,S}
    set_constant = MOI.constant(set)
    variables, vector_constraint = MOI.add_constrained_variables(model, S(1))
    return VectorizeBridge{T,S}(variables[1], vector_constraint, set_constant)
end

function supports_constrained_variable(
    ::Type{VectorizeBridge{T}},
    ::Type{<:MOI.Utilities.ScalarLinearSet{T}},
) where {T}
    return true
end

function MOI.Bridges.added_constrained_variable_types(
    ::Type{VectorizeBridge{T,S}},
) where {T,S}
    return Tuple{Type}[(S,)]
end

function MOI.Bridges.added_constraint_types(::Type{<:VectorizeBridge})
    return Tuple{Type,Type}[]
end

function concrete_bridge_type(
    ::Type{<:VectorizeBridge{T}},
    S::Type{<:MOI.Utilities.ScalarLinearSet{T}},
) where {T}
    return VectorizeBridge{T,MOI.Utilities.vector_set_type(S)}
end

# Attributes, Bridge acting as a model
MOI.get(::VectorizeBridge, ::MOI.NumberOfVariables)::Int64 = 1

function MOI.get(bridge::VectorizeBridge, ::MOI.ListOfVariableIndices)
    return [bridge.variable]
end

function MOI.get(
    ::VectorizeBridge{T,S},
    ::MOI.NumberOfConstraints{MOI.VectorOfVariables,S},
)::Int64 where {T,S}
    return 1
end

function MOI.get(
    bridge::VectorizeBridge{T,S},
    ::MOI.ListOfConstraintIndices{MOI.VectorOfVariables,S},
) where {T,S}
    return [bridge.vector_constraint]
end

# References
function MOI.delete(model::MOI.ModelLike, bridge::VectorizeBridge)
    MOI.delete(model, bridge.variable)
    return
end

# Attributes, Bridge acting as a constraint

function MOI.get(
    ::MOI.ModelLike,
    ::MOI.ConstraintSet,
    bridge::VectorizeBridge{T,S},
) where {T,S}
    return MOI.Utilities.scalar_set_type(S, T)(bridge.set_constant)
end

function MOI.set(
    ::MOI.ModelLike,
    attr::MOI.ConstraintSet,
    bridge::VectorizeBridge,
    ::MOI.Utilities.ScalarLinearSet,
)
    # This would require modifing any constraint which uses the bridged
    # variable.
    return throw(
        MOI.SetAttributeNotAllowed(
            attr,
            "The variable `$(bridge.variable)` is bridged by the " *
            "`VectorizeBridge`.",
        ),
    )
end

function MOI.get(
    model::MOI.ModelLike,
    attr::MOI.ConstraintPrimal,
    bridge::VectorizeBridge,
)
    x = MOI.get(model, attr, bridge.vector_constraint)
    @assert length(x) == 1
    y = x[1]
    status = MOI.get(model, MOI.PrimalStatus(attr.result_index))
    if !MOI.Utilities.is_ray(status)
        # If it is an infeasibility certificate, it is a ray and satisfies the
        # homogenized problem, see https://github.com/jump-dev/MathOptInterface.jl/issues/433
        # Otherwise, we need to add the set constant since the ConstraintPrimal
        # is defined as the value of the function and the set_constant was
        # removed from the original function
        y += bridge.set_constant
    end
    return y
end

function MOI.get(
    model::MOI.ModelLike,
    attr::MOI.ConstraintDual,
    bridge::VectorizeBridge,
)
    x = MOI.get(model, attr, bridge.vector_constraint)
    @assert length(x) == 1
    return x[1]
end

function MOI.get(
    model::MOI.ModelLike,
    attr::MOI.VariablePrimal,
    bridge::VectorizeBridge,
)
    value = MOI.get(model, attr, bridge.variable)
    status = MOI.get(model, MOI.PrimalStatus(attr.result_index))
    if MOI.Utilities.is_ray(status)
        return value  # Return the homogenized value.
    end
    return value + bridge.set_constant
end

function MOI.get(
    model::MOI.ModelLike,
    attr::MOI.VariablePrimalStart,
    bridge::VectorizeBridge,
)
    return MOI.get(model, attr, bridge.variable) + bridge.set_constant
end

function MOI.supports(
    model::MOI.ModelLike,
    attr::MOI.VariablePrimalStart,
    ::Type{<:VectorizeBridge},
)
    return MOI.supports(model, attr, MOI.VariableIndex)
end

function MOI.set(
    model::MOI.ModelLike,
    attr::MOI.VariablePrimalStart,
    bridge::VectorizeBridge,
    value,
)
    MOI.set(model, attr, bridge.variable, value - bridge.set_constant)
    return
end

function MOI.Bridges.bridged_function(bridge::VectorizeBridge{T}) where {T}
    return MOI.Utilities.operate(+, T, bridge.variable, bridge.set_constant)
end

function unbridged_map(
    bridge::VectorizeBridge{T},
    vi::MOI.VariableIndex,
) where {T}
    func = MOI.Utilities.operate(-, T, vi, bridge.set_constant)
    return (bridge.variable => func,)
end


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    AbstractGeometry{T,N}

A geometry in a `N`-dimensional space with coordinates of type `T`.
"""
abstract type AbstractGeometry{T,N} end

"""
    x ∈ geometry

Check if coordinates `x` are in the `geometry`.
"""
in

#------------------
# IMPLEMENTATIONS
#------------------
include("geometries/rectangle.jl")


text:
module LinearMaps

export AbstractLinearMap, LinearMap

import Base: +, -, *, \, /, ==

abstract AbstractLinearMap{T}
Base.eltype{T}(::AbstractLinearMap{T})=T
Base.eltype{T}(::Type{AbstractLinearMap{T}})=T
Base.eltype{L<:AbstractLinearMap}(::Type{L})=eltype(super(L))

Base.isreal{T<:Real}(::AbstractLinearMap{T})=true
Base.isreal(::AbstractLinearMap)=false # standard assumptions
Base.issym(::AbstractLinearMap)=false # standard assumptions
Base.ishermitian{T<:Real}(A::AbstractLinearMap{T})=issym(A)
Base.ishermitian(::AbstractLinearMap)=false # standard assumptions
Base.isposdef(::AbstractLinearMap)=false # standard assumptions

Base.size(A::AbstractLinearMap,n)=(n==1 || n==2 ? size(A)[n] : error("AbstractLinearMap objects have only 2 dimensions"))

# any AbstractLinearMap subtype will have to overwrite at least one of the two following methods to avoid running in circles
*(A::AbstractLinearMap,x::AbstractVector)=Base.A_mul_B!(similar(x,promote_type(eltype(A),eltype(x)),size(A,1)),A,x)
Base.A_mul_B!(y::AbstractVector,A::AbstractLinearMap,x::AbstractVector)=begin
    length(y)==size(A,1) || throw(DimensionMismatch("A_mul_B!"))
    copy!(y,A*x)
end

# the following for multiplying with transpose and ctranspose map are optional:
# subtypes can overwrite nonmutating methods, implement mutating methods or do nothing
Base.At_mul_B(A::AbstractLinearMap,x::AbstractVector)=(@which Base.At_mul_B!(x,A,x))!=methods(Base.At_mul_B!,(AbstractVector,AbstractLinearMap,AbstractVector))[end] ?
    Base.At_mul_B!(similar(x,promote_type(eltype(A),eltype(x)),size(A,2)),A,x) : throw(MethodError(Base.At_mul_B,(A,x)))
Base.At_mul_B!(y::AbstractVector,A::AbstractLinearMap,x::AbstractVector)=begin
    length(y)==size(A,2) || throw(DimensionMismatch("At_mul_B!"))
    (@which Base.At_mul_B(A,x))!=methods(Base.At_mul_B,(AbstractLinearMap,AbstractVector))[end] ? copy!(y,Base.At_mul_B(A,x)) : throw(MethodError(Base.At_mul_B!,(y,A,x)))
end
Base.Ac_mul_B(A::AbstractLinearMap,x::AbstractVector)=(@which Base.Ac_mul_B!(x,A,x))!=methods(Base.Ac_mul_B!,(AbstractVector,AbstractLinearMap,AbstractVector))[end] ?
    Base.Ac_mul_B!(similar(x,promote_type(eltype(A),eltype(x)),size(A,2)),A,x) : throw(MethodError(Base.Ac_mul_B,(A,x)))
Base.Ac_mul_B!(y::AbstractVector,A::AbstractLinearMap,x::AbstractVector)=begin
    length(y)==size(A,2) || throw(DimensionMismatch("At_mul_B!"))
    (@which Base.Ac_mul_B(A,x))!=methods(Base.Ac_mul_B,(AbstractLinearMap,AbstractVector))[end] ? copy!(y,Base.Ac_mul_B(A,x)) : throw(MethodError(Base.Ac_mul_B!,(y,A,x)))
end

# full: create matrix representation of AbstractLinearMap
function Base.full{T}(A::AbstractLinearMap{T})
    M,N=size(A)
    mat=zeros(T,(M,N))
    for i=1:N
        v=zeros(T,N)
        v[i]=one(T)
        w=pointer_to_array(pointer(mat,(i-1)*M+1),M)
        A_mul_B!(w,A,v)
    end
    return mat
end

include("transpose.jl") # transposing linear maps
include("linearcombination.jl") # defining linear combinations of linear maps
include("composition.jl") # composition of linear maps
include("wrappedmap.jl") # wrap a matrix of linear map in a new type, thereby allowing to alter its properties
include("identitymap.jl") # the identity map, to be able to make linear combinations of AbstractLinearMap objects and I
include("functionmap.jl") # using a function as linear map

LinearMap{T}(A::Union{AbstractMatrix{T},AbstractLinearMap{T}};isreal::Bool=Base.isreal(A),issym::Bool=Base.issym(A),ishermitian::Bool=Base.ishermitian(A),isposdef::Bool=Base.isposdef(A)) =
    WrappedMap(A;isreal=isreal,issym=issym,ishermitian=ishermitian,isposdef=isposdef)
LinearMap(f::Function,M::Int,N::Int=M;ismutating::Bool=false,isreal::Bool=true,issym::Bool=false,ishermitian::Bool=(isreal && issym),isposdef::Bool=false,ftranspose::OptionalFunction=nothing,fctranspose::OptionalFunction=nothing) =
    FunctionMap(f,M,N;ismutating=ismutating,isreal=isreal,issym=issym,ishermitian=ishermitian,isposdef=isposdef,ftranspose=ftranspose,fctranspose=fctranspose)
LinearMap(f::Function,eltype::Type,M::Int,N::Int=M;ismutating::Bool=false,issym::Bool=false,ishermitian::Bool=(eltype<:Real && issym),isposdef::Bool=false,ftranspose::OptionalFunction=nothing,fctranspose::OptionalFunction=nothing) =
    FunctionMap{eltype}(f,M,N;ismutating=ismutating,issym=issym,ishermitian=ishermitian,isposdef=isposdef,ftranspose=ftranspose,fctranspose=fctranspose)

end # module


text:
using Revise
using FlippingModel
using DataFrames, Optim
# using BrowseTables
using Test
import Distributions: pdf

##
pokes
names(pokes)
pokes[!,:Leave] = Vector{Union{Missing,Bool}}(missing,length(pokes.Side))
pokes[!,:Bout] .= 0
bouts = by(pokes, :Session) do dd
    process_bouts(dd)
end

b = filter(bouts) do row
    !ismissing(row.Leave) &&
    row.Correct_start &&
    row.Start_Reward &&
    row.Streak > 10 &&
    3 < row.Streak_within_Block < 23 &&
    row.ReverseStreak > 20 &&
    (ismissing(row.Pre_Interpoke) || row.Pre_Interpoke < 1) &&
    !ismissing(row.Leave) &&
    row.Protocol_Day > 1 &&
    row.Day > 20191215
end
disallowmissing!(b, :Leave)
##
using CSV
file = "/Volumes/GoogleDrive/My Drive/Flipping/Datasets/Stimulations/WT_new_protocols/Results/bouts.csv"
CSV.write(file,b)
# open_html_table(b[end-100:end,:])


text:
#####################################################
# Helper functions for vectorize/reconstruct values #
#####################################################

vectorize(d::UnivariateDistribution, r::Real) = [r]
vectorize(d::MultivariateDistribution, r::AbstractVector{<:Real}) = copy(r)
vectorize(d::MatrixDistribution, r::AbstractMatrix{<:Real}) = copy(vec(r))

# NOTE:
# We cannot use reconstruct{T} because val is always Vector{Real} then T will be Real.
# However here we would like the result to be specifric type, e.g. Array{Dual{4,Float64}, 2},
# otherwise we will have error for MatrixDistribution.
# Note this is not the case for MultivariateDistribution so I guess this might be lack of
# support for some types related to matrices (like PDMat).
reconstruct(d::UnivariateDistribution, val::AbstractVector) = val[1]
reconstruct(d::MultivariateDistribution, val::AbstractVector) = copy(val)
function reconstruct(d::MatrixDistribution, val::AbstractVector)
    return copy(reshape(val, size(d)))
end
function reconstruct!(r, d::Distribution, val::AbstractVector)
    return reconstruct!(r, d, val)
end
function reconstruct!(r, d::MultivariateDistribution, val::AbstractVector)
    r .= val
    return r
end
function reconstruct(d::Distribution, val::AbstractVector, n::Int)
    return reconstruct(d, val, n)
end
function reconstruct(d::UnivariateDistribution, val::AbstractVector, n::Int)
    return copy(val)
end
function reconstruct(d::MultivariateDistribution, val::AbstractVector, n::Int)
    return copy(reshape(val, size(d)[1], n))
end
function reconstruct(d::MatrixDistribution, val::AbstractVector, n::Int)
    tmp = reshape(val, size(d)[1], size(d)[2], n)
    orig = [tmp[:, :, i] for i in 1:size(tmp, 3)]
    return orig
end
function reconstruct!(r, d::Distribution, val::AbstractVector, n::Int)
    return reconstruct!(r, d, val, n)
end
function reconstruct!(r, d::MultivariateDistribution, val::AbstractVector, n::Int)
    r .= val
    return r
end


text:
module ProximalSEM

    using StructuralEquationModels
    import StructuralEquationModels: print_type_name, print_field_types, sem_fit, start_val
    import ProximalCore, ProximalAlgorithms, ProximalOperators
    import ProximalCore: prox!

    ProximalCore.prox!(y, f, x, gamma) = ProximalOperators.prox!(y, f, x, gamma)

    include("diff/Proximal.jl")
    include("optimizer/ProximalAlgorithms.jl")

    export SemDiffProximal, sem_fit

end


text:
function plot_all(mha,mhb;hidden=false)
    colors=["g","r"]
    figure()
    for (i,mh) in enumerate((mha, mhb))
        bestparams = mh.mapvalue.curr
        bestpred = pred(bestparams.theta, bestparams.init1, bestparams.init2)
        for j in 1:size(mh.obj.data, 1)
            if j == 1
                plot(mh.obj.data[j,:]', colors[i], alpha=0.3, linewidth=2, label="Class $i")
            elseif j == 2
                plot(mh.obj.data[j,:]', colors[i], alpha=0.3, linewidth=2, label="Class $i")
            else
                plot(mh.obj.data[j,:]', colors[i], alpha=0.3, linewidth=2)
            end
        end
        if !hidden
            plot(bestpred, colors[i], linewidth=5)
        end
    end
end

function plot_trace(mha, mhb, name)
    colors=["g","r"]
    figure()
    for (i,mh) in enumerate((mha, mhb))
        accum = Float64[]
        for params in mh.db
            push!(accum, params.(symbol(name)))
        end
        plot(accum, colors[i])
        println("Mean: $(mean(accum)), stddev: $(std(accum))")
    end
end

function plot_posterior(mh)
    figure()
    for params in mh.db
        plot(pred(params.theta, params.init1, params.init2), "k")
    end
end


text:
# Basic tests of poses at zero

using RoME
using Test
using TensorCast

@testset "basic pose2 trivial case without forcing autoinit..." begin

fg = initfg()

cov = Matrix(Diagonal(1e-4.*[1;1;1]))

# first pose position
addVariable!(fg, :x0, Pose2)
addFactor!(fg, [:x0], PriorPose2(MvNormal(zeros(3),cov)))

# second Pose position with odo
addVariable!(fg, :x1, Pose2)
addFactor!(fg, [:x0;:x1], Pose2Pose2(MvNormal(zeros(3),cov)))

# third pose position with odo
addVariable!(fg, :x2, Pose2)
addFactor!(fg, [:x1;:x2], Pose2Pose2(MvNormal(zeros(3),cov)))

#

Xc_badval = 0.01.*randn(3,100)
Xc_badval[1,:] .-= 5.0
Xc_badval[2,:] .-= 2.0
Xc_badval[3,:] .+= 0.5

badval = map((Xc)->DFG.getPoint(Pose2, Xc), eachcol(Xc_badval))

setValKDE!(fg, :x2, manikde!(badval, Pose2))

N = 100
# batchSolve!(fg, N=N)
tree = solveTree!(fg)

for xx in [:x0; :x1; :x2]
  pts = getVal(fg, xx)
  pts_μ = mean(getManifold(Pose2), pts)
  isapprox(getManifold(Pose2), pts_μ, ProductRepr([0.,0], [1. 0; 0 1]), atol=0.01)

  Xpts = getCoordinates.(Ref(Pose2), pts)
  @cast Xptsarr[j,i] := Xpts[i][j]

  @test 0.95*N < sum( -1.0 .< Xptsarr[1,:] .< 1.0 )
  @test 0.95*N < sum( -1.0 .< Xptsarr[2,:] .< 1.0 )
  @test 0.95*N < sum( -0.5 .< Xptsarr[3,:] .< 0.5 )
end

end

@testset "basic pose2 with forcing bad initialization..." begin

fg = initfg()

cov = Matrix(Diagonal(1e-4.*[1;1;1]))

# first pose position
addVariable!(fg, :x0, Pose2)
addFactor!(fg, [:x0], PriorPose2(MvNormal(zeros(3),cov)))

# second Pose position with odo
addVariable!(fg, :x1, Pose2)
addFactor!(fg, [:x0;:x1], Pose2Pose2(MvNormal(zeros(3),cov)))

# third pose position with odo
addVariable!(fg, :x2, Pose2)
addFactor!(fg, [:x1;:x2], Pose2Pose2(MvNormal(zeros(3),cov)))

ensureAllInitialized!(fg)

Xc_badval = 0.000001.*randn(3,100)
Xc_badval[1,:] .-= 5.0
Xc_badval[2,:] .-= 2.0
Xc_badval[3,:] .+= 0.5
badval = map((Xc)->DFG.getPoint(Pose2, Xc), eachcol(Xc_badval))

setValKDE!(fg, :x2, manikde!(badval, Pose2))

# tree = wipeBuildNewTree!(fg, drawpdf=true, show=true)

N = 100
getSolverParams(fg).N = N
solveTree!(fg)

for xx in [:x0; :x1; :x2]
  manipts = getVal(fg, xx)
  manicrd = getCoordinates.(Ref(Pose2), manipts)
  @cast pts[j,i] := manicrd[i][j]

  @test 0.95*N < sum( -0.5 .< pts[1,:] .< 0.5 )
  @test 0.95*N < sum( -1.0 .< pts[2,:] .< 1.0 )
  @test 0.95*N < sum( -1.0 .< pts[3,:] .< 1.0 )
end

end

# using RoMEPlotting
#
# drawPoses(fg, spscale=0.5)
#
# spyCliqMat(tree, :x2)

@testset "test basic banana (split)..." begin

fg = initfg()

addVariable!(fg, :x0, Pose2())
addVariable!(fg, :x1, Pose2())

addFactor!(fg,
           [:x0],
           PriorPose2(MvNormal(zeros(3), Matrix(Diagonal([0.01;0.01;1.0].^2)))),
           graphinit=false )
#

addFactor!(fg,
           [:x0;:x1],
           Pose2Pose2(MvNormal([1.0;0;0], Matrix(Diagonal([0.01;0.01;0.01].^2)))),
           graphinit=false )
#

solveTree!(fg)

manipts = getPoints(getKDE(fg, :x0))
manicrd = getCoordinates.(Ref(Pose2), manipts)
@cast pts[j,i] := manicrd[i][j]

N = size(pts,2)

@test 0.7*N < sum(abs.(pts[1,:]) .< 0.1)
@test 0.7*N < sum(abs.(pts[2,:]) .< 0.1)
@test 0.7*N < sum(abs.(pts[3,:]) .< 2.0)

manipts = getPoints(getKDE(fg, :x1))
manicrd = getCoordinates.(Ref(Pose2), manipts)
@cast pts[j,i] := manicrd[i][j]

@test 0.7*N < sum(0.0 .< pts[1,:])
@test 0.7*N < sum( 0.5 .< sqrt.(sum(pts[1:2,:].^2, dims=1)) .< 1.5)
@test 0.7*N < sum(abs.(pts[3,:]) .< 2.0)

#
# using RoMEPlotting
# using Gadfly
# Gadfly.set_default_plot_size(35cm,25cm)
#
# plotKDE(fg, ls(fg))

end


text:
mutable struct Vertex
    pos::Array{Float64}
end

Base.copy(s::Vertex) = Vertex(copy(s.pos))

mutable struct Plane
    normal::Array{Float64}
    w::Float64
end

mutable struct Polygon
    vertices::Array{Vertex,1}
    plane::Plane
end

mutable struct Node
    plane::Union{Nothing,Plane}
    front::Union{Nothing,Node}
    back::Union{Nothing,Node}
    polygons::Array{Polygon,1}
end

mutable struct Solid
    vertices::Array{Vertex,1}
    indices::Array{Int64,1}
end


text:
using Base.Test
using MicroLogging

configure_logging(min_level=:warn)

@testset "DeformableRegistration" begin

@time @testset "submodules" begin
    include("ImageProcessing.jl")
    include("Transformation.jl")
    include("Interpolation.jl")
    include("Distance.jl")
    include("Regularizer.jl")
end

# test full parametric and nonparametric image registration
@time @testset "example registration cases" begin
    include("ParametricRegistration.jl")
    include("NonparametricRegistration.jl")
end

end


text:
#=
    Utility functions to handle I/O graph from edgefiles.
=#

"""
    Create graph from CSV file.

    The file must contain as header the columns: 'Source', 'Target'
    and (optional)'Type'. The indexes of the vertices must be
    integers and withing the range [1,N], where N is the total
    number of vertices in the graph.

    Args:
        file_path:
            String containing the path to the CSV file.
        is_directed:
            {true, false} boolean value to signal whether the
            network is directed(true) or undirected(false).
    Return:
        graph:
            'Graph' structure representing the CSV edgelist file. This
            object contains at least the integer edge property called
            'edgetype'.
"""
function graph_from_csv(file_path::String, is_directed::Bool)
    df = CSV.File(file_path) |> DataFrame
    graph = Graph(is_directed)
    N = maximum([maximum(df.Source), maximum(df.Target)])

    # Create the vertices
    for i in 1:N
        add_node(graph)
    end
    # Add the connections
    for j in 1:length(df.Source)
        add_edge(df.Source[j], df.Target[j], graph)
    end
    if "Type" in names(df)
        if typeof(df.Type[1])==Int
            set_edges_properties("edgetype", df.Type, graph)
        # Converts the string types to integer type.
        elseif typeof(df.Type[1])==String
            unique_types = collect(Set(df.Type))
            label_map = Dict{String, Int}()
            for j in 1:length(unique_types)
                label_map[unique_types[j]] = j
            end
            new_int_edgeprop = [label_map[j] for j in df.Type]
            types = [ df.Type[j] for j in 1:length(df.Type)]
            set_edges_properties("edgetype", new_int_edgeprop, graph)
            set_edges_properties("edgetype", types, graph)
        else
            error("Format of 'Type' column is not accepted.")
        end
    else
        set_edges_properties("edgetype", [1 for i in 1:graph.M], graph)
    end
    return graph
end

"""
    Create graph from edgelist.
    edgelist is an array of dimension M x 2 containing
    the source and target of each edge.
    The function assumes that the vertices are labeled
    from 1 to N.
"""
function graph_from_edgelist(edgelist::Array{Int, 2}, is_directed::Bool)
    graph = Graph(is_directed)
    N = maximum([maximum(edgelist[:,1]), maximum(edgelist[:,2])])

    # -- Create the vertices --
    for j in 1:N
        add_node(graph)
    end
    # -- Then, build the connections
    for j in 1:length(edgelist[:,1])
        add_edge(edgelist[j,1], edgelist[j,2], graph)
    end
    return graph
end

"""
    Given a formatted edgelist file, it returns two objects: a N x 2
    array containing the edgelist with format [[src,tgt], ...] and
    a array of string arrays containing all the other columns in the file.

    Example:
    file -> "1 2 prop1 prop2
             2 3 prop1 prop2
             ... ..... ....."

    returns ->
    [1 2; 2 3; ...] and [["prop1", "prop2"], ["prop1", "prop2"], ...]
"""
function process_edgefile(fname::String, convert_int=false)
    src_tgt = Array{String}[]
    edge_prop = Array{String}[]

    open(fname, "r") do f
        for line in eachline(f)
            line_elem = split(line)
            elem1 = [ line_elem[1], line_elem[2] ]
            #elem1 = Int[parse(Int, line_elem[1]), parse(Int, line_elem[2])]

            append!(src_tgt, [elem1])
            append!(edge_prop, [line_elem[3:length(line_elem)]])
        end
    end

    src_tgt = reduce(hcat, src_tgt)
    edges = permutedims(src_tgt)
    if convert_int
        edges = parse.(Int, edges)
    end
    return edges, edge_prop
end

"""
    Process the edge properties returned by function 'process_edgefile'.

    As returned from 'process_edgefile', eprops is an array containing
    string arrays with the properties' values for each edge.

    For example, 'eprops' holds two string arrays if the edgefile has two
    extra columns. The size of these arrays are equal to the number of edges.
    For this, 'names' represents the arrays with the string names for each
    column.

    Returns a dictionary where the keys are the names in 'names' and the
    values are arrays with the edge properties indexed from 1 to M.
"""
function process_eprops(eprops::Array{Array{String}}, names::Array{String})
    container = Array{String}[]
    n_props = length(eprops[1])
    for j in 1:n_props
        append!(container, [String[]])
    end
    for j in 1:length(eprops)
        for k in 1:n_props
            append!(container[k], [eprops[j][k]])
        end
    end

    holder = Dict{String, Array{String}}()
    for (j, name) in enumerate(names)
        holder[name] = container[j]
    end
    return holder
end

"""

"""
function create_indexing(edges::Array{String,2})
    unroll = reduce(vcat, edges)
    unroll = collect(Set(unroll))
    N = length(unroll)

    hash = Dict{String, Int}()
    for (j, name) in enumerate(unroll)
        hash[name] = j
    end

    M, dummy = size(edges)
    new_edges = Array{Int}[]
    for j in 1:M
        u = hash[edges[j,1]]
        v = hash[edges[j,2]]
        append!(new_edges, [[u, v]])
    end
    new_edges = copy(transpose(reduce(hcat, new_edges)))
    return new_edges, hash
end

"""
    Read edgelist, create and prepare the graph for the application of the fast
    fibration algorithm.

    The expected edgelist has 3 columns: the source column, the target column and
    the edgetype column.

    The source and target must be integers ranging from 1 to N, where N is the total
    number of nodes of the graph. The third column can be strings with the unique
    values for each type. These edgetypes will be set as edge integer properties
    ranging from 1 to K, where K is the total number of unique types.

    Returns a graph object 'graph' containing the integer edge property called
    'edgetype'.
"""
function load_net(fname::String, is_directed::Bool, convert_int=false, etype_col="Column 1")
    edges, eprops = process_edgefile(fname, convert_int)
    # -- if the source and target columns need to be transformed to the interval 1 to
    # -- N, then we create a mapping 'name_map' for each unique value of the nodes.
    if !convert_int
        edges, name_map = create_indexing(edges)
    end
    # -- Save the original string indexes of nodes as 'node_name' vertex property --
    graph = graph_from_edgelist(edges, is_directed)
    N = graph.N
    nodes_name = [ "" for j in 1:N ]
    if !convert_int
        keys_name = collect(keys(name_map))
        for key in keys_name
            nodes_name[name_map[key]]*=key
        end
        set_vertices_properties("node_name", nodes_name, graph)
    end

    if length(eprops[1])!=0
        col_names = [ "Column $j" for j in 1:length(eprops[1]) ]
        fmt_eprops = process_eprops(eprops, col_names)
    end

    edgetype_str = collect(String, Set(fmt_eprops[etype_col]))
    edgetype_map = Dict{String, Int}()
    for (j, value) in enumerate(edgetype_str)
        edgetype_map[value] = j
    end

    int_edgetypes = Int[]
    for (j, etype_str) in enumerate(fmt_eprops[etype_col])
        push!(int_edgetypes, edgetype_map[etype_str])
    end
    set_edges_properties("edgetype", int_edgetypes, graph)

    return graph, fmt_eprops
end


text:
module PACKAGE_NAME

function hello(name="world")
    println("Hello, ", name, "!")
end

end # module PACKAGE_NAME

max_stars_repo_name: fredrikekre/test-julia-pkg-template

max_issues_repo_name: fredrikekre/test-julia-pkg-template

max_forks_repo_name: fredrikekre/test-julia-pkg-template


text:
export Data

"""
    Data( m :: Matrix )

Data structure to store the solution of an initial-value problem along time.

`data=Data(m)` is of parametric type and offers
- `data.U`, a 1-element vector with a copy of the matrix `m`;
- `(data.datalength,data.datasize)=size(m)`  where \
`datalength` is the number of computed modes, and `datasize` the number of involved equations, typically 2.
"""
mutable struct Data

    U :: Vector{Array{ComplexF64,2}}
    datasize :: Int
    datalength :: Int

    function Data( v :: AbstractArray )

        (datalength , datasize ) = size(v)

        U = [copy(v)]

        new(U, datasize, datalength)

    end

end

Base.:length( data :: Data ) = length( data.U )

Base.:size( data :: Data ) = size( first(data.U) )

Base.:(==)(d1::Data, d2::Data) =
    d1.datalength == d2.datalength &&
    d1.datasize == d2.datasize &&
    d1.U == d2.U


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons, Josef Kircher
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

using DifferentialEquations, DiffEqCallbacks, Zygote

############ Model-Exchange ############

"""
Handles events and returns the values and nominals of the changed continuous states.
"""
function handleEvents(c::fmi2Component, initialEventMode)

    if initialEventMode == false
        fmi2EnterEventMode(c)
    end

    eventInfo = fmi2NewDiscreteStates(c)

    eventInfo.newDiscreteStatesNeeded = fmi2True
    valuesOfContinuousStatesChanged = fmi2False
    nominalsOfContinuousStatesChanged = fmi2False

    #set inputs here
    #fmiSetReal(myFMU, InputRef, Value)

    while eventInfo.newDiscreteStatesNeeded == fmi2True
        # update discrete states
        eventInfo = fmi2NewDiscreteStates(c)
        valuesOfContinuousStatesChanged = eventInfo.valuesOfContinuousStatesChanged
        nominalsOfContinuousStatesChanged = eventInfo.nominalsOfContinuousStatesChanged

        if eventInfo.terminateSimulation == fmi2True
            @error "Event info returned error!"
        end
    end

    fmi2EnterContinuousTimeMode(c)

    return valuesOfContinuousStatesChanged, nominalsOfContinuousStatesChanged

end

"""
Returns the event indicators for an FMU.
"""
function condition(c::fmi2Component, out, x, t, integrator) # Event when event_f(u,t) == 0

    fmi2SetTime(c, t)
    fmi2SetContinuousStates(c, x)
    indicators = fmi2GetEventIndicators(c)

    copy!(out, indicators)

end

"""
Handles the upcoming events.
"""
function affect!(c::fmi2Component, integrator, idx)
    # Event found - handle it
    continuousStatesChanged, nominalsChanged = handleEvents(c, false)

    if continuousStatesChanged == fmi2True
        integrator.u = fmi2GetContinuousStates(c)
    end

    if nominalsChanged == fmi2True
        x_nom = fmi2GetNominalsOfContinuousStates(c)
    end
end

"""
Does one step in the simulation.
"""
function stepCompleted(c::fmi2Component, x, t, integrator)

     (status, enterEventMode, terminateSimulation) = fmi2CompletedIntegratorStep(c, fmi2True)
     if enterEventMode == fmi2True
        affect!(c, integrator, 0)
     end

end

"""
Returns the state derivatives of the FMU.
"""
function fx(c::fmi2Component, x, p, t)
    fmi2SetTime(c, t)
    fmi2SetContinuousStates(c, x)
    dx = fmi2GetDerivatives(c)
end

"""
Source: FMISpec2.0.2[p.90 ff]: 3.2.4 Pseudocode Example

Simulates an fmu instance for the set simulation time.
Events are handled correctly.

Returns an ODESolution.
"""
function fmi2SimulateME(c::fmi2Component, t_start::Real = 0.0, t_stop::Real = 1.0;
                        solver = nothing,
                        customFx = nothing,
                        recordValues::fmi2ValueReferenceFormat = nothing,
                        saveat = [],
                        setup = true)

    recordValues = prepareValueReference(c, recordValues)

    if length(recordValues) > 0
        @warn "fmi2SimulateME(...): Recording values via `recordValues` not implemented, yet."
    end

    if solver == nothing
        solver = Tsit5()
    end

    if customFx == nothing
        customFx = (x, p, t) -> fx(c, x, p, t)
    end

    if setup
        fmi2Reset(c)
        fmi2SetupExperiment(c, t_start, t_stop)
        fmi2EnterInitializationMode(c)
        fmi2ExitInitializationMode(c)
    end

    eventCb = VectorContinuousCallback((out, x, t, integrator) -> condition(c, out, x, t, integrator),
                                       (integrator, idx) -> affect!(c, integrator, idx),
                                       Int64(c.fmu.modelDescription.numberOfEventIndicators))

    stepCb = FunctionCallingCallback((x, t, integrator) -> stepCompleted(c, x, t, integrator);
                                     func_everystep=true,
                                     func_start=true)

     # First evaluation of the FMU
     x0 = fmi2GetContinuousStates(c)
     x0_nom = fmi2GetNominalsOfContinuousStates(c)

     fmi2SetContinuousStates(c, x0)
     handleEvents(c, true)

     # Get states of handling initial Events
     x0 = fmi2GetContinuousStates(c)
     x0_nom = fmi2GetNominalsOfContinuousStates(c)

     p = []
     problem = ODEProblem(customFx, x0, (t_start, t_stop), p,)
     solution = solve(problem, solver, callback=CallbackSet(eventCb, stepCb), saveat=saveat)
end

############ Co-Simulation ############

"""
Starts a simulation of the Co-Simulation FMU instance.
"""
function fmi2SimulateCS(c::fmi2Component, t_start::Real, t_stop::Real;
                        recordValues::fmi2ValueReferenceFormat = nothing,
                        saveat = [],
                        setup = true)

    recordValues = prepareValueReference(c, recordValues)
    variableSteps = c.fmu.modelDescription.CScanHandleVariableCommunicationStepSize

    # default setup
    if length(saveat) == 0
        saveat = LinRange(t_start, t_stop, 100)
    end
    dt = (t_stop - t_start) / length(saveat)

    # setup if no variable steps
    if variableSteps == false
        if length(saveat) >= 2
            dt = saveat[2] - saveat[1]
        end
    end

    if setup
        fmi2Reset(c)
        fmi2SetupExperiment(c, t_start, t_stop)
        fmi2EnterInitializationMode(c)
        fmi2ExitInitializationMode(c)
    end

    t = t_start

    sd = nothing
    record = length(recordValues) > 0

    #numDigits = length(string(round(Integer, 1/dt)))

    if record
        sd = fmi2SimulationResult()
        sd.valueReferences = recordValues
        sd.dataPoints = []
        sd.fmu = c.fmu

        values = fmi2GetReal(c, sd.valueReferences)
        push!(sd.dataPoints, (t, values...))

        i = 1
        while t < t_stop
            if variableSteps
                dt = saveat[i+1] - saveat[i]
            end

            fmi2DoStep(c, t, dt)
            t = t + dt #round(t + dt, digits=numDigits)
            i += 1

            values = fmi2GetReal(c, sd.valueReferences)
            push!(sd.dataPoints, (t, values...))
        end
    else
        i = 1
        while t < t_stop
            if variableSteps
                dt = saveat[i+1] - saveat[i]
            end

            fmi2DoStep(c, t, dt)
            t = t + dt #round(t + dt, digits=numDigits)
            i += 1
        end
    end

    sd
end

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

"""
Starts a simulation of the fmu instance for the matching fmu type, if both types are available the CoSimulation Simulation is started.
"""
function fmi2Simulate(c::fmi2Component, t_start::Real = 0.0, t_stop::Real = 1.0;
                      recordValues::fmi2ValueReferenceFormat = nothing,
                      saveat = [],
                      setup = true)

    if fmi2IsCoSimulation(c.fmu)
        return fmi2SimulateCS(c, t_start, t_stop; recordValues=recordValues, saveat=saveat, setup=setup)
    elseif fmi2IsModelExchange(c.fmu)
        return fmi2SimulateME(c, t_start, t_stop; recordValues=recordValues, saveat=saveat, setup=setup)
    else
        error(unknownFMUType)
    end
end


text:
module SVGDraw

include("style.jl")
include("shape.jl")

export R², R³, R⁴, BoundingBox, setproperty!, getproperty,
       VGElement, Rect, Circ, Polyline, Brace,
       lefttop, leftbottom, righttop, rightbottom,
       midtop, midbottom, midleft, midright,
       CssRules,
       make_svg_root, make_svg_ele
end


text:
################################################################################
##### CoveredRange
################################################################################

# a range covered by subranges
struct CoveredRange{R,S} <: AbstractVector{UnitRange{Int}}
    offsets::R
    stopping::S
end

struct FixedLength
    length::Int
end

struct LengthAtMost
    maxlength::Int
    maxstop::Int
end

function compute_stop(offset, stopping::FixedLength)
    return offset+stopping.length
end
function compute_stop(offset, stopping::LengthAtMost)
    return min(offset+stopping.maxlength, stopping.maxstop)
end

compute_range(offset, stopping)::UnitRange{Int} = (offset+1):compute_stop(offset, stopping)

Base.size(o::CoveredRange,args...) = size(o.offsets, args...)
Base.@propagate_inbounds function Base.getindex(o::CoveredRange, inds...)
    offset = o.offsets[inds...]
    return compute_range(offset, o.stopping)
end

################################################################################
##### RoundedRange
################################################################################
struct RoundedRange{R} <: AbstractVector{Int}
    range::R
end

Base.@propagate_inbounds function Base.getindex(r::RoundedRange, i)
    rough = r.range[i]
    return round(Int, rough)
end

function roundedrange(start; stop, length)
    inner = LinRange(start, stop, length)
    return RoundedRange(inner)
end

Base.size(o::RoundedRange, args...) = size(o.range, args...)

################################################################################
##### TileIterator
################################################################################
struct TileIterator{N,C} <: AbstractArray{NTuple{N, UnitRange{Int}}, N}
    covers1d::C
end

function TileIterator(covers1d::NTuple{N, AbstractVector{UnitRange{Int}}}) where {N}
    C = typeof(covers1d)
    return TileIterator{N, C}(covers1d)
end

function TileIterator(axes::Indices{N}, tilesize::Dims{N}) where {N}
    TileIterator(axes, RelaxLastTile(tilesize))::TileIterator{N}
end

"""
    titr = TileIterator(axes::NTuple{N, AbstractUnitRange}, strategy)

Decompose `axes` into an iterator `titr` of smaller axes according to `strategy`.

The `strategy` argument controls the details of the tiling. For instance
if the length of an axis is not divisible by the tile size, what should happen?
One approach would be to relax the size requirement for the last tile.
Another possibility to relax the `stride` so that all tiles are of the requested size,
but tiles may be slightly overlapping.
These two possibilities are implemented by [`RelaxLastTile`](@ref) and [`RelaxStride`](@ref).

# Examples
```jldoctest
julia> using TiledIteration

julia> collect(TileIterator((1:3, 0:5), RelaxLastTile((2, 3))))
2×2 Array{Tuple{UnitRange{Int64},UnitRange{Int64}},2}:
 (1:2, 0:2)  (1:2, 3:5)
 (3:3, 0:2)  (3:3, 3:5)

julia> collect(TileIterator((1:3, 0:5), (2, 3))) # defaults to RelaxLastTile
2×2 Array{Tuple{UnitRange{Int64},UnitRange{Int64}},2}:
 (1:2, 0:2)  (1:2, 3:5)
 (3:3, 0:2)  (3:3, 3:5)

julia> collect(TileIterator((1:3, 0:5), RelaxStride((2, 3))))
2×2 Array{Tuple{UnitRange{Int64},UnitRange{Int64}},2}:
 (1:2, 0:2)  (1:2, 3:5)
 (2:3, 0:2)  (2:3, 3:5)
```
"""
function TileIterator(axes, strategy)
    covers1d = map(cover1d, axes, split(strategy))
    return TileIterator(covers1d)
end

# strategies
"""
    RelaxStride(tilesize)

Tiling strategy, that guarantees each tile of size `tilesize`.
If the needed tiles will slightly overlap, to cover everything.

# Examples
```jldoctest
julia> using TiledIteration

julia> collect(TileIterator((1:4,), RelaxStride((2,))))
2-element Array{Tuple{UnitRange{Int64}},1}:
 (1:2,)
 (3:4,)

julia> collect(TileIterator((1:4,), RelaxStride((3,))))
2-element Array{Tuple{UnitRange{Int64}},1}:
 (1:3,)
 (2:4,)
```

See also [`TileIterator`](@ref).
"""
struct RelaxStride{N}
    tilesize::Dims{N}
end

"""
    RelaxLastTile(tilesize)

Tiling strategy, that permits the size of the last tiles along each dimension to be smaller
than `tilesize` if needed. All other tiles are of size `tilesize`.

# Examples
```jldoctest
julia> using TiledIteration

julia> collect(TileIterator((1:4,), RelaxLastTile((2,))))
2-element Array{Tuple{UnitRange{Int64}},1}:
 (1:2,)
 (3:4,)

julia> collect(TileIterator((1:7,), RelaxLastTile((2,))))
4-element Array{Tuple{UnitRange{Int64}},1}:
 (1:2,)
 (3:4,)
 (5:6,)
 (7:7,)
```

See also [`TileIterator`](@ref).
"""
struct RelaxLastTile{N}
    tilesize::Dims{N}
end

"""
    split(strategy)

Split an N dimensional strategy into an NTuple of 1 dimensional strategies.
"""
function split end

function split(strategy::RelaxStride)
    map(strategy.tilesize) do s
        RelaxStride((s,))
    end
end

function split(strategy::RelaxLastTile)
    map(strategy.tilesize) do s
        RelaxLastTile((s,))
    end
end

function cover1d(ax, strategy::RelaxStride{1})
    tilelen = stride = first(strategy.tilesize)
    firstoffset = first(ax)-1
    lastoffset = last(ax) - tilelen
    stepcount = ceil(Int, (lastoffset - firstoffset) / stride) + 1
    offsets = roundedrange(firstoffset, stop=lastoffset, length=stepcount)
    @assert last(ax) - tilelen <= last(offsets) <= last(ax)

    stopping = FixedLength(tilelen)
    return CoveredRange(offsets, stopping)
end

function cover1d(ax, strategy::RelaxLastTile{1})
    tilelen = stride = first(strategy.tilesize)
    maxstop = last(ax)
    stopping = LengthAtMost(tilelen, maxstop)

    lo = first(ax)
    hi = last(ax)
    stepcount = if tilelen <= stride
        floor(Int, (hi - lo) / stride) + 1
    else
        ceil(Int, (hi + 1 - lo - tilelen) / stride) + 1
    end
    firstoffset = lo - 1
    offsets = range(firstoffset, step=tilelen, length=stepcount)
    return CoveredRange(offsets, stopping)
end

Base.@propagate_inbounds function Base.getindex(o::TileIterator, inds::Integer...)
    cis = CartesianIndices(o)[inds...]
    o[cis]
end
Base.@propagate_inbounds function Base.getindex(o::TileIterator, ci::CartesianIndex)
    map(getindex, o.covers1d, Tuple(ci))
end

Base.size(o::TileIterator) = map(length, o.covers1d)
Base.IndexStyle(o::TileIterator) = IndexCartesian()


text:
struct NNODE{C,O} <: NeuralNetDiffEqAlgorithm
    chain::C
    opt::O
end
NNODE(chain;opt=Flux.ADAM(0.1)) = NNODE(chain,opt)

function DiffEqBase.solve(
    prob::DiffEqBase.AbstractODEProblem,
    alg::NeuralNetDiffEqAlgorithm,
    args...;
    dt,
    timeseries_errors = true,
    save_everystep=true,
    adaptive=false,
    abstol = 1f-6,
    verbose = false,
    maxiters = 100)

    DiffEqBase.isinplace(prob) && error("Only out-of-place methods are allowed!")

    u0 = prob.u0
    tspan = prob.tspan
    f = prob.f
    p = prob.p
    t0 = tspan[1]

    #hidden layer
    chain  = alg.chain
    opt    = alg.opt
    ps     = Flux.params(chain)
    data   = Iterators.repeated((), maxiters)

    #train points generation
    ts = tspan[1]:dt:tspan[2]

    #The phi trial solution
    phi(t) = u0 .+ (t .- tspan[1]).*chain(Tracker.collect([t]))

    if u0 isa Number
        dfdx = t -> Tracker.gradient(t -> sum(phi(t)), t; nest = true)[1]
        loss = () -> sum(abs2,sum(abs2,dfdx(t) .- f(phi(t)[1],p,t)[1]) for t in ts)
    else
        dfdx = t -> (phi(t+sqrt(eps(typeof(dt)))) - phi(t)) / sqrt(eps(typeof(dt)))
        #dfdx(t) = Flux.Tracker.forwarddiff(phi,t)
        #dfdx(t) = Tracker.collect([Flux.Tracker.gradient(t->phi(t)[i],t, nest=true) for i in 1:length(u0)])
        #loss function for training
        loss = () -> sum(abs2,sum(abs2,dfdx(t) - f(phi(t),p,t)) for t in ts)
    end

    cb = function ()
        l = loss()
        verbose && println("Current loss is: $l")
        l < abstol && Flux.stop()
    end
    Flux.train!(loss, ps, data, opt; cb = cb)

    #solutions at timepoints
    if u0 isa Number
        u = [phi(t)[1].data for t in ts]
    else
        u = [phi(t).data for t in ts]
    end

    sol = DiffEqBase.build_solution(prob,alg,ts,u,calculate_error = false)
    DiffEqBase.has_analytic(prob.f) && DiffEqBase.calculate_solution_errors!(sol;timeseries_errors=true,dense_errors=false)
    sol
end #solve


text:
# Autogenerated wrapper script for FastTransforms_jll for armv6l-linux-gnueabihf
export libfasttransforms

using CompilerSupportLibraries_jll
using FFTW_jll
using MPFR_jll
using OpenBLAS_jll
JLLWrappers.@generate_wrapper_header("FastTransforms")
JLLWrappers.@declare_library_product(libfasttransforms, "libfasttransforms.so")
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll, FFTW_jll, MPFR_jll, OpenBLAS_jll)
    JLLWrappers.@init_library_product(
        libfasttransforms,
        "lib/libfasttransforms.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using BlockBandedMatrices
using LinearAlgebra

function der_blocks(::Type{T}, basis::Basis, a, b) where T
    (a ∉ [0,1] || b ∉ [0,1]) &&
        error("Can only calculate derivative operators of orders 0–2!")

    g = basis.grid
    elrange = elems(g)
    n = order(g)

    f0 = zeros(elcount(g), n, n)
    for i = elrange
        f0[i,:,:] = sparse(I, n, n)
    end
    fa = a == 1 ? -basis.L′ : f0 # ∂ᴴ = -∂
    fb = b == 1 ? basis.L′ : f0

    d̃ = [zeros(T,n,n) for i = elrange]
    indices = Tuple{Integer,Integer}[]
    for i = elrange
        ii = (i-1)*(n-1) .+ 1
        push!(indices, (ii,ii))
        sel = (1:n) .+ (i-1)*(n-1)
        fael = view(fa, i, :, :)
        fbel = view(fb, i, :, :)
        d̃el = d̃[i]
        for m = 1:n
            fm = fael[m,:]
            for m′ = 1:n
                fm′ = fbel[m′,:]
                d̃el[m,m′] = dot(fm,fm′.*basis.grid.W[i,:])
            end
        end
    end

    for i = 2:elcount(g)
        d̃⁻ = d̃[i-1][end,end]
        d̃⁺ = d̃[i][1,1]
        d̃[i-1][end,end] += d̃⁺
        d̃[i][1,1] += d̃⁻
    end
    indices,d̃
end

der_blocks(::Type{T},basis::Basis,o) where T =
    der_blocks(T,basis::Basis,fld(o,2),cld(o,2))

function triderop(indices, D::AbstractVector{AbstractMatrix{T}}) where T
    tmp = Vector{T}(undef, maximum(maximum.(indices)) + size(last(D),1)-1)
    Dm = Tridiagonal(tmp[2:end], tmp, tmp[2:end])
    for (i,d) in zip(indices,D)
        s = i[1]:i[1]+size(d,1)-1
        @view(Dm[s,s]) .= d
    end
    Dm
end

function block_banded_derop(basis::Basis, D::Vector{M}) where {T,M<:AbstractMatrix{T}}
    rows = block_structure(basis)
    lu = length(D) > 1 ? vcat(1,repeat([2,1],length(D)-1)) : [0]
    m = sum(rows)
    Dm = BlockSkylineMatrix(Zeros{T}(m,m), (rows,rows), (lu,lu))
    for i in 1:length(D)
        s = 1+(i>1)
        e = size(D[i],1)-(i<length(D))
        @view(Dm[Block(2i-1,2i-1)]) .= @view((D[i])[s:e,s:e])
        if i < length(D)
            @view(Dm[Block(2i,2i)]) .= @view((D[i])[end,end])
            @view(Dm[Block(2i-1,2i)]) .= @view((D[i])[s:e,end])
            @view(Dm[Block(2i,2i-1)]) .= reshape(@view((D[i])[end,s:e]), 1, length(s:e))
        end
        if i > 1
            @view(Dm[Block(2i-1,2i-2)]) .= @view((D[i])[s:e,1])
            @view(Dm[Block(2i-2,2i-1)]) .= reshape(@view((D[i])[1,s:e]), 1, length(s:e))
        end
        if i > 1 && i < length(D)
            @view(Dm[Block(2i-2,2i)]) .= @view((D[i])[1,end])
            @view(Dm[Block(2i,2i-2)]) .= @view((D[i])[end,1])
        end
    end
    Dm
end

function derop(::Type{T}, basis::Basis, o) where T
    g = basis.grid
    n = order(g)
    elrange = elems(g)
    indices,d̃ = der_blocks(T, basis, o)
    D = [zeros(T,n,n) for i = elrange]
    for i = elrange
        Del = D[i]
        d̃el = d̃[i]
        for m = 1:n
            for m′ = 1:n
                Del[m,m′] = g.N[i,m]*g.N[i,m′]*d̃el[m,m′]
            end
        end
    end

    if g.bl == :dirichlet0
        D[1] = D[1][2:end,2:end]
        for i in 2:length(indices)
            indices[i] = (indices[i][1]-1,indices[i][2]-1)
        end
    end
    if g.br == :dirichlet0
        D[end] = D[end][1:end-1,1:end-1]
    end

    n == 2 && return triderop(indices, D)

    block_banded_derop(basis, D)
end

derop(basis::Basis, o) =
    derop(eltype(basis), basis, o)

function kinop(::Type{T}, basis::Basis) where T
    D2 = derop(T, basis, 2)
    D2 ./= -2
    D2
end
kinop(basis::Basis) = kinop(eltype(basis), basis)

export derop, kinop