text:
#===============================================================================

    verify_response.jl

    Author: Tom Price
    Date: Dec 2018

    This Julia script tests the JSON data structures
    that are returned in the body of responses to GET requests

    *     based on documentation of REST API from juliaapi_new.txt

    *     Julia listens on http://localhost:8081 via HttpServer.jl and dispatch.jl
          and sends responses to the Rails app on http://localhost:3000

    *     the relevant controller in the Rails app is
          chaipcr/web/app/controllers/experiments_controller.rb

    *     currently calls are GET-ed to
          http://127.0.0.1:8081/experiments/#{experiment.id}/standard_curve
          http://127.0.0.1:8081/experiments/#{experiment.id}/amplification
          http://127.0.0.1:8081/experiments/#{experiment.id}/meltcurve
          http://127.0.0.1:8081/experiments/#{experiment.id}/analyze

===============================================================================#

import JSON, DataStructures.OrderedDict

#===============================================================================

    call: experiments/:experiment_id/standard_curve

===============================================================================#

function verify_response(
    ::Type{Val{standard_curve}},
    response ::Any
)
    facts() do
        context("Verifying response body") do
            @fact (isa(response,OrderedDict)) --> true
            @fact (haskey(response,"valid")) --> true
            if (haskey(response,"error"))
                @fact (length(response)) --> 2
                @fact (isa(response["error"],String)) --> true
                @fact (response["valid"]) --> false
            else
                @fact (response["valid"]) --> true
                @fact (haskey(response,"targets")) --> true
                # @fact (haskey(response,"groups")) --> true # do not verify groups
                array=response["targets"]
                @fact (isa(array,Vector)) --> true
                for i in range(1,length(array))
                    dict=array[i]
                    @fact (isa(dict,OrderedDict)) --> true
                    @fact (haskey(dict,"target_id")) --> true
                    if (length(dict)==2)
                        @fact (haskey(dict,"error")) --> true
                        @fact (isa(dict["error"],String)) --> true
                    else
                        @fact (length(dict)) --> 5
                        @fact (haskey(dict,"slope")) --> true
                        @fact (isa(dict["slope"],Number)) --> true
                        @fact (haskey(dict,"offset")) --> true
                        @fact (isa(dict["offset"],Number)) --> true
                        @fact (haskey(dict,"efficiency")) --> true
                        @fact (isa(dict["efficiency"],Number)) --> true
                        @fact (haskey(dict,"r2")) --> true
                        @fact (isa(dict["r2"],Number)) --> true
                    end
                end
            end
        end
    end
    FactCheck.exitstatus()
end

#===============================================================================

    call: experiments/:experiment_id/amplification

===============================================================================#

function verify_response(
    ::Type{Val{amplification}},
    response ::Any
)
    facts() do
        context("Verifying response body") do
            @fact (isa(response,OrderedDict)) --> true
            @fact (haskey(response,"valid")) --> true
            if (haskey(response,"error"))
                @fact (length(response)) --> 2
                @fact (isa(response["error"],String)) --> true
                @fact (response["valid"]) --> false
            else
                # @fact (length(response)) --> 8
                @fact (response["valid"]) --> true
                measurements=["rbbs_ary3","blsub_fluos","dr1_pred","dr2_pred"]
                n_channels=length(response["rbbs_ary3"])
                @fact (n_channels in QpcrAnalysis.CHANNELS) --> true
                n_wells=length(response["rbbs_ary3"][1])
                n_steps=length(response["rbbs_ary3"][1][1])
                n_pred=length(response["dr1_pred"][1][1])
                for m in measurements
                    @fact (haskey(response,m)) --> true
                    @fact (isa(response[m],Vector)) --> true
                    @fact (length(response[m])) --> n_channels
                    for c in range(1,n_channels)
                        @fact (isa(response[m][c],Vector)) --> true
                        @fact (length(response[m][c])) --> n_wells
                        for i in range(1,n_wells)
                            @fact (isa(response[m][c][i],Vector)) --> true
                            if (m=="rbbs_ary3" || m=="blsub_fluos")
                                @fact (length(response[m][c][i])) --> n_steps
                                for j in range(1,n_steps)
                                    @fact (isa(response[m][c][i][j],Number) ||
                                        (response[m][c][i][j]==nothing)) --> true
                                end
                            else # dr1_pred, dr2_pred
                                @fact (length(response[m][c][i])) --> n_pred
                                for j in range(1,n_pred)
                                    @fact (isa(response[m][c][i][j],Number) ||
                                        (response[m][c][i][j]==nothing)) --> true
                                end
                            end
                        end
                    end
                end
                statistics=["cq","d0"]
                for s in statistics
                    @fact (haskey(response,s)) --> true
                    @fact (isa(response[s],Vector)) --> true
                    @fact (length(response[s])) --> n_channels
                    for c in range(1,n_channels)
                        @fact (isa(response[s][c],Vector)) --> true
                        @fact (length(response[s][c])) --> n_wells
                        for i in range(1,n_wells)
                            @fact (isa(response[s][c][i],Number) ||
                                response[s][c][i]==nothing) --> true
                        end
                    end
                end
                @fact (haskey(response,"ct_fluos")) --> true
                @fact (isa(response["ct_fluos"],Vector)) --> true
                @fact (length(response["ct_fluos"])) --> n_channels
                for c in range(1,n_channels)
                    @fact (isa(response["ct_fluos"][c],Number) ||
                        response["ct_fluos"][c]==nothing) --> true
                end
                variables=["rbbs_ary3","blsub_fluos","cq","d0"]
                @fact (haskey(response,"assignments_adj_labels_dict")) --> true
                @fact (isa(response["assignments_adj_labels_dict"],OrderedDict)) --> true
                # @fact (length(response["assignments_adj_labels_dict"])) --> n_genotypes
                for g in range(1,length(response["assignments_adj_labels_dict"]))
                    @fact (isa(response["assignments_adj_labels_dict"][variables[g]],Vector)) --> true
                    @fact (length(response["assignments_adj_labels_dict"][variables[g]])) --> n_wells
                    for i in range(1,n_wells)
                        @fact (isa(response["assignments_adj_labels_dict"][variables[g]][i],String)) --> true
                    end
                end
            end
        end
    end
    FactCheck.exitstatus()
end

#===============================================================================

    call: experiments/:experiment_id/meltcurve

===============================================================================#

function verify_response(
    ::Type{Val{meltcurve}},
    response ::Any
)
    facts() do
        context("Verifying response body") do
            @fact (isa(response,OrderedDict)) --> true
            @fact (haskey(response,"valid")) --> true
            if (haskey(response,"error"))
                @fact (length(response)) --> 2
                @fact (isa(response["error"],String)) --> true
                @fact (response["valid"]) --> false
            else
                @fact (response["valid"]) --> true
                variables=["melt_curve_data","melt_curve_analysis"]
                # @fact length(response) --> length(variables)
                n_channels=length(response["melt_curve_data"])
                n_wells=length(response["melt_curve_data"][1])
                n_grid=length(response["melt_curve_data"][1][1][1])
                for v in variables
                    n = (v=="melt_curve_data") ? 3 : 2
                    @fact (haskey(response,v)) --> true
                    @fact (isa(response[v],Vector)) --> true
                    @fact (length(response[v])) --> n_channels
                    for i in range(1,n_channels) # channel
                        @fact (isa(response[v][i],Vector)) --> true
                        for j in range(1,n_wells) # well
                            @fact (length(response[v][i][j])) --> n
                            for k in range(1,n) # temperature, fluorescence, slope / Tm, area
                                if (v=="melt_curve_data")
                                    @fact abs(length(response[v][i][j][k]) - n_grid) --> less_than_or_equal(1)
                                end
                                @fact (isa(response[v][i][j][k],Vector)) --> true
                                for m in range(1,length(response[v][i][j][k])) # prediction locations / temp maxima
                                    @fact (isa(response[v][i][j][k][m],Number)) --> true
                                end
                            end
                        end
                    end
                end
            end
        end
    end
    FactCheck.exitstatus()
end

#===============================================================================

    call: system/loadscript?script=path%2Fto%2Fanalyze.jl

===============================================================================#

function verify_response(
    ::Type{Val{loadscript}},
    response ::Any
)
    facts() do
        context("Verifying response body") do
            @fact (isa(response,OrderedDict)) --> true
            @fact (haskey(response,"valid")) --> true
            if (haskey(response,"error"))
                @fact (length(response)) --> 2
                @fact (isa(response["error"],String)) --> true
                @fact (response["valid"]) --> false
            else
                @fact (response["valid"]) --> true
                @fact (haskey(response,"script")) --> true
                @fact (isa(response["script"],String)) --> true
            end
        end
    end
    FactCheck.exitstatus()
end

#===============================================================================

    call: experiments/:experiment_id/thermal_performance_diagnostic

===============================================================================#

function verify_response(
    ::Type{Val{thermal_performance_diagnostic}},
    response ::Any
)
    facts() do
        context("Verifying response body") do
            @fact (isa(response,OrderedDict)) --> true
            @fact (haskey(response,"valid")) --> true
            if (haskey(response,"error"))
                @fact (length(response)) --> 2
                @fact (isa(response["error"],String)) --> true
                @fact (response["valid"]) --> false
            else
                @fact (length(response)) --> 4
                @fact (response["valid"]) --> true
                @fact (haskey(response,"Heating")) --> true
                @fact (haskey(response,"Cooling")) --> true
                @fact (haskey(response,"Lid")) --> true
                @fact (isa(response["Heating"],OrderedDict)) --> true
                @fact (isa(response["Cooling"],OrderedDict)) --> true
                @fact (isa(response["Lid"],OrderedDict)) --> true
                @fact (length(response["Heating"])) --> 3
                @fact (length(response["Cooling"])) --> 3
                @fact (length(response["Lid"])) --> 2
                @fact (haskey(response["Heating"],"AvgRampRate")) --> true
                @fact (haskey(response["Cooling"],"AvgRampRate")) --> true
                @fact (haskey(response["Lid"],"HeatingRate")) --> true
                @fact (haskey(response["Heating"],"TotalTime")) --> true
                @fact (haskey(response["Cooling"],"TotalTime")) --> true
                @fact (haskey(response["Lid"],"TotalTime")) --> true
                @fact (haskey(response["Heating"],"MaxBlockDeltaT")) --> true
                @fact (haskey(response["Cooling"],"MaxBlockDeltaT")) --> true
                @fact (isa(response["Heating"]["AvgRampRate"],Vector)) --> true
                @fact (isa(response["Cooling"]["AvgRampRate"],Vector)) --> true
                @fact (isa(response["Lid"]["HeatingRate"],Vector)) --> true
                @fact (isa(response["Heating"]["TotalTime"],Vector)) --> true
                @fact (isa(response["Cooling"]["TotalTime"],Vector)) --> true
                @fact (isa(response["Lid"]["TotalTime"],Vector)) --> true
                @fact (isa(response["Heating"]["MaxBlockDeltaT"],Vector)) --> true
                @fact (isa(response["Cooling"]["MaxBlockDeltaT"],Vector)) --> true
                @fact (length(response["Heating"]["AvgRampRate"])) --> 2
                @fact (length(response["Cooling"]["AvgRampRate"])) --> 2
                @fact (length(response["Lid"]["HeatingRate"])) --> 2
                @fact (length(response["Heating"]["TotalTime"])) --> 2
                @fact (length(response["Cooling"]["TotalTime"])) --> 2
                @fact (length(response["Lid"]["TotalTime"])) --> 2
                @fact (length(response["Heating"]["MaxBlockDeltaT"])) --> 2
                @fact (length(response["Cooling"]["MaxBlockDeltaT"])) --> 2
                @fact (isa(response["Heating"]["AvgRampRate"][1],Number)) --> true
                @fact (isa(response["Cooling"]["AvgRampRate"][1],Number)) --> true
                @fact (isa(response["Lid"]["HeatingRate"][1],Number)) --> true
                @fact (isa(response["Heating"]["TotalTime"][1],Number)) --> true
                @fact (isa(response["Cooling"]["TotalTime"][1],Number)) --> true
                @fact (isa(response["Lid"]["TotalTime"][1],Number)) --> true
                @fact (isa(response["Heating"]["MaxBlockDeltaT"][1],Number)) --> true
                @fact (isa(response["Cooling"]["MaxBlockDeltaT"][1],Number)) --> true
                @fact (isa(response["Heating"]["AvgRampRate"][2],Bool)) --> true
                @fact (isa(response["Cooling"]["AvgRampRate"][2],Bool)) --> true
                @fact (isa(response["Lid"]["HeatingRate"][2],Bool)) --> true
                @fact (isa(response["Heating"]["TotalTime"][2],Bool)) --> true
                @fact (isa(response["Cooling"]["TotalTime"][2],Bool)) --> true
                @fact (isa(response["Lid"]["TotalTime"][2],Bool)) --> true
                @fact (isa(response["Heating"]["MaxBlockDeltaT"][2],Bool)) --> true
                @fact (isa(response["Cooling"]["MaxBlockDeltaT"][2],Bool)) --> true
            end
        end
    end
    FactCheck.exitstatus()
end

#===============================================================================

    call: experiments/:experiment_id/thermal_consistency

===============================================================================#

function verify_response(
    ::Type{Val{thermal_consistency}},
    response ::Any
)
    facts() do
        context("Verifying response body") do
            @fact (isa(response,OrderedDict)) --> true
            @fact (haskey(response,"valid")) --> true
            if (haskey(response,"error"))
                @fact (length(response)) --> 2
                @fact (isa(response["error"],String)) --> true
                @fact (response["valid"]) --> false
            else
                @fact (length(response)) --> 3
                @fact (response["valid"]) --> true
                @fact (haskey(response,"tm_check")) --> true
                @fact (haskey(response,"delta_Tm")) --> true
                @fact (isa(response["tm_check"],Vector)) --> true
                for i in range(1,length(response["tm_check"]))
                    @fact (isa(response["tm_check"][i],OrderedDict)) --> true
                    @fact (length(response["tm_check"][i])) --> 2
                    @fact (haskey(response["tm_check"][i],"Tm")) --> true
                    @fact (haskey(response["tm_check"][i],"area")) --> true
                    @fact (isa(response["tm_check"][i]["Tm"],Vector)) --> true
                    @fact (length(response["tm_check"][i]["Tm"])) --> 2
                    @fact (isa(response["tm_check"][i]["Tm"][1],Number)) --> true
                    @fact (isa(response["tm_check"][i]["Tm"][2],Bool)) --> true
                    @fact (isa(response["tm_check"][i]["area"],Number)) --> true
                end
                @fact (isa(response["delta_Tm"],Vector)) --> true
                @fact (length(response["delta_Tm"])) --> 2
                @fact (isa(response["delta_Tm"][1],Number)) --> true
                @fact (isa(response["delta_Tm"][2],Bool)) --> true
            end
        end
    end
    FactCheck.exitstatus()
end

#===============================================================================

    call: experiments/:experiment_id/optical_cal

===============================================================================#

# success response body (optical_cal):

function verify_response(
    ::Type{Val{optical_cal}},
    response ::Any
)
    facts() do
        context("Verifying response body") do
            @fact (isa(response,OrderedDict)) --> true
            if (response["valid"]==true)
                @fact (length(response)) --> 1
            else
                @fact (response["valid"]) --> false
                @fact (length(response)) --> 2
                @fact (haskey(response,"error")) --> true
                @fact (isa(response["error"],String)) --> true
            end
        end
    end
    FactCheck.exitstatus()
end

#===============================================================================

    call: experiments/:experiment_id/optical_test_single_channel

===============================================================================#

function verify_response(
    ::Type{Val{optical_test_single_channel}},
    response ::Any
)
    facts() do
        context("Verifying response body") do
            @fact (isa(response,OrderedDict)) --> true
            @fact (haskey(response,"valid")) --> true
            if (haskey(response,"error"))
                @fact (length(response)) --> 2
                @fact (isa(response["error"],String)) --> true
                @fact (response["valid"]) --> false
            else
                @fact (response["valid"]) --> true
                @fact (haskey(response,"optical_data")) --> true
                @fact (isa(response["optical_data"],Vector)) --> true
                for i in range(1,length(response["optical_data"]))
                    @fact (isa(response["optical_data"][i],OrderedDict)) --> true
                    @fact (length(response["optical_data"][i])) --> 3
                    @fact (haskey(response["optical_data"][i],"baseline")) --> true
                    @fact (haskey(response["optical_data"][i],"excitation")) --> true
                    @fact (haskey(response["optical_data"][i],"valid")) --> true
                    @fact (isa(response["optical_data"][i]["baseline"],Number)) --> true
                    @fact (isa(response["optical_data"][i]["excitation"],Number)) --> true
                    @fact (isa(response["optical_data"][i]["valid"],Bool)) --> true
                end
            end
        end
    end
    FactCheck.exitstatus()
end

#===============================================================================

    call: experiments/:experiment_id/optical_test_dual_channel

===============================================================================#

function verify_response(
    ::Type{Val{optical_test_dual_channel}},
    response ::Any
)
    facts() do
        context("Verifying response body") do
            @fact (isa(response,OrderedDict)) --> true
            @fact (haskey(response,"valid")) --> true
            if (haskey(response,"error"))
                @fact (length(response)) --> 2
                @fact (isa(response["error"],String)) --> true
                @fact (response["valid"]) --> false
            else
                signals=["baseline","water","HEX","FAM"]
                @fact (response["valid"]) --> true
                @fact (length(response)) --> 3
                @fact (haskey(response,"optical_data")) --> true
                @fact (isa(response["optical_data"],Vector)) --> true
                n_wells=length(response["optical_data"])
                for i in range(1,n_wells) # well
                    @fact (isa(response["optical_data"][i],OrderedDict)) --> true
                    @fact (length(response["optical_data"][i])) --> length(signals)
                    for signal in signals
                        @fact (haskey(response["optical_data"][i],signal)) --> true
                        @fact (isa(response["optical_data"][i][signal],Vector)) --> true
                        @fact (length(response["optical_data"][i][signal])) --> 2
                        for j in range(1,2) # channel
                            @fact (isa(response["optical_data"][i][signal][j],Vector)) --> true
                            @fact (length(response["optical_data"][i][signal][j])) --> 2
                            @fact (isa(response["optical_data"][i][signal][j][1],Number)) --> true
                            @fact (isa(response["optical_data"][i][signal][j][2],Bool)) --> true
                        end
                    end
                end
                @fact (haskey(response,"Ch1:Ch2")) --> true
                @fact (isa(response["Ch1:Ch2"],OrderedDict)) --> true
                @fact (length(response["Ch1:Ch2"])) --> 2
                for signal in signals[3:4]
                    @fact (haskey(response["Ch1:Ch2"],signal)) --> true
                    @fact (isa(response["Ch1:Ch2"][signal],Vector)) --> true
                    @fact (length(response["Ch1:Ch2"][signal])) --> n_wells
                    for i in range(1,n_wells) # well
                        @fact (isa(response["Ch1:Ch2"][signal][i],Number) ||
                            response["Ch1:Ch2"][signal][i] == nothing) --> true # division by zero: JSON parser converts Inf, -Inf, NaN -> "null"
                    end
                end
            end
        end
    end
    FactCheck.exitstatus()
end

max_stars_repo_name: MakerButt/chaipcr

max_issues_repo_name: MakerButt/chaipcr

max_forks_repo_name: MakerButt/chaipcr


text:
# Functions for sampling from populations

### draw a pair of distinct integers in [1:n]

function samplepair(n::Int)
    i1 = randi(n)
    i2 = randi(n-1)
    return (i1, i2 == i1 ? n : i2)
end

function samplepair(a::AbstractArray)
    i1, i2 = samplepair(length(a))
    return a[i1], a[i2]
end

### internal sampling algorithms

# Fisher-Yates sampling
immutable FisherYatesSampler
    n::Int
    seq::Vector{Int}   # Internal sequence for shuffling

    FisherYatesSampler(n::Int) = new(n, [1:n])
end

function rand!(s::FisherYatesSampler, a::AbstractArray, x::AbstractArray)
    # draw samples without-replacement to x

    n::Int = s.n
    k::Int = length(x)
    if k > n
        throw(ArgumentError("Cannot draw more than n samples without replacement."))
    end

    seq::Vector{Int} = s.seq
    for i = 1:k
        j = randi(i, n)
        sj = seq[j]
        x[i] = a[sj]
        seq[j] = seq[i]
        seq[i] = sj
    end
    x
end

fisher_yates_sample!(a::AbstractArray, x::AbstractArray) = rand!(FisherYatesSampler(length(a)), a, x)

# self-avoiding sampling
function self_avoid_sample!{T}(a::AbstractArray{T}, x::AbstractArray)
    # This algorithm is suitable when length(x) << length(a)

    s = Set{T}()
    # sizehint(s, length(x))
    rgen = RandIntSampler(length(a))

    # first one
    idx = rand(rgen)
    x[1] = a[idx]
    push!(s, idx)

    # remaining
    for i = 2:length(x)
        idx = rand(rgen)
        while (idx in s)
            idx = rand(rgen)
        end
        x[i] = a[idx]
        push!(s, idx)
    end
    x
end

# Ordered sampling without replacement
# Original author: Mike Innes

function ordered_sample!(a::AbstractArray, x::AbstractArray)
    n = length(a)
    k = length(x)
    offset = 0
    i = 1

    while offset < k
        rk = k - offset
        if i == n
            for j = 1 : rk
                @inbounds x[offset + j] = a[i]
            end
            offset = k
        else
            m = rand_binom(rk, 1.0 / (n - i + 1))
            for j = 1 : m
                @inbounds x[offset + j] = a[i]
            end
            i += 1
            offset += m
        end
    end
    x
end

###########################################################
#
#   Interface functions
#
###########################################################

sample(a::AbstractArray) = a[randi(length(a))]

function sample!(a::AbstractArray, x::AbstractArray; replace::Bool=true, ordered::Bool=false)
    n = length(a)
    k = length(x)

    if isempty(x)
        return x
    end

    if replace
        if ordered
            ordered_sample!(a, x)
        else
            s = RandIntSampler(n)
            for i = 1:k
                @inbounds x[i] = a[rand(s)]
            end
        end

    else  # without replacement
        k <= n || error("Cannot draw more samples without replacement.")

        if k == 1
            @inbounds x[1] = sample(a)

        elseif k == 2
            @inbounds (x[1], x[2]) = samplepair(a)

        elseif n < k * max(k, 100)
            fisher_yates_sample!(a, x)

        else
            self_avoid_sample!(a, x)
        end

        if ordered
            sort!(x)
        end
    end
    return x
end

function sample{T}(a::AbstractArray{T}, n::Integer; replace::Bool=true, ordered::Bool=false)
    sample!(a, Array(T, n); replace=replace, ordered=ordered)
end

function sample{T}(a::AbstractArray{T}, dims::Dims; replace::Bool=true, ordered::Bool=false)
    sample!(a, Array(T, dims); replace=replace, ordered=ordered)
end

################################################################
#
#  Weighted sampling
#
################################################################

function sample(wv::WeightVec)
    t = rand() * sum(wv)
    w = values(wv)
    n = length(w)
    i = 1
    cw = w[1]
    while cw < t && i < n
        i += 1
        @inbounds cw += w[i]
    end
    return i
end

sample(a::AbstractArray, wv::WeightVec) = a[sample(wv)]

# Original author: Mike Innes
function ordered_sample!(a::AbstractArray, wv::WeightVec, x::AbstractArray)
    n = length(a)
    k = length(x)
    offset = 0
    i = 1
    wsum = sum(wv)
    w = values(wv)

    while offset < k
        rk = k - offset
        wi = w[i]

        if i == n || wi >= wsum
            for j = 1 : rk
                @inbounds x[offset + j] = a[i]
            end
            offset = k
        else
            m = rand_binom(rk, wi / wsum)
            for j = 1 : m
                @inbounds x[offset + j] = a[i]
            end
            i += 1
            wsum -= wi
            offset += m
        end
    end
    x
end

function sample!(a::AbstractArray, wv::WeightVec, x::AbstractArray; ordered::Bool=false)
    n = length(a)
    k = length(x)

    if ordered
        ordered_sample!(a, wv, x)
    else
        if k > 100 * n
            shuffle!(ordered_sample!(a, wv, x))
        else
            for i = 1 : k
                @inbounds x[i] = sample(a, wv)
            end
        end
    end
    return x
end

sample{T}(a::AbstractArray{T}, wv::WeightVec, n::Integer; ordered::Bool=false) =
    sample!(a, wv, Array(T, n); ordered=ordered)

sample{T}(a::AbstractArray{T}, wv::WeightVec, dims::Dims; ordered::Bool=false) =
    sample!(a, wv, Array(T, dims); ordered=ordered)

# wsample interface

wsample(w::RealVector) = sample(weights(w))
wsample(a::AbstractArray, w::RealVector) = sample(a, weights(w))

wsample!(a::AbstractArray, w::RealVector, x::AbstractArray; ordered::Bool=false) =
    sample!(a, weights(w), x; ordered=ordered)

wsample{T}(a::AbstractArray{T}, w::RealVector, n::Integer; ordered::Bool=false) =
    wsample!(a, w, Array(T, n); ordered=ordered)

wsample{T}(a::AbstractArray{T}, w::RealVector, dims::Dims; ordered::Bool=false) =
    wsample!(a, w, Array(T, dims); ordered=ordered)

ext: jl


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

mutable struct IoK8sApiExtensionsV1beta1Scale <: SwaggerModel
    apiVersion::Any # spec type: Union{ Nothing, String } # spec name: apiVersion
    kind::Any # spec type: Union{ Nothing, String } # spec name: kind
    metadata::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1ObjectMeta } # spec name: metadata
    spec::Any # spec type: Union{ Nothing, IoK8sApiExtensionsV1beta1ScaleSpec } # spec name: spec
    status::Any # spec type: Union{ Nothing, IoK8sApiExtensionsV1beta1ScaleStatus } # spec name: status

    function IoK8sApiExtensionsV1beta1Scale(;apiVersion=nothing, kind=nothing, metadata=nothing, spec=nothing, status=nothing)
        o = new()
        validate_property(IoK8sApiExtensionsV1beta1Scale, Symbol("apiVersion"), apiVersion)
        setfield!(o, Symbol("apiVersion"), apiVersion)
        validate_property(IoK8sApiExtensionsV1beta1Scale, Symbol("kind"), kind)
        setfield!(o, Symbol("kind"), kind)
        validate_property(IoK8sApiExtensionsV1beta1Scale, Symbol("metadata"), metadata)
        setfield!(o, Symbol("metadata"), metadata)
        validate_property(IoK8sApiExtensionsV1beta1Scale, Symbol("spec"), spec)
        setfield!(o, Symbol("spec"), spec)
        validate_property(IoK8sApiExtensionsV1beta1Scale, Symbol("status"), status)
        setfield!(o, Symbol("status"), status)
        o
    end
end # type IoK8sApiExtensionsV1beta1Scale

const _property_map_IoK8sApiExtensionsV1beta1Scale = Dict{Symbol,Symbol}(Symbol("apiVersion")=>Symbol("apiVersion"), Symbol("kind")=>Symbol("kind"), Symbol("metadata")=>Symbol("metadata"), Symbol("spec")=>Symbol("spec"), Symbol("status")=>Symbol("status"))
const _property_types_IoK8sApiExtensionsV1beta1Scale = Dict{Symbol,String}(Symbol("apiVersion")=>"String", Symbol("kind")=>"String", Symbol("metadata")=>"IoK8sApimachineryPkgApisMetaV1ObjectMeta", Symbol("spec")=>"IoK8sApiExtensionsV1beta1ScaleSpec", Symbol("status")=>"IoK8sApiExtensionsV1beta1ScaleStatus")
Base.propertynames(::Type{ IoK8sApiExtensionsV1beta1Scale }) = collect(keys(_property_map_IoK8sApiExtensionsV1beta1Scale))
Swagger.property_type(::Type{ IoK8sApiExtensionsV1beta1Scale }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiExtensionsV1beta1Scale[name]))}
Swagger.field_name(::Type{ IoK8sApiExtensionsV1beta1Scale }, property_name::Symbol) =  _property_map_IoK8sApiExtensionsV1beta1Scale[property_name]

function check_required(o::IoK8sApiExtensionsV1beta1Scale)
    true
end

function validate_property(::Type{ IoK8sApiExtensionsV1beta1Scale }, name::Symbol, val)
end


text:
using Distributions
import CovSel
import HD

# simple model that differes in only one edge
srand(123)

p = 20
Sigmax = eye(p,p)
Sigmay = zeros(p,p)
rho = 0.7
for i=1:p
    for j=1:p
        Sigmay[i,j]=rho^abs(i-j)
    end
end
sqmy = sqrtm(Sigmay)
n = 1000
X = randn(n,p)
Y = randn(n,p) * sqmy;
hSx = cov(X)
hSy = cov(Y)

tmp = ones(p,p)

# λ = 0.2
λ = rand(Uniform(0.15, 0.2))
λvec = λ*ones(div(p*(p+1), 2))

reload("HD")
reload("CovSel")

@time o1 = CovSel.Alt.differencePrecisionNaive(hSx, hSy, λ, tmp)
@time o2 = CovSel.Alt.differencePrecisionActiveShooting(hSx, hSy, λ, tmp)
@time o3 = CovSel.differencePrecisionActiveShooting(hSx, hSy, λvec)

@show maximum(abs.(o1 - o2))
@show maximum(abs.(CovSel.vec2tril(o3, p) - tril(o1)))

@show trace(hSx * o1 * hSy * o1) / 2 - trace(o1*(hSy-hSx))
o3 = CovSel.tril2symmetric(CovSel.vec2tril(o3, p))
@show trace(hSx * o3 * hSy * o3) / 2 - trace(o3*(hSy-hSx))

# non_zero_set = find( abs(precM) .> 1e-4 )
# @time CovSel.covsel_refit!(X, Z, U, S, non_zero_set)
#
# # find connected components
# import Images
# nonZero = zeros(S)
# for i in eachindex(S)
#   if abs(S[i]) > λ
#     nonZero[ i ] = 1
#   end
# end
# labels = Images.label_components(nonZero)
#
# function collect_groups(labels)
#     groups = [Int[] for i = 1:maximum(labels)]
#     for (i,l) in enumerate(labels)
#         if l != 0
#             push!(groups[l], i)
#         end
#     end
#     groups
# end
#
# groups = collect_groups(labels)
# S[ groups[1] ]
#
# S[1:2, 1:2]


text:
# Triangular distribution
immutable TriangularDist <: ContinuousUnivariateDistribution
    a::Float64
    b::Float64
    c::Float64
    function TriangularDist(a::Real, b::Real, c::Real)
        a < b || error("a<b must be true")
        a <= c <= b || error("a<=c<=b must be true")
        new(float64(a), float64(b), float64(c))
    end
end

## Support
isupperbounded(::Union(TriangularDist, Type{TriangularDist})) = true
islowerbounded(::Union(TriangularDist, Type{TriangularDist})) = true
isbounded(::Union(TriangularDist, Type{TriangularDist})) = true

minimum(d::TriangularDist) = d.a
maximum(d::TriangularDist) = d.b
insupport(d::TriangularDist, x::Real) = minimum(d) <= x <= maximum(d)

## Properties
mean(d::TriangularDist) = (d.a + d.b + d.c) / 3.0
median(d::TriangularDist) = d.c >= (d.a+d.b)/2.0 ?
    d.a + sqrt((d.b-d.a)*(d.c-d.a))/sqrt(2.0) :
    d.b - sqrt((d.b-d.a)*(d.b-d.c))/sqrt(2.0)
mode(d::TriangularDist) = d.c

var(d::TriangularDist) = (d.a^2.0 + d.b^2.0 + d.c^2.0-d.a*d.b-d.a*d.c-d.b*d.c)/18.0
skewness(d::TriangularDist) = sqrt(2.0)*(d.a+d.b-2.0d.c)*(2d.a-d.b-d.c)*(d.a-2d.b+d.c)/5.0/(d.a^2.0 + d.b^2.0 + d.c^2.0-d.a*d.b-d.a*d.c-d.b*d.c)^(3.0/2.0)
kurtosis(d::TriangularDist) = -0.6

entropy(d::TriangularDist) = 0.5 + log((d.b-d.a)/2.0)

## Functions
function pdf(d::TriangularDist, x::Real)
    if d.a<=x<=d.c
        return 2.0*(x-d.a)/(d.b-d.a)/(d.c-d.a)
    elseif d.c<x<=d.b
        return 2.0*(d.b-x)/(d.b-d.a)/(d.b-d.c)
    else
        return zero(x)
    end
end

function cdf(d::TriangularDist, x::Real)
    if x<d.a
        return zero(x)
    elseif d.a<=x<=d.c
        return (x-d.a)^2/(d.b-d.a)/(d.c-d.a)
    elseif d.c<x<=d.b
        return 1.0-(d.b-x)^2/(d.b-d.a)/(d.b-d.c)
    else
        return one(x)
    end
end

function quantile(d::TriangularDist, p::Real)
    @checkquantile p begin
        if p <= (d.c-d.a)/(d.b-d.a)
            return d.a + sqrt((d.b-d.a)*(d.c-d.a)*p)
        else
            return d.b - sqrt((d.b-d.a)*(d.b-d.c)*(1-p))
        end
    end
end

function mgf(d::TriangularDist, t::Real)
    if t==zero(t)
        return one(t)
    else
        nominator = (d.b-d.c)*exp(d.a*t)-(d.b-d.a)*exp(d.c*t)+(d.c-d.a)*exp(d.b*t)
        denominator = (d.b-d.a)*(d.c-d.a)*(d.b-d.c)*t^2
        return 2*nominator/denominator
    end
end

function cf(d::TriangularDist, t::Real)
    # Is this correct?
    if t==zero(t)
        return one(t)
    else
        nominator = (d.b-d.c)*exp(im*d.a*t)-(d.b-d.a)*exp(im*d.c*t)+(d.c-d.a)*exp(im*d.b*t)
        denominator = (d.b-d.a)*(d.c-d.a)*(d.b-d.c)*t^2
        return -2*nominator/denominator
    end
end

## Sampling
function rand(d::TriangularDist)
    u = rand()

    if u < (d.c-d.a)/(d.b-d.a)
        return d.a + sqrt(u*(d.b-d.a)*(d.c-d.a))
    else
        return d.b - sqrt((1.0-u)*(d.b-d.a)*(d.b-d.c))
    end
end


text:
"""
reorder!(X,V)
Reorder Array `X` according to its correlation with `V`. Uses Hungarian
Algo to maximise diagonal of correlation matrix.
"""
function reorder!(X::AbstractArray{T,2},V::AbstractArray{T,2}) where T
    cormat = cor(V, X)
    sur_perm = munkres(-cormat)
    X .= X[:,sur_perm]
end

"""
function normweightvectors!(V)
Normalise the weight vector of each pattern so that the maximum absolute value
is one.
"""
function normweightvectors!(V::AbstractArray{T,2}) where T
    for i in 1:size(V,2)
        weightvec = view(V,:,i)
        # flip if two of three largest absolute weights are negative
        if mean(sign.(weightvec[findnmax(abs.(weightvec),3)])) < 0
            weightvec .= -weightvec
        end
        weightvec ./= norm(weightvec)
    end
    V .= V[:,sortperm(abs.(vec(sum(V,dims=1))))]
end

"""
function track(Z, V, M)
Track zscored firing rates Z using weight vectors V.
Tracking uses quadratic form:
R(t) = z(t)' * Pk(t) * (zt)
where Pk is masked by matrix M.
"""
function track(Z::AbstractArray{T,2}, V::AbstractArray{T,2},
               M::AbstractArray{Bool,2}) where T<:Float64
    K = size(V,2)
    R = zeros(Float64, size(Z,2), K)
    for k in 1:K
        Rk = view(R,:,k)
        Pk = V[:,k]*V[:,k]'
        Pk[M] .= 0.
        trackK!(Rk, Pk, Z)
    end
    return R
end

function trackK!(Rk,Pk,Z)
    for t in 1:size(Z,2)
        Rk[t] = Z[:,t]'*Pk*Z[:,t]
    end
    nothing
end

function strength(Z::AbstractArray{T,2}, V::AbstractArray{T,2}) where T<:Float64
    K = size(V,2)
    R = zeros(Float64, K)
    for k in 1:K
        Pk = V[:,k]*V[:,k]'
        Pk[diagind(Pk)] .= 0.
        for t in 1:size(Z,2)
            R[k] += Z[:,t]'*Pk*Z[:,t]
        end
    end
    return R
end

"""
track_partial(Z::AbstractArray{Float64,2}, V::AbstractArray{Float64,2}, multipliers)
As `track` but with mask on the outer product of `V` (`Pk`). Useful for assessing contribution of specific sets of neurons,
or specific interaction between sets of neurons.
"""
function track_partial(Z::AbstractArray{Float64,2}, V::AbstractArray{Float64,2}, multipliers)
    npatterns = size(V, 2)
    R = zeros(size(Z,1), npatterns, length(multipliers))
    Z = Z'
    for k in 1:npatterns
        for (mi,m) in enumerate(multipliers)
            Pk = V[:,k]*V[:,k]'
            Pk[diagind(Pk)] .= 0.;
            Pk .*= m
            trackK!(view(R,:,k, mi), Pk, Z)
        end
    end
    return R
end

"""
marchenko_thresh(n,B)
Get threshold for eigenvalue distribution using Marchenko-Pasur Law. See wiki: https://en.wikipedia.org/wiki/Marchenko%E2%80%93Pastur_distribution
"""
marchenko_thresh(n,B) = ((1-sqrt(n/B))^2,(1+sqrt(n/B))^2)

"""
bootstrap(model::PopulationModel, Z, n=100, dim=2)
Fit model `n` times sampling (with replacement) from dimension `dim`.
"""
function bootstrap(model::PopulationModel, Z::AbstractArray{T,N}, n=100, dim=2) where {T,N}
    nsamples = size(Z,dim)
    sample_inds = zeros(Int,nsamples)
    f(_) = begin
        m_i = deepcopy(model)
        sample_inds .= sample(1:nsamples, nsamples, replace=true)
        fit!(m_i, Z[to_indices(Z,Tuple(x != dim ? Colon() : sample_inds for x in 1:N))...])
        return m_i
    end
    return pmap(f,1:n)
end

"""
function confint{T<:PopulationModel}(f::Function,
model::T,
models::AbstractArray{T,1};
α=0.05)
Get confidence interval of parameter from sample of `models` where the
parameter is defined by `f(model)`
"""
function confint(f::Function,
                 model::T,
                 models::AbstractArray{T,1};
                 α=0.05) where T<:PopulationModel
    params = f.(models)
    reorder!.(params, [f(model)])
    lo = zeros(size(params[1])...)
    hi = copy(lo)
    α2 = (α/2)*100
    for i in eachindex(params[1])
        y = cat(1,[x[i] for x in params])
        lo[i] = StatsBase.percentile(y, α2)
        hi[i] = StatsBase.percentile(y, 100-α2)
    end
    return lo,hi
end

function confint(model::T, models::AbstractArray{T,1}; kwargs...) where T<:PopulationModel
    confint(weights, model, models; kwargs...)
end

stability(f::Function, model, models) = mean(similarity(f, model, models),2 )

"""
similarity(f::Function=weights, model, models)
Compare a model to other models using the function `f` (defaults to weights).
"""
function similarity(f::Function, model, models)
    params = f.(models)
    orig = f(model)
    PopulationDynamics.reorder!.(params, [orig])
    cors = VectorOfArray(diag.(cor.([orig], params)))
    return cors
end

similarity(model,models) = similarity(weights, model, models)

"""
get_sig_membs(model::P, btstrp::AbstractArray{P,1}, sigf = (lo,hi) -> (lo .> 0);
α=0.05)
Find the members of the model that consistently contribute to the model under
subsampling (usually bootstrap).
"""
function get_sig_membs(model::P, btstrp::AbstractArray{P,1}, sigf = (lo,hi) -> (lo .> 0);
                       α=0.05) where P<:PopulationModel
    V = weights(model)
    lo,hi=confint(weights, model, btstrp, α=α)
    sign_membs = falses(size(V,2))
    for i in indices(V,2)
        sign_membs[i] = sigf(lo[:,i],hi[:,i])
    end
    sign_membs
end


text:
# http://stronglyconvex.com/blog/accelerated-gradient-descent.html
# TODO: Need to specify alphamax on each iteration
# Flip notation relative to Duckworth
# Start with x_{0}
# y_{t} = x_{t - 1} - alpha g(x_{t - 1})
# If converged, return y_{t}
# x_{t} = y_{t} + (t - 1.0) / (t + 2.0) * (y_{t} - y_{t - 1})

# L should be function or any other callable
immutable AcceleratedGradientDescent{L} <: Optimizer
    linesearch!::L
end

Base.summary(::AcceleratedGradientDescent) = "Accelerated Gradient Descent"

function AcceleratedGradientDescent(; linesearch = LineSearches.HagerZhang())
    AcceleratedGradientDescent(linesearch)
end

type AcceleratedGradientDescentState{T,N}
    x::Array{T,N}
    x_previous::Array{T,N}
    f_x_previous::T
    iteration::Int
    y::Array{T,N}
    y_previous::Array{T,N}
    s::Array{T,N}
    @add_linesearch_fields()
end

function initial_state{T}(method::AcceleratedGradientDescent, options, d, initial_x::Array{T})
    value_gradient!(d, initial_x)
    AcceleratedGradientDescentState(copy(initial_x), # Maintain current state in state.x
                         copy(initial_x), # Maintain previous state in state.x_previous
                         T(NaN), # Store previous f in state.f_x_previous
                         0, # Iteration
                         copy(initial_x), # Maintain intermediary current state in state.y
                         similar(initial_x), # Maintain intermediary state in state.y_previous
                         similar(initial_x), # Maintain current search direction in state.s
                         @initial_linesearch()...) # Maintain a cache for line search results in state.lsr
end

function update_state!{T}(d, state::AcceleratedGradientDescentState{T}, method::AcceleratedGradientDescent)
    n = length(state.x)
    state.iteration += 1
    # Search direction is always the negative gradient
    @simd for i in 1:n
        @inbounds state.s[i] = -gradient(d, i)
    end

    # Determine the distance of movement along the search line
    lssuccess = perform_linesearch!(state, method, d)

    # Record previous state
    copy!(state.x_previous, state.x)

    # Make one move in the direction of the gradient
    copy!(state.y_previous, state.y)
    @simd for i in 1:n
        @inbounds state.y[i] = state.x[i] + state.alpha * state.s[i]
    end

    # Update current position with Nesterov correction
    scaling = (state.iteration - 1) / (state.iteration + 2)
    @simd for i in 1:n
        @inbounds state.x[i] = state.y[i] + scaling * (state.y[i] - state.y_previous[i])
    end

    lssuccess == false # break on linesearch error
end


text:
export iszero_row, modular_hnf, howell_form, _hnf_modular, kernel_mod

import Nemo.matrix

import Base.vcat

#

function matrix(A::Array{fmpz, 2})
  m = matrix(FlintZZ, A)
  return m
end

function matrix(A::Array{T, 2}) where T <: RingElem
  r, c = size(A)
  (r < 0 || c < 0) && error("Array must be non-empty")
  m = matrix(parent(A[1, 1]), A)
  return m
end

function matrix(A::Array{T, 1}) where T <: RingElem
  return matrix(reshape(A,length(A),1))
end

function Array(a::fmpz_mat; S::Type{T} = fmpz) where T
  A = Array{T}(rows(a), cols(a))
  for i = 1:rows(a)
    for j = 1:cols(a)
      A[i,j] = T(a[i,j])
    end
  end
  return A
end

function iszero_row(M::fmpz_mat, i::Int)
  for j = 1:cols(M)
    if M[i,j] != 0
      return false
    end
  end
  return true
end

function iszero_row(M::nmod_mat, i::Int)
  zero = UInt(0)
  for j in 1:cols(M)
    t = ccall((:nmod_mat_get_entry, :libflint), Base.GMP.Limb, (Ref{nmod_mat}, Int, Int), M, i - 1, j - 1)
    if t != zero
      return false
    end
  end
  return true
end

function iszero_row(M::MatElem{T}, i::Int) where T
  for j in 1:cols(M)
    if !iszero(M[i,j])
      return false
    end
  end
  return true
end

function iszero_row(M::Array{T, 2}, i::Int) where T <: Integer
  for j = 1:Base.size(M, 2)
    if M[i,j] != 0
      return false
    end
  end
  return true
end

function iszero_row(M::Array{fmpz, 2}, i::Int)
  for j = 1:Base.size(M, 2)
    if M[i,j] != 0
      return false
    end
  end
  return true
end

function iszero_row(M::Array{T, 2}, i::Int) where T <: RingElem
  for j in 1:Base.size(M, 2)
    if !iszero(M[i,j])
      return false
    end
  end
  return true
end

function divexact!(a::fmpz_mat, b::fmpz_mat, d::fmpz)
  ccall((:fmpz_mat_scalar_divexact_fmpz, :libflint), Void,
               (Ref{fmpz_mat}, Ref{fmpz_mat}, Ref{fmpz}), a, a, d)
end

function mul!(a::fmpz_mat, b::fmpz_mat, c::fmpz)
  ccall((:fmpz_mat_scalar_mul_fmpz, :libflint), Void,
                  (Ref{fmpz_mat}, Ref{fmpz_mat}, Ref{fmpz}), a, b, c)
end

#computes (hopefully) the hnf for vcat(a*I, m) and returns ONLY the
#non-singular part. By definition, the result wil have full rank
#
#Should be rewritten to use Howell and lifting rather the big HNF
#
function modular_hnf(m::fmpz, a::fmpz_mat, shape::Symbol = :upperright)
  c = vcat(parent(a)(m), a)
  n = cols(a)
  w = view(c, n+1, 1, 2*n, n)
  ccall((:fmpz_mat_scalar_mod_fmpz, :libflint), Void, (Ref{fmpz_mat}, Ref{fmpz_mat}, Ref{fmpz}), w, w, m)
  if shape == :lowerleft
    c = _hnf(c, shape)
    return sub(c, n+1:2*n, 1:n)
  else
    c = hnf(c)
    c = sub(c, 1:n, 1:n)
  end
end

function _hnf(x::fmpz_mat, shape::Symbol = :upperright)
  if shape == :lowerleft
    h = hnf(_swapcols(x))
    _swapcols!(h)
    _swaprows!(h)
    return h::fmpz_mat
  end
  return hnf(x)::fmpz_mat
end

function _hnf_modular_eldiv(x::fmpz_mat, m::fmpz, shape::Symbol = :upperright)
  if shape == :lowerleft
    h = hnf_modular_eldiv!(_swapcols(x), m)
    _swapcols!(h)
    _swaprows!(h)
    return h
  elseif shape == :upperright
    return hnf_modular_eldiv(x, m)
  else
    error("shape $shape not supported")
  end
end

function hnf_modular_eldiv!(x::fmpz_mat, d::fmpz)
   (rows(x) < cols(x)) &&
                error("Matrix must have at least as many rows as columns")
   ccall((:fmpz_mat_hnf_modular_eldiv, :libflint), Void,
                (Ref{fmpz_mat}, Ref{fmpz}), x, d)
   return x
end

function ishnf(x::fmpz_mat, shape::Symbol)
  if shape == :upperright
    return ishnf(x)
  elseif shape == :lowerleft
    r = rows(x)
    i = 0
    j_old = cols(x) + 1

    for i in rows(x):-1:1

      if iszero_row(x, i)
        break
      end

      j = cols(x)
      while iszero(x[i, j])
        j = j - 1
      end
      x[i, j] < 0 && return false
      j >= j_old && return false
      for k in i+1:r
        x[k, j] < 0 && return false
        x[k, j] >= x[i, j] && return false
      end
      j_old = j
      i = i - 1
    end

    for k in i:-1:1
      !iszero_row(x, k) && return false
    end
    return true
  end
end

function _swaprows(x::fmpz_mat)
  y = deepcopy(x)
  _swaprows!(y)
  return y
end

function _swapcols(x::fmpz_mat)
  y = deepcopy(x)
  _swapcols!(y)
  return y
end

function _swaprows!(x::fmpz_mat)
  r = rows(x)
  c = cols(x)

  if r == 1
    return x
  end

  if r % 2 == 0
    for i in 1:div(r,2)
      for j = 1:c
        # we swap x[i,j] <-> x[r-i+1,j]
        s = ccall((:fmpz_mat_entry, :libflint), Ptr{fmpz}, (Ref{fmpz_mat}, Int, Int), x, i - 1, j - 1)
        t = ccall((:fmpz_mat_entry, :libflint), Ptr{fmpz}, (Ref{fmpz_mat}, Int, Int), x, (r - i + 1) - 1, j - 1)
        ccall((:fmpz_swap, :libflint), Void, (Ptr{fmpz}, Ptr{fmpz}), t, s)
      end
    end
  else
    for i in 1:div(r-1,2)
      for j = 1:c
        # we swap x[i,j] <-> x[r-i+1,j]
        s = ccall((:fmpz_mat_entry, :libflint), Ptr{fmpz}, (Ref{fmpz_mat}, Int, Int), x, i - 1, j - 1)
        t = ccall((:fmpz_mat_entry, :libflint), Ptr{fmpz}, (Ref{fmpz_mat}, Int, Int), x, (r - i + 1) - 1, j - 1)
        ccall((:fmpz_swap, :libflint), Void, (Ptr{fmpz}, Ptr{fmpz}), t, s)
      end
    end
  end
  nothing
end

function _swaprows!(x::fmpz_mat, i::Int, j::Int)
  ccall((:_fmpz_mat_swap_rows, :libflint), Void, (Ref{fmpz_mat}, Int, Int), x, i-1, j-1)
  nothing
end

function _swaprows!(x::nmod_mat, i::Int, j::Int)
  ccall((:_nmod_mat_swap_rows, :libflint), Void, (Ref{nmod_mat}, Int, Int), x, i-1, j-1)
  nothing
end

function _swaprows!(x::nmod_mat)
  r = rows(x)
  c = cols(x)

  if r == 1
    return nothing
  end

  if r % 2 == 0
    for i in 1:div(r,2)
      for j = 1:c
        # we swap x[i,j] <-> x[r-i+1,j]
        s = ccall((:nmod_mat_get_entry, :libflint), Base.GMP.Limb, (Ref{nmod_mat}, Int, Int), x, i - 1, j - 1)
        t = ccall((:nmod_mat_get_entry, :libflint), Base.GMP.Limb, (Ref{nmod_mat}, Int, Int), x, (r - i + 1) - 1, j - 1)
        set_entry!(x, r - i + 1, j, s)
        set_entry!(x, i, j, t)
      end
    end
  else
    for i in 1:div(r-1,2)
      for j = 1:c
        # we swap x[i,j] <-> x[r-i+1,j]
        s = ccall((:nmod_mat_get_entry, :libflint), Base.GMP.Limb, (Ref{nmod_mat}, Int, Int), x, i - 1, j - 1)
        t = ccall((:nmod_mat_get_entry, :libflint), Base.GMP.Limb, (Ref{nmod_mat}, Int, Int), x, (r - i + 1) - 1, j - 1)
        set_entry!(x, i, j, t)
        set_entry!(x, r - i + 1, j, s)
      end
    end
  end
  nothing
end

function _swapcols!(x::nmod_mat)
  r = rows(x)
  c = cols(x)

  if c == 1
    return nothing
  end

  if c % 2 == 0
    for i in 1:div(c,2)
      for j = 1:r
        # swap x[j,i] <-> x[j,c-i+1]
        s = ccall((:nmod_mat_get_entry, :libflint), Base.GMP.Limb, (Ref{nmod_mat}, Int, Int), x, j - 1, i - 1)
        t = ccall((:nmod_mat_get_entry, :libflint), Base.GMP.Limb, (Ref{nmod_mat}, Int, Int), x, j - 1, (c - i + 1 ) - 1)
        set_entry!(x, j, i, t)
        set_entry!(x, j, c - i + 1, s)
      end
    end
  else
    for i in 1:div(c-1,2)
      for j = 1:r
        # swap x[j,i] <-> x[j,c-i+1]
        s = ccall((:nmod_mat_get_entry, :libflint), Base.GMP.Limb, (Ref{nmod_mat}, Int, Int), x, j - 1, i - 1)
        t = ccall((:nmod_mat_get_entry, :libflint), Base.GMP.Limb, (Ref{nmod_mat}, Int, Int), x, j - 1, (c - i + 1 ) - 1)
        set_entry!(x, j, i, t)
        set_entry!(x, j, c - i + 1, s)
      end
    end
  end
  nothing
end

function _swapcols!(x::fmpz_mat)
  r = rows(x)
  c = cols(x)

  if c == 1
    return x
  end

  if c % 2 == 0
    for i in 1:div(c,2)
      for j = 1:r
        # swap x[j,i] <-> x[j,c-i+1]
        s = ccall((:fmpz_mat_entry, :libflint), Ptr{fmpz}, (Ref{fmpz_mat}, Int, Int), x, j - 1, i - 1)
        t = ccall((:fmpz_mat_entry, :libflint), Ptr{fmpz}, (Ref{fmpz_mat}, Int, Int), x, j - 1, (c - i + 1 ) - 1)
        ccall((:fmpz_swap, :libflint), Void, (Ptr{fmpz}, Ptr{fmpz}), t, s)
      end
    end
  else
    for i in 1:div(c-1,2)
      for j = 1:r
        # swap x[j,i] <-> x[j,c-i+1]
        s = ccall((:fmpz_mat_entry, :libflint), Ptr{fmpz}, (Ref{fmpz_mat}, Int, Int), x, j - 1, i - 1)
        t = ccall((:fmpz_mat_entry, :libflint), Ptr{fmpz}, (Ref{fmpz_mat}, Int, Int), x, j - 1, (c - i + 1 ) - 1)
        ccall((:fmpz_swap, :libflint), Void, (Ptr{fmpz}, Ptr{fmpz}), t, s)
      end
    end
  end
  nothing
end

function _swapcols(x::Generic.Mat)
  z = deepcopy(x)
  _swapcols!(z)
  return z
end

function _swapcols!(x::Generic.Mat)
  r = rows(x)
  c = cols(x)
  t = base_ring(x)(0)

  if c == 1
    return x
  end

  if c % 2 == 0
    for i in 1:div(c,2)
      for j = 1:r
        # swap x[j,i] <-> x[j,c-i+1]
        x[j, i], x[j, c - i + 1] = x[j, c - i + 1], x[j, i]
      end
    end
  else
    for i in 1:div(c-1,2)
      for j = 1:r
        # swap x[j,i] <-> x[j,c-i+1]
        x[j, i], x[j, c - i + 1] = x[j, c - i + 1], x[j, i]
      end
    end
  end
  nothing
end

function _swaprows(x::Generic.Mat)
  z = deepcopy(x)
  _swaprows(z)
  return z
end

function _swaprows!(x::Generic.Mat)
  r = rows(x)
  c = cols(x)

  if r == 1
    return x
  end

  if r % 2 == 0
    for i in 1:div(r,2)
      for j = 1:c
        # we swap x[i,j] <-> x[r-i+1,j]
        x[i, j], x[r - i + 1, j] = x[r - i + 1, j], x[i, j]
      end
    end
  else
    for i in 1:div(r-1,2)
      for j = 1:c
        x[i, j], x[r - i + 1, j] = x[r - i + 1, j], x[i, j]
        # we swap x[i,j] <-> x[r-i+1,j]
      end
    end
  end
  nothing
end
################################################################################
#
################################################################################

function maximum(f::typeof(abs), a::fmpz_mat)
  m = ccall((:fmpz_mat_entry, :libflint), Ptr{fmpz}, (Ref{fmpz_mat}, Int, Int), a, 0,0)
  for i=1:rows(a)
    for j=1:cols(a)
      z = ccall((:fmpz_mat_entry, :libflint), Ptr{fmpz}, (Ref{fmpz_mat}, Int, Int), a, i-1, j-1)
      if ccall((:fmpz_cmpabs, :libflint), Cint, (Ptr{fmpz}, Ptr{fmpz}), m, z) < 0
        m = z
      end
    end
  end
  r = fmpz()
  ccall((:fmpz_abs, :libflint), Void, (Ref{fmpz}, Ptr{fmpz}), r, m)
  return r
end

function maximum(a::fmpz_mat)
  m = ccall((:fmpz_mat_entry, :libflint), Ptr{fmpz}, (Ref{fmpz_mat}, Int, Int), a, 0,0)
  for i=1:rows(a)
    for j=1:cols(a)
      z = ccall((:fmpz_mat_entry, :libflint), Ptr{fmpz}, (Ref{fmpz_mat}, Int, Int), a, i-1, j-1)
      if ccall((:fmpz_cmp, :libflint), Cint, (Ptr{fmpz}, Ptr{fmpz}), m, z) < 0
        m = z
      end
    end
  end
  r = fmpz()
  ccall((:fmpz_set, :libflint), Void, (Ref{fmpz}, Ptr{fmpz}), r, m)
  return r
end

function minimum(a::fmpz_mat)
  m = ccall((:fmpz_mat_entry, :libflint), Ptr{fmpz}, (Ref{fmpz_mat}, Int, Int), a, 0,0)
  for i=1:rows(a)
    for j=1:cols(a)
      z = ccall((:fmpz_mat_entry, :libflint), Ptr{fmpz}, (Ref{fmpz_mat}, Int, Int), a, i-1, j-1)
      if ccall((:fmpz_cmp, :libflint), Cint, (Ptr{fmpz}, Ptr{fmpz}), m, z) > 0
        m = z
      end
    end
  end
  r = fmpz()
  ccall((:fmpz_set, :libflint), Void, (Ref{fmpz}, Ptr{fmpz}), r, m)
  return r
end

function lift_unsigned(a::nmod_mat)
  z = zero_matrix(FlintZZ, rows(a), cols(a))
  ccall((:fmpz_mat_set_nmod_mat_unsigned, :libflint), Void,
          (Ref{fmpz_mat}, Ref{nmod_mat}), z, a)
  return z
end

################################################################################
# possibly a slice or window in fmpz_mat?
# the nr x nc matrix starting in (a,b)
################################################################################

function sub(A::fmpz_mat, r::UnitRange, c::UnitRange)
  return deepcopy(view(A, r, c))
end

################################################################################
#
#  misc.jl : At some point this should migrate to Nemo
#
################################################################################

# compute basis for the right kernel space by calling flint
# look at flint documentation of nmod_mat_nullspace

function _right_kernel(x::nmod_mat)
  z = zero_matrix(base_ring(x), cols(x), max(rows(x),cols(x)))
  n = ccall((:nmod_mat_nullspace, :libflint), Int, (Ref{nmod_mat}, Ref{nmod_mat}), z, x)
  return z,n
end

# compute basis for the left kernel space
# output is array of arrays, which span the kernel

function kernel(a)
  x = transpose(a)
  R = base_ring(a)
  z, n = _right_kernel(x)
  z = transpose(z)
  T = elem_type(base_ring(a))
  ar = typeof(Array{T}(cols(z)))[]
  for i in 1:n
    t = Array{T}(cols(z))
    for j in 1:cols(z)
      t[j] = R(z[i, j])
    end
    push!(ar,t)
  end
  return ar
end

function right_kernel(a::MatElem)
  R = base_ring(a)
  if typeof(R) <: Generic.ResRing{fmpz}
    nn, zz = _right_kernel(a)
  else
    zz, nn = nullspace(a)
  end
  #TODO: There is some inconsistency between returning the nullity or the kernel first...
  if typeof(nn) <: Integer
    n = nn
    z = zz
  else
    n = zz
    z = nn
  end
  T = elem_type(base_ring(a))
  ar = typeof(Array{T}(rows(z)))[]
  for i in 1:n
    t = Array{T}(rows(z))
    for j in 1:rows(z)
      t[j] = R(z[j, i])
    end
    push!(ar,t)
  end
  return ar
end

left_kernel(a::MatElem) = right_kernel(transpose(a))

function lift(a::Generic.Mat{Generic.Res{fmpz}})
  z = zero_matrix(FlintZZ, rows(a), cols(a))
  for i in 1:rows(a)
    for j in 1:cols(a)
      z[i, j] = lift(a[i, j])
    end
  end
  return z
end

function _rref(a::Generic.Mat{Generic.Res{fmpz}})
  m = modulus(base_ring(a))
  b = zero_matrix(FlintZZ, rows(a), cols(a))
  # I actually don't know if this is necessary
  for i in 1:rows(b)
    for j in 1:cols(b)
      b[i,j] = lift(a[i,j]) % m
    end
  end

  # fmpz_mat_rref_mod assumes that input is reduced modulo m
  r = ccall((:fmpz_mat_rref_mod, :libflint), Int, (Ptr{Void}, Ref{fmpz_mat}, Ref{fmpz}), C_NULL, b, m)
  return r, parent(a,false)(b)
end

function _rref(a::nmod_mat)
  b = rref(a)
  # TODO: Clean up one we use new Nemo version.
  if length(b) == 1
    return rank(b), b
  else
    return b
  end
end

function _right_kernel(a::Generic.Mat{Generic.Res{fmpz}})
  r, b = _rref(a)
  pivots = Array{Int}(r)
  nonpivots = Array{Int}(cols(b) - r)
  X = zero_matrix(base_ring(a),cols(b),cols(b) - r)

  if r == 0
    return vcat(identity_matrix(FlintZZ, cols(b) - r), zero_matrix(FlintZZ, r, cols(b) - r)), cols(b)
  elseif !((cols(b) - r) == 0)
    i = 1
    j = 1
    k = 1
    for i in 1:r
      while b[i,j] == 0
        nonpivots[k] = j
        k += 1
        j += 1
      end
      pivots[i] = j
      j += 1
    end
    while k <= cols(b) - r
      nonpivots[k] = j
      k += 1
      j += 1
    end

    for i in 1:cols(b) - r
      for j in 1:r
        X[pivots[j],i] = - b[j,nonpivots[i]]
      end
      X[nonpivots[i],i] = 1
    end
  end
  return X, cols(b) - r
end

function kernel_mod(a::fmpz_mat, m::fmpz)
  b = deepcopy(a)
  for i in 1:rows(b)
    for j in 1:cols(b)
      b[i,j] = b[i,j] % m
    end
  end

  # fmpz_mat_rref_mod assumes that input is reduced modulo m
  r = ccall((:fmpz_mat_rref_mod, :libflint), Int, (Ptr{Void}, Ref{fmpz_mat}, Ref{fmpz}), C_NULL, b, m)
  pivots = Array{Int}(r)
  nonpivots = Array{Int}(cols(b) - r)
  X = zero_matrixSpace(FlintZZ,cols(b),cols(b))

  if r == 0
    return identity_matrix(FlintZZ, cols(b))
  elseif !((cols(b) - r) == 0)
    i = 1
    j = 1
    k = 1
    for i in 1:r
      while b[i,j] == 0
        nonpivots[k] = j
        k += 1
        j += 1
      end
      pivots[i] = j
      j += 1
    end
    while k <= cols(b) - r
      nonpivots[k] = j
      k += 1
      j += 1
    end

    for i in 1:cols(b) - r
      for j in 1:r
        X[pivots[j],i] = - FlintZZ(b[j,nonpivots[i]])
      end
      X[nonpivots[i],i] = FlintZZ(1)
    end
  end
  return X, r
end

# Another kernel function
function _kernel(x::fmpz_mat)
  H, U = hnf_with_transform(x)
  i = 1
  for i in 1:rows(H)
    if iszero_row(H, i)
      break
    end
  end
  return sub(U, i:rows(U), 1:cols(U))
end

################################################################################
#
#  Copy matrix into another matrix
#
################################################################################

# Copy B into A at position (i, j)
function _copy_matrix_into_matrix(A::fmpz_mat, i::Int, j::Int, B::fmpz_mat)
  for k in 0:rows(B) - 1
    for l in 0:cols(B) - 1
      d = ccall((:fmpz_mat_entry, :libflint),
                Ptr{fmpz}, (Ref{fmpz_mat}, Int, Int), B, k, l)
      t = ccall((:fmpz_mat_entry, :libflint),
                Ptr{fmpz}, (Ref{fmpz_mat}, Int, Int), A, i - 1 + k, j - 1 + l)
      ccall((:fmpz_set, :libflint), Void, (Ptr{fmpz}, Ptr{fmpz}), t, d)
    end
  end
end

doc"""
    isposdef(a::fmpz_mat) -> Bool

> Tests if $a$ positive definite by testing if all principal minors
> have positive determinant.
"""
function isposdef(a::fmpz_mat)
  for i=1:rows(a)
    if det(sub(a, 1:i, 1:i)) <= 0
      return false
    end
  end
  return true
end

#scales the i-th column of a by 2^d[1,i]
function mult_by_2pow_diag!(a::Array{BigFloat, 2}, d::fmpz_mat)
  s = size(a)
  R = RealRing()::_RealRing
  tmp_mpz::BigInt = R.z1
  for i = 1:s[1]
    for j = 1:s[2]
      e = ccall((:mpfr_get_z_2exp, :libmpfr), Clong, (Ref{BigInt}, Ref{BigFloat}), tmp_mpz, a[i,j])
      ccall((:mpfr_set_z_2exp, :libmpfr), Void, (Ref{BigFloat}, Ref{BigInt}, Clong, Int32), a[i,j], tmp_mpz, e+Clong(Int(d[1,j])), __get_rounding_mode())
    end
  end
end

#converts BigFloat -> fmpz via round(a*2^l), in a clever(?) way
function round_scale(a::Array{BigFloat, 2}, l::Int)
  s = size(a)
  b = zero_matrix(FlintZZ, s[1], s[2])
  return round_scale!(b, a, l)
end

function round_scale!(b::fmpz_mat, a::Array{BigFloat, 2}, l::Int)
  s = size(a)
  R = RealRing()::_RealRing

  local tmp_mpz::BigInt, tmp_fmpz::fmpz
  tmp_mpz = R.z1
  tmp_fmpz = R.zz1
  tmp_mpfr = deepcopy(a[1,1])  #cannot use the R.?? tmp variable as it may/will
                               #have the wrong precision

  rd = __get_rounding_mode()
  for i = 1:s[1]
    for j = 1:s[2]
      e = a[i,j].exp
      a[i,j].exp += l
      ccall((:mpfr_round, :libmpfr), Int32, (Ref{BigFloat}, Ref{BigFloat}, Int32), tmp_mpfr, a[i,j], rd)
      a[i,j].exp = e
      f = ccall((:mpfr_get_z_2exp, :libmpfr), Clong, (Ref{BigInt}, Ref{BigFloat}),
        tmp_mpz, tmp_mpfr)
      ccall((:fmpz_set_mpz, :libflint), Void, (Ref{fmpz}, Ref{BigInt}), tmp_fmpz, tmp_mpz)
      if f > 0
        ccall((:fmpz_mul_2exp, :libflint), Void, (Ref{fmpz}, Ref{fmpz}, UInt), tmp_fmpz, tmp_fmpz, f)
      else
        ccall((:fmpz_tdiv_q_2exp, :libflint), Void, (Ref{fmpz}, Ref{fmpz}, UInt), tmp_fmpz, tmp_fmpz, -f);
      end
      setindex!(b, tmp_fmpz, i, j)
    end
  end
  return b
end

function shift!(g::fmpz_mat, l::Int)
  for i=1:rows(g)
    for j=1:cols(g)
      z = ccall((:fmpz_mat_entry, :libflint), Ptr{fmpz}, (Ref{fmpz_mat}, Int, Int), g, i-1, j-1)
      if l > 0
        ccall((:fmpz_mul_2exp, :libflint), Void, (Ptr{fmpz}, Ptr{fmpz}, Int), z, z, l)
      else
        ccall((:fmpz_tdiv_q_2exp, :libflint), Void, (Ptr{fmpz}, Ptr{fmpz}, Int), z, z, -l)
      end
    end
  end
  return g
end

doc"""
***
    mod!(M::fmpz_mat, p::fmpz)
> Reduces every entry modulo $p$ in-place, ie. applies the mod function to every entry.
"""
function mod!(M::fmpz_mat, p::fmpz)
  for i=1:rows(M)
    for j=1:cols(M)
      z = ccall((:fmpz_mat_entry, :libflint), Ptr{fmpz}, (Ref{fmpz_mat}, Int, Int), M, i - 1, j - 1)
      ccall((:fmpz_mod, :libflint), Void, (Ptr{fmpz}, Ptr{fmpz}, Ref{fmpz}), z, z, p)
    end
  end
  nothing
end

doc"""
***
    mod(M::fmpz_mat, p::fmpz) -> fmpz_mat
> Reduces every entry modulo $p$, ie. applies the mod function to every entry.
"""
function mod(M::fmpz_mat, p::fmpz)
  N = deepcopy(M)
  mod!(N, p)
  return N
end

doc"""
***
    vcat(A::Array{Generic.Mat, 1}) -> Generic.Mat
    vcat(A::Array{fmpz_mat}, 1}) -> fmpz_mat
> Forms a big matrix my vertically concatenating the matrices in $A$.
> All component matrices need to have the same number of columns.
"""
function vcat(A::Array{Generic.Mat{T}, 1}) where T
  if any(x->cols(x) != cols(A[1]), A)
    error("Matrices must have same number of columns")
  end
  M = zero_matrix(base_ring(A[1]), sum(rows, A), cols(A[1]))
  s = 0
  for i=A
    for j=1:rows(i)
      for k=1:cols(i)
        M[s+j, k] = i[j,k]
      end
    end
    s += rows(i)
  end
  return M
end

function vcat(A::Array{fmpz_mat, 1})
  if any(x->cols(x) != cols(A[1]), A)
    error("Matrices must have same number of columns")
  end
  M = zero_matrix(base_ring(A[1]), sum(rows, A), cols(A[1]))
  s = 0
  for i=A
    for j=1:rows(i)
      for k=1:cols(i)
        M[s+j, k] = i[j,k]
      end
    end
    s += rows(i)
  end
  return M
end

function vcat(A::Array{nmod_mat, 1})
  if any(x->cols(x) != cols(A[1]), A)
    error("Matrices must have same number of columns")
  end
  M = zero_matrix(base_ring(A[1]), sum(rows, A), cols(A[1]))
  s = 0
  for i=A
    for j=1:rows(i)
      for k=1:cols(i)
        M[s+j, k] = i[j,k]
      end
    end
    s += rows(i)
  end
  return M
end

################################################################################
#
#  Smith normal form with trafo
#
################################################################################

#=
g, e,f = gcdx(a, b)
U = [1 0 ; -divexact(b, g)*f 1]*[1 1; 0 1];
V = [e -divexact(b, g) ; f divexact(a, g)];

then U*[ a 0; 0 b] * V = [g 0 ; 0 l]
=#
doc"""
***
  snf_with_transform(A::fmpz_mat, l::Bool = true, r::Bool = true) -> fmpz_mat, fmpz_mat, fmpz_mat

> Given some integer matrix A, compute the Smith normal form (elementary
> divisor normal form) of A. If l and/ or r are true, then the corresponding
> left and/ or right transformation matrices are computed as well.
"""
function snf_with_transform(A::fmpz_mat, l::Bool = true, r::Bool = true)
  if r
    R = identity_matrix(FlintZZ, cols(A))
  end

  if l
    L = identity_matrix(FlintZZ, rows(A))
  end
  # TODO: if only one trafo is required, start with the HNF that does not
  #       compute the trafo
  #       Rationale: most of the work is on the 1st HNF..
  #
  S = deepcopy(A)
  while !isdiag(S)
    if l
      S, T = hnf_with_transform(S)
      L = T*L
    else
      S = hnf(S)
    end

    if isdiag(S)
      break
    end
    if r
      S, T = hnf_with_transform(S')
      R = T*R
    else
      S = hnf(S')
    end
    S = S'
  end
  #this is probably not really optimal...
  for i=1:min(rows(S), cols(S))
    if S[i,i] == 1
      continue
    end
    for j=i+1:min(rows(S), cols(S))
      if S[j,j] == 0
        continue
      end
      if S[i,i] != 0 && S[j,j] % S[i,i] == 0
        continue
      end
      g, e,f = gcdx(S[i,i], S[j,j])
      a = divexact(S[i,i], g)
      S[i,i] = g
      b = divexact(S[j,j], g)
      S[j,j] *= a
      if l
        # U = [1 0; -b*f 1] * [ 1 1; 0 1] = [1 1; -b*f -b*f+1]
        # so row i and j of L will be transformed. We do it naively
        # those 2x2 transformations of 2 rows should be a c-primitive
        # or at least a Nemo/Hecke primitive
        for k=1:cols(L)
          x = -b*f
#          L[i,k], L[j,k] = L[i,k]+L[j,k], x*L[i,k]+(x+1)*L[j,k]
          L[i,k], L[j,k] = L[i,k]+L[j,k], x*(L[i,k]+L[j,k])+L[j,k]
        end
      end
      if r
        # V = [e -b ; f a];
        # so col i and j of R will be transformed. We do it naively
        # careful: at this point, R is still transposed
        for k=1:rows(R)
          R[i, k], R[j, k] = e*R[i,k]+f*R[j,k], -b*R[i,k]+a*R[j,k]
        end
      end
    end
  end

  if l
    if r
      return S, L, R'
    else
      # last is dummy
      return S, L, L
    end
  elseif r
    # second is dummy
    return S, R, R'
  else
    # last two are dummy
    return S, S, S
  end
end

function Base.nullspace(M::nmod_mat)
  R = base_ring(M)
  if isprime(fmpz(modulus(R)))
    k = zero_matrix(R, cols(M), cols(M))
    n = ccall((:nmod_mat_nullspace, :libflint), Int, (Ref{nmod_mat}, Ref{nmod_mat}), k, M)
    return (k, n)
  end

  N = hcat(M', identity_matrix(R, cols(M)))
  ex = 0
  if rows(N) < cols(N)
    ex = cols(N) - rows(N)
    N = vcat(N, zero_matrix(R, ex, cols(N)))
  end
  H = howell_form(N)
  nr = 1
  while nr <= rows(H) && !iszero_row(H, nr)
    nr += 1
  end
  nr -= 1
  h = sub(H, 1:nr, 1:rows(M))
  for i=1:rows(h)
    if iszero_row(h, i)
      k = sub(H, i:rows(h), rows(M)+1:cols(H))
      return k', rows(k)
    end
  end
  return zero_matrix(R,rows(M),0),0
end

function lift(M::FmpzMatSpace, Mp::Union{nmod_mat,Generic.Mat{Generic.Res{fmpz}}})
  @assert M.cols == cols(Mp) && M.rows == rows(Mp)
  N = M()
  for i=1:M.rows
    for j=1:M.cols
      N[i,j] = lift(Mp[i,j])
    end
  end
  return N
end

################################################################################
#
#  Is diagonal
#
################################################################################

doc"""
***
  isdiag(A::fmpz_mat)

> Tests if A is diagonal.
"""
function isdiag(A::fmpz_mat)
  for i = 1:cols(A)
    for j = 1:rows(A)
      if i != j && !iszero(A[j, i])
        return false
      end
    end
  end
  return true
end

################################################################################
#
#  Triangular solving
#
################################################################################

# Solves A x = b for A upper triangular m\times n matrix and b m\times 1.
function solve_ut(A::MatElem{T}, b::MatElem{T}) where T
  m = rows(A)
  n = cols(A)
  @assert m == rows(b)
  @assert m <= n
  x = similar(A, n, 1)
  pivot_cols = Vector{Int}()
  r = 0
  last_pivot = n + 1
  for i = m:-1:1
    for j = 1:last_pivot - 1
      if iszero(A[i, j])
        continue
      end
      x[j, 1] = b[i, 1]
      for k = 1:r
        x[j, 1] -= A[i, pivot_cols[k]]*x[pivot_cols[k], 1]
      end
      x[j, 1] *= inv(A[i, j])
      last_pivot = j
      r += 1
      push!(pivot_cols, j)
      break
    end
  end
  return x
end

# Solves A x = b for A lower triangular m\times n matrix and b m\times 1.
function solve_lt(A::MatElem{T}, b::MatElem{T}) where T
  m = rows(A)
  n = cols(A)
  @assert m == rows(b)
  @assert m <= n
  x = similar(A, n, 1)
  pivot_cols = Vector{Int}()
  r = 0
  last_pivot = 0
  for i = 1:m
    for j = n:-1:last_pivot + 1
      if iszero(A[i, j])
        continue
      end
      x[j, 1] = b[i, 1]
      for k = 1:r
        x[j, 1] -= A[i, pivot_cols[k]]*x[pivot_cols[k], 1]
      end
      x[j, 1] *= inv(A[i, j])
      last_pivot = j
      r += 1
      push!(pivot_cols, j)
      break
    end
  end
  return x
end

# =======================================
# Array interface for MatElem
# =======================================
#
# TODO: re-write for special types (fmpz_mat e.g.) to gain efficiency
#

length(A::Nemo.MatElem) = rows(A) * cols(A)
Base.ndims(A::Nemo.MatElem) = 2

function Base.size(A::Nemo.MatElem, n::Int)
  if n == 1
    return rows(A)
  elseif n == 2
    return cols(A)
  elseif n < 1
    error("arraysize: dimension out of range")
  else
    return 1
  end
end

function Base.indices(A::Nemo.MatElem)
  return (Base.OneTo(rows(A)), Base.OneTo(cols(A)))
end

function Base.indices(A::Nemo.MatElem, n::Int)
  return Base.OneTo(size(A, n))
end

function Base.eachindex(A::Nemo.MatElem)
  return Base.OneTo(length(A))
end

function Base.stride(A::Nemo.MatElem, n::Int)
  if n <= 1
    return 1
  elseif n == 2
    return rows(A)
  else
    return length(A)
  end
end

Base.eltype(A::Nemo.MatElem{T}) where T <: Nemo.RingElem = T

getindex(A::Nemo.MatElem, n::Int) = A[1 + ((n-1) % rows(A)), 1 + div((n-1), rows(A))]

function setindex!(A::Nemo.MatElem{T}, n::Int, s::T) where T <: RingElem
  A[1 + ((n-1) % rows(A)), 1 + div((n-1), rows(A))] = s
end

Base.start(A::Nemo.MatElem) = 1
Base.next(A::Nemo.MatElem, i::Int) = A[i], i+1
Base.done(A::Nemo.MatElem, i::Int) = i > length(A)

function setindex!(A::Nemo.MatElem{T}, b::Nemo.MatElem{T}, ::Colon, i::Int) where T
  @assert cols(b) == 1 && rows(b) == rows(A)
  for j=1:rows(A)
    A[j,i] = b[j]
  end
  b
end

function setindex!(A::Nemo.MatElem{T}, b::Nemo.MatElem{T}, i::Int, ::Colon) where T
  @assert rows(b) == 1 && cols(b) == cols(A)
  for j=1:cols(A)
    A[i,j] = b[j]
  end
  b
end

function setindex!(A::Nemo.MatElem, b::Array{<: Any, 1}, ::Colon, i::Int)
  @assert length(b) == rows(A)
  for j=1:rows(A)
    A[j,i] = b[j]
  end
  b
end

function setindex!(A::Nemo.MatElem, b::Array{ <: Any, 1}, i::Int, ::Colon)
  @assert length(b) == cols(A)
  for j=1:cols(A)
    A[i,j] = b[j]
  end
  b
end

function setindex!(A::Nemo.MatElem, b, ::Colon, i::Int)
  for j=1:rows(A)
    A[j,i] = b
  end
  b
end

function setindex!(A::Nemo.MatElem, b, i::Int, ::Colon)
  for j=1:cols(A)
    A[i,j] = b
  end
  b
end

getindex(A::Nemo.MatElem, i::Int, ::Colon) = A[i:i, 1:cols(A)]
getindex(A::Nemo.MatElem, ::Colon, i::Int) = A[1:rows(A), i:i]

function Base.hcat(A::Nemo.MatElem...)
  r = rows(A[1])
  c = cols(A[1])
  R = base_ring(A[1])
  for i=2:length(A)
    @assert rows(A[i]) == r
    @assert base_ring(A[i]) == R
    c += cols(A[i])
  end
  X = zero_matrix(R, r, c)
  o = 1
  for i=1:length(A)
    for j=1:cols(A[i])
      X[:, o] = A[i][:, j]
      o += 1
    end
  end
  return X
end

function Base.vcat(A::Nemo.MatElem...)
  r = rows(A[1])
  c = cols(A[1])
  R = base_ring(A[1])
  for i=2:length(A)
    @assert cols(A[i]) == c
    @assert base_ring(A[i]) == R
    r += rows(A[i])
  end
  X = zero_matrix(R, r, c)
  o = 1
  for i=1:length(A)
    for j=1:rows(A[i])
      X[o, :] = A[i][j, :]
      o += 1
    end
  end
  return X
end

function Base.cat(n::Int, A::Nemo.MatElem...)
  if n==1
    return vcat(A...)
  elseif n==2
    return hcat(A...)
  else
    error("does not make sense here")
  end
end

function Base.cat(dims::Tuple{Int, Int}, A::Nemo.MatElem...)
  @assert dims == (1,2)

  z = [similar(x) for x = A]
  X = z[1]
  for i=1:length(A)
    if i==1
      X = hcat(A[1], z[2:end]...)
    else
      X = vcat(X, hcat(z[1:i-1]..., A[i], z[i+1:end]...))
    end
  end
  return X
end

Base.cat(dims, A::Nemo.MatElem...) = cat(Tuple(dims), A...)

doc"""
    reduce_mod!(A::Nemo.MatElem{T}, B::Nemo.MatElem{T}) where T <: Nemo.FieldElem

> For a reduced row echelon matrix $B$, reduce $A$ modulo $B$, ie. all the pivot
> columns will be zero afterwards.
"""
function reduce_mod!(A::Nemo.MatElem{T}, B::Nemo.MatElem{T}) where T <: Nemo.FieldElem
  @assert isrref(B)
  for h=1:rows(A)
    j = 1
    for i=1:rows(B)
      while iszero(B[i, j])
        j += 1
      end
      A[h, :] -= A[h, j] * B[i, :]
    end
  end
  return A
end

doc"""
    reduce_mod(A::Nemo.MatElem{T}, B::Nemo.MatElem{T}) where T <: Nemo.FieldElem -> MatElem

> For a reduced row echelon matrix $B$, reduce $A$ modulo $B$, ie. all the pivot
> columns will be zero afterwards.
"""
function Nemo.reduce_mod(A::Nemo.MatElem{T}, B::Nemo.MatElem{T}) where T <: Nemo.FieldElem
  C = deepcopy(A)
  reduce_mod!(C, B)
  return C
end

doc"""
    find_pivot(A::Nemo.MatElem{<:Nemo.RingElem}) -> Array{Int, 1}

> Find the pivot-columns of the reduced row echelon matrix $A$
"""
function find_pivot(A::Nemo.MatElem{<:Nemo.RingElem})
  @assert isrref(A)
  p = Int[]
  j = 0
  for i=1:rows(A)
    j += 1
    if j > cols(A)
      return p
    end
    while iszero(A[i,j])
      j += 1
      if j > cols(A)
        return p
      end
    end
    push!(p, j)
  end
  return p
end

doc"""
    cansolve(A::Nemo.MatElem{T}, B::Nemo.MatElem{T}) where T <: Nemo.FieldElem -> Bool, MatElem
> Tries to solve $Ax = B$
"""
function Nemo.cansolve(A::Nemo.MatElem{T}, B::Nemo.MatElem{T}) where T <: Nemo.FieldElem
  R = base_ring(A)
  @assert R == base_ring(B)
  @assert rows(A) == rows(B)
  mu = [A B]
  rk, mu = rref(mu)
  p = find_pivot(mu)
  if any(i->i>cols(A), p)
    return false, B
  end
  sol = zero_matrix(R, cols(A), cols(B))
  for i = 1:length(p)
    for j = 1:cols(B)
      sol[p[i], j] = mu[i, cols(A) + j]
    end
  end
  return true, sol
end

doc"""
    cansolve_with_nullspace(A::Nemo.MatElem{T}, B::Nemo.MatElem{T}) where T <: Nemo.FieldElem -> Bool, MatElem, MatElem
> Tries to solve $Ax = B$, returns a solution and the nullspace.
"""
function Nemo.cansolve_with_nullspace(A::Nemo.MatElem{T}, B::Nemo.MatElem{T}) where T <: Nemo.FieldElem
  R = base_ring(A)
  @assert R == base_ring(B)
  @assert rows(A) == rows(B)
  mu = [A B]
  rk, mu = rref(mu)
  p = find_pivot(mu)
  if any(i->i>cols(A), p)
    return false, B, B
  end
  sol = zero_matrix(R, cols(A), cols(B))
  for i = 1:length(p)
    for j = 1:cols(B)
      sol[p[i], j] = mu[i, cols(A) + j]
    end
  end
  n = zero_matrix(R, cols(A), cols(A) - length(p))
  np = sort(setdiff(1:cols(A), p))
  i = 0
  push!(p, cols(A)+1)
  for j = 1:length(np)
    if np[j] >= p[i+1]
      i += 1
    end
    if i > 0
      n[p[i], j] = -mu[i, np[j]]
    end
    n[np[j], j] = 1
  end
  return true, sol, n
end

#TODO: different to cansolve*(fmpz_mat) is hnf_with_tranformation -> hnf_with_trafo
#maybe (definitely!) agree on one name and combine?

doc"""
    cansolve(A::Nemo.MatElem{T}, B::Nemo.MatElem{T}) where T <: Nemo.RingElem -> Bool, MatElem
> Tries to solve $Ax = B$ where the matrices are defined over a euclidean ring.
"""
function Nemo.cansolve(a::Nemo.MatElem{S}, b::Nemo.MatElem{S}) where S <: Nemo.RingElem
  R = base_ring(a)
  @assert R == base_ring(b)
  @assert rows(a) == rows(b)

  H, T = hnf_with_trafo(transpose(a))
  b = deepcopy(b)
  z = similar(a, cols(b), cols(a))
  l = min(rows(a), cols(a))
  for i = 1:cols(b)
    for j = 1:l
      k = 1
      while k <= cols(H) && iszero(H[j, k])
        k += 1
      end
      if k > cols(H)
        continue
      end
      q, r = divrem(b[k, i], H[j, k])
      if !iszero(r)
        return false, b
      end
      for h = k:cols(H)
        b[h, i] -= q*H[j, h]
      end
      z[i, j] = q
    end
  end
  if !iszero(b)
    return false, b
  end
  return true, transpose(z*T)
end

doc"""
    cansolve_with_nullspace(A::Nemo.MatElem{T}, B::Nemo.MatElem{T}) where T <: Nemo.RingElem -> Bool, MatElem, MatElem
> Tries to solve $Ax = B$ where the matrices are defined over a euclidean ring. If successful,
> a basis for the nullspace is computed as well.
"""
function Nemo.cansolve_with_nullspace(a::Nemo.MatElem{S}, b::Nemo.MatElem{S}) where S <: Nemo.RingElem
  R = base_ring(a)
  @assert R == base_ring(b)
  @assert rows(a) == rows(b)

  H, T = hnf_with_trafo(transpose(a))
  z = similar(a, cols(b), cols(a))
  l = min(rows(a), cols(a))
  for i=1:cols(b)
    for j=1:l
      k = 1
      while k <= cols(H) && iszero(H[j, k])
        k += 1
      end
      if k > cols(H)
        continue
      end
      q, r = divrem(b[k, i], H[j, k])
      if !iszero(r)
        return false, b, b
      end
      for h=k:cols(H)
        b[h, i] -= q*H[j, h]
      end
      z[i, k] = q
    end
  end
  if !iszero(b)
    return false, b, b
  end

  for i = rows(H):-1:1
    for j = 1:cols(H)
      if !iszero(H[i,j])
        N = similar(a, cols(a), rows(H) - i)
        for k = 1:rows(N)
          for l = 1:cols(N)
            N[k,l] = T[rows(T) - l + 1, k]
          end
        end
        return true, transpose(z*T), N
      end
    end
  end
  N =  similar(a, cols(a), 0)

  return true, (z*T), N
end

doc"""
   elementary_divisors(A::fmpz_mat) -> Array{fmpz, 1}
> The elementary divisors of $A$, ie. the diagonal entries of the Smith normal form of $A$.
"""
function elementary_divisors(A::fmpz_mat)
  s = snf(A)
  return [s[i,i] for i=1:min(cols(s), rows(s))]
end

doc"""
   maximal_elementary_divisors(A::fmpz_mat) -> fmpz
> The largest elementary divisor of $A$, ie. the last diagonal entry of the Smith normal form of $A$.
"""
function maximal_elementary_divisor(A::fmpz_mat)
  return elementary_divisors(A)[end]
end


text:
using IPIL8

# --------------------------------------------------------------------------------------------------
@testset "Тесты на регрессию прямой задачи                  " begin

    # Получим стандартные аргументы
    eval(IPIL8.default_parameters());

    # Загрузим зафиксированные данные регрессионных тестов
    include("direct_regression_tests_data.jl");

    # Тестируем внутреннюю функцию, поэтому входные и выходные векторы размерности N-2;
    y = u₀[2:end-1];
    RP = zero(y);

# --------------------------------------------------------------------------------------------------
@testset "Вектор правой части                               " begin

    IPIL8.directRP!(RP, y, 1, h, N, ε, ulₘ, urₘ);

    @test(isapprox(known_RP, RP));

end # testset Вектор правой части
# --------------------------------------------------------------------------------------------------

# --------------------------------------------------------------------------------------------------
@testset "Три диагонали якобиана правой части               " begin

    dl = zeros(N-2);
    dd = zeros(N-1);
    du = zeros(N-2);

    IPIL8.DRP!(dl, dd, du, y, 1, h, N, ε, ulₘ, urₘ)

    @test(isapprox(known_dl, dl));
    @test(isapprox(known_dd, dd));
    @test(isapprox(known_du, du));

end
# --------------------------------------------------------------------------------------------------

# --------------------------------------------------------------------------------------------------
@testset "Решение прямой задачи                             " begin

    solution = zeros(N+1, M+1);

    IPIL8.solve!(solution, u₀, N, h, Tₘ, M, τ, ε, ulₘ, urₘ);

    @test(isapprox(known_u, solution));

end
# --------------------------------------------------------------------------------------------------

end # testset Регрессия прямой задачи
# --------------------------------------------------------------------------------------------------


text:
abstract type AbstractOracle{T<:AbstractTransitionSystem} end

initconfig(oracle::AbstractOracle, args...) = initconfig(system(oracle), args...)

xys(oracle, trees) = Base.Iterators.flatten(xys(oracle, tree) for tree in trees)
xys(oracle, tree::DependencyTree; kwargs...) =
    ((state.cfg, state.G) for state in oracle(tree; kwargs...))

"""
   UnparsableTree

Represents trees that an oracle can't parse.
"""
struct UnparsableTree
    err
end
Base.IteratorSize(pairs::UnparsableTree) = Base.HasLength()
Base.iterate(pairs::UnparsableTree, state...) = nothing
Base.length(::UnparsableTree) = 0

"""
    OracleState

Temporary state for building a transition parse.
"""
struct OracleState{C}
    cfg::C
    A::Vector{TransitionOperator}
    G::Vector{TransitionOperator}
end

include("oracles/exploration.jl")

"""
    TreeOracle

A gold tree for training.

Collects an `oracle`, a gold `tree`, and an exploration `policy`.
"""
struct TreeOracle{O<:AbstractOracle}
    oracle::O
    tree::DependencyTree
    policy

    function TreeOracle(oracle::AbstractOracle, tree::DependencyTree, policy=NeverExplore())
        if projective_only(system(oracle)) && !isprojective(tree)
            UnparsableTree(NonProjectiveGraphError(tree))
        else
            new{typeof(oracle)}(oracle, tree, policy)
        end
    end
end

Base.IteratorSize(o::TreeOracle) = Base.SizeUnknown()

function Base.iterate(o::TreeOracle)
    system, tree, policy = o.oracle.system, o.tree, o.policy
    cfg = initconfig(system, tree)
    state = oracle_state(o, cfg)
    t = policy(state)
    next = t(cfg)
    return (state, next)
end

function Base.iterate(o::TreeOracle, cfg)
    if isfinal(cfg)
        return nothing
    else
        state, policy = oracle_state(o, cfg), o.policy
        t = policy(state)
        next = t(cfg)
        return (state, next)
    end
end

include("oracles/static.jl")
include("oracles/dynamic.jl")


text:
export
    wigner,
    WignerFunction

#=
    Wigner function by Laguerre Polynominal
=#

function wigner(m::Integer, n::Integer, x::Real, p::Real)
    w = gaussian_function(x, p)
    w *= coefficient_of_wave_function(m, n)
    w *= z_to_power(m, n, x, p)
    w *= laguerre(m, n, x, p)

    return w
end

wigner(m::Integer, n::Integer) = (x, p)->wigner(m, n, x, p)

mutable struct WignerFunction{T<:Integer}
    m_dim::T
    n_dim::T
    xs
    ps
    W::Array{ComplexF64,4}

    function WignerFunction(m_dim::T, n_dim::T, xs, ps) where {T<:Integer}
        if m_dim != 0 && n_dim != 0 && !isempty(xs) && !isempty(ps)
            W = Array{ComplexF64,4}(undef, m_dim, n_dim, length(xs), length(ps))
            for m = 1:m_dim, n = 1:n_dim, (x_i, x) = enumerate(xs), (p_j, p) = enumerate(ps)
                W[m, n, x_i, p_j] = wigner(m ,n, x, p)
            end
        else
            W = Array{ComplexF64,4}(undef, 0, 0, 0, 0)
        end
        new{T}(m_dim, n_dim, xs, ps, W)
    end
end

function WignerFunction(m_dim::T, n_dim::T) where {T<:Integer}
    return WignerFunction(m_dim, n_dim, [], [])
end

function WignerFunction(xs::Vector, ps::Vector)
    return WignerFunction(0, 0, xs, ps)
end

function WignerFunction(xs::StepRangeLen, ps::StepRangeLen; ρ_size=35)
    return WignerFunction(ρ_size, ρ_size, xs, ps)
end

function (wf::WignerFunction)(ρ::AbstractMatrix)
    reshape(real(sum(ρ .* wf.W, dims=(1, 2))), length(wf.xs), length(wf.ps))
end

#=
    Wigner function by Fourier transform
=#

function wigner(x, p)
end


text:
using GaussianProcesses
using BenchmarkLite

# Define Benchmark test

type GP_UpdateBenchmark <: Proc
    d::Int
    op::Function
    GP_UpdateBenchmark(d::Int, op::Function) = new(d, op)
end

AbstractString(proc::GP_UpdateBenchmark) = "Dim: $(proc.d)"
Base.length(proc::GP_UpdateBenchmark, cfg) = cfg
Base.isvalid(proc::GP_UpdateBenchmark, cfg) = (isa(cfg, Int) && cfg > 0)
function Base.start(proc::GP_UpdateBenchmark, cfg)
    n = cfg
    x = 2π * rand(proc.d,n)                              # Training set
    y = Float64[sum(sin(x[:,i])) for i in 1:n]/proc.d    # y = 1/d Σᵢ sin(xᵢ)
    GP(x, y, MeanZero(), SEIso(0.0,0.0))
end

function Base.run(proc::GP_UpdateBenchmark, cfg, s)
    proc.op(s)
end
Base.done(proc::GP_UpdateBenchmark, cfg, s) = nothing

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

# Inputs parameters
dims = [10, 30, 50, 100] # Dimensions
num_train = [20, 50, 100, 1000]  # Num training points

println("Benchmarking update_mll! function....")
mll_procs = vec(Proc[ GP_UpdateBenchmark(d, GaussianProcesses.update_mll_and_dmll!) for d in dims ])
mll_table = run(mll_procs, num_train)
show(mll_table)


text:
export InteractiveBackend

# *************************************************************************** #
#                             Interactive Backend
# *************************************************************************** #

"The Interactive backend gets functions to act on a TensorNetworkCircuit
interactively."
mutable struct InteractiveBackend{T} <: AbstractBackend
    tensors::Dict{Symbol, T}
    metrics::Metrics

    function InteractiveBackend{T}() where {T}
        tensors = Dict{Symbol, T}()
        new(tensors, Metrics())
    end
    # Default construct with ComplexF32
    function InteractiveBackend()
        tensors = Dict{Symbol, Array{ComplexF32}}()
        return InteractiveBackend{Array{ComplexF32}}()
    end
end

"""
    save_tensor_data(backend::InteractiveBackend{T},
                              tensor_label::Symbol,
                              tensor_data::Array{<:Number}) where T <: AbstractArray

Function to save tensor data to a dictionary in an interactive backend.
"""
function save_tensor_data(backend::InteractiveBackend{T},
                          tensor_label::Symbol,
                          tensor_data::U) where {T <: AbstractArray, U <: AbstractArray}
    backend.tensors[tensor_label] = convert(T, tensor_data)
end

"""
    load_tensor_data(backend::InteractiveBackend{T},
                              tensor_label::Symbol) where T <: AbstractArray

Function to load tensor data from backend storage (if present)
"""
function load_tensor_data(backend::InteractiveBackend{T},
                          tensor_label::Symbol) where T <: AbstractArray
    if tensor_label in keys(backend.tensors)
        return backend.tensors[tensor_label]
    end
end

"""
    contract_tensors(backend::InteractiveBackend{T},
                              A_label::Symbol, A_ncon_indices::Array{Int, 1},
                              B_label::Symbol, B_ncon_indices::Array{Int, 1},
                              C_label::Symbol) where T <: AbstractArray

Function to interactively contract two tensors A and B to create a tensor C.
"""
function contract_tensors(backend::InteractiveBackend{T},
                          A_label::Symbol, A_ncon_indices::Array{Int, 1},
                          B_label::Symbol, B_ncon_indices::Array{Int, 1},
                          C_label::Symbol) where T <: AbstractArray

    # Get the tensor data for A and B and compute the tensor data for C.
    A_data = backend.tensors[A_label]
    B_data = backend.tensors[B_label]
    C_data = contract_tensors((A_data, B_data),
                              (A_ncon_indices, B_ncon_indices))

    # Save the new tensor and delete the tensors that were contracted.
    save_tensor_data(backend, C_label, C_data)
    delete_tensor!(backend, A_label)
    delete_tensor!(backend, B_label)
end

"""
    save_output(backend::InteractiveBackend{T},
                         node::Symbol, name::String) where T <: AbstractArray

Function to save the result of fully contracting a network under the given
name.
"""
function save_output(backend::InteractiveBackend{T},
                     node::Symbol, name::String="result") where T <: AbstractArray
    name = Symbol(name)
    save_tensor_data(backend, name, backend.tensors[node])
end

"""
    reshape_tensor(backend::InteractiveBackend{T},
                            tensor::Symbol,
                            groups::Array{Array{Int, 1}, 1}) where T <: AbstractArray

Function to reshape a given tensor.
"""
function reshape_tensor(backend::InteractiveBackend{T},
                        tensor::Symbol,
                        groups::Array{Array{Int, 1}, 1}) where T <: AbstractArray
    tensor_dims = size(backend.tensors[tensor])
    backend.tensors[tensor] = reshape_tensor(backend.tensors[tensor], [prod(tensor_dims[x]) for x in groups])
end

"""
    permute_tensor(backend::InteractiveBackend{T},
                            tensor::Symbol, axes::Array{Int, 1}) where T <: AbstractArray

Function to permute the axes of the given tensor
"""
function permute_tensor(backend::InteractiveBackend{T},
                        tensor::Symbol,
                        axes::Array{Int, 1}) where T <: AbstractArray
    backend.tensors[tensor] = permute_tensor(backend.tensors[tensor], axes)
end

"""
    decompose_tensor!(backend::InteractiveBackend{T},
                               tensor::Symbol,
                               left_positions::Array{Int, 1},
                               right_positions::Array{Int, 1};
                               threshold::AbstractFloat=1e-13,
                               max_rank::Int=0,
                               left_label::Symbol,
                               right_label::Symbol) where T <: AbstractArray

Function to decompose a single tensor into two tensors and return the dimension
of the newly created virtual edge.
"""
function decompose_tensor!(backend::InteractiveBackend{T},
                           tensor::Symbol,
                           left_positions::Array{Int, 1},
                           right_positions::Array{Int, 1};
                           threshold::AbstractFloat=1e-13,
                           max_rank::Int=0,
                           left_label::Symbol,
                           right_label::Symbol) where T <: AbstractArray

    node_data = backend.tensors[tensor]

    (B, C, chi) = decompose_tensor(node_data,
                              left_positions,
                              right_positions,
                              threshold=threshold,
                              max_rank=max_rank)

    backend.tensors[left_label] = B
    backend.tensors[right_label] = C

    delete_tensor!(backend, tensor)
    chi
end

"""
    delete_tensor!(backend::InteractiveBackend{T}, tensor_label::Symbol) where T <: AbstractArray

Mark tensor for deletion
"""
function delete_tensor!(backend::InteractiveBackend{T}, tensor_label::Symbol) where T <: AbstractArray
    delete!(backend.tensors, tensor_label)
end

"""
    view_tensor!(backend::InteractiveBackend{T}, view_node, node, bond_idx, bond_range) where T <: AbstractArray

Create a view on a tensor
"""
function view_tensor!(backend::InteractiveBackend{T}, view_node, node, bond_idx, bond_range) where T <: AbstractArray
    node_data = backend.tensors[node]
    backend.tensors[view_node] = tensor_view(node_data, bond_idx, bond_range)
end


text:
"""
    Replace(; value)

Replace `missing`s with one of the specified constant values, depending on the input type.
If multiple values of the same type are provided then the first one will be used.
If the input data is of a different type then the no replacement will be performed.

# Keyword Arguments
* `values::Tuple`: A scalar or tuple of different values that should be used to replace
  missings. Typically, one value per type you're considering imputing for.

# Example
```jldoctest
julia> using Impute: Replace, impute

julia> M = [1.0 2.0 missing missing 5.0; 1.1 2.2 3.3 missing 5.5]
2×5 Matrix{Union{Missing, Float64}}:
 1.0  2.0   missing  missing  5.0
 1.1  2.2  3.3       missing  5.5

julia> impute(M, Replace(; values=0.0); dims=2)
2×5 Matrix{Union{Missing, Float64}}:
 1.0  2.0  0.0  0.0  5.0
 1.1  2.2  3.3  0.0  5.5
```
"""
struct Replace <: Imputor
    values::Tuple
end

Replace(; values) = isa(values, Tuple) ? Replace(values) : Replace(tuple(values))

function _impute!(data::AbstractArray{Union{T, Missing}}, imp::Replace) where T
    i = findfirst(x -> isa(x, T), imp.values)
    i === nothing && return data
    return Base.replace!(data, missing => imp.values[i])
end


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

mutable struct VirtualMachineScaleSetExtensionListResult <: SwaggerModel
    value::Any # spec type: Union{ Nothing, Vector{VirtualMachineScaleSetExtension} } # spec name: value
    nextLink::Any # spec type: Union{ Nothing, String } # spec name: nextLink

    function VirtualMachineScaleSetExtensionListResult(;value=nothing, nextLink=nothing)
        o = new()
        validate_property(VirtualMachineScaleSetExtensionListResult, Symbol("value"), value)
        setfield!(o, Symbol("value"), value)
        validate_property(VirtualMachineScaleSetExtensionListResult, Symbol("nextLink"), nextLink)
        setfield!(o, Symbol("nextLink"), nextLink)
        o
    end
end # type VirtualMachineScaleSetExtensionListResult

const _property_map_VirtualMachineScaleSetExtensionListResult = Dict{Symbol,Symbol}(Symbol("value")=>Symbol("value"), Symbol("nextLink")=>Symbol("nextLink"))
const _property_types_VirtualMachineScaleSetExtensionListResult = Dict{Symbol,String}(Symbol("value")=>"Vector{VirtualMachineScaleSetExtension}", Symbol("nextLink")=>"String")
Base.propertynames(::Type{ VirtualMachineScaleSetExtensionListResult }) = collect(keys(_property_map_VirtualMachineScaleSetExtensionListResult))
Swagger.property_type(::Type{ VirtualMachineScaleSetExtensionListResult }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_VirtualMachineScaleSetExtensionListResult[name]))}
Swagger.field_name(::Type{ VirtualMachineScaleSetExtensionListResult }, property_name::Symbol) =  _property_map_VirtualMachineScaleSetExtensionListResult[property_name]

function check_required(o::VirtualMachineScaleSetExtensionListResult)
    (getproperty(o, Symbol("value")) === nothing) && (return false)
    true
end

function validate_property(::Type{ VirtualMachineScaleSetExtensionListResult }, name::Symbol, val)
end


text:
#ДАНО: Где-то на неограниченном со всех сторон поле и без внутренних перегородок имеется единственный маркер.
#Робот - в произвольной клетке поля.
#РЕЗУЛЬТАТ: Робот - в клетке с тем маркером.

function glavnaia(r::Robot)
    num_steps=1
    side=Nord
    while ismarker(r)==false
        for _ in 1:2
            check_mark(r,side,num_steps)
            side=next(side)
        end
        num_steps+=1
    end
end

function check_mark(r::Robot, side::HorizonSide, num_steps::Int)
    for _ in 1:num_steps
        move!(r,side)
        if ismarker(r)==true
            break
        end
    end
end

next(side::HorizonSide)=HorizonSide(mod(Int(side)+3,4))

max_stars_repo_name: Aminulik/mirea-progs

max_issues_repo_name: Aminulik/mirea-progs

max_forks_repo_name: Aminulik/mirea-progs


text:
# * Inner products

_inner_product(a, B::AbstractFiniteDifferences, b) = step(B)*(a*b)
_inner_product(a,  ::FEDVR, b) = (a*b)

# u or v is a matrix and the result should be an array
function _inner_product(u, Ac, B, v)
    A = parent(Ac)
    lsel,rsel = combined_restriction_selection(A, B)
    _inner_product(view(u, :, lsel), B, view(v, rsel, :))
end

# Left and right functions are vectors and their inner product should
# yield a scalar
function _inner_product(u::Adjoint{<:Any,<:AbstractVector},
                        Ac, B,
                        v::AbstractVector)
    A = parent(Ac)
    lsel,rsel = combined_restriction_selection(A, B)
    _inner_product(transpose(view(u, lsel)), B, view(v, rsel))
end

# vector'vector -> scalar
function _inner_product(u::Adjoint{<:Any, <:AbstractVector},
                        Ac::QuasiAdjoint{T,BB}, B::BB,
                        v::AbstractVector) where {T,BB<:BSpline{T}}
    A = parent(Ac)
    S = A === B ? B.S : (Ac*B)

    dot(u', S, v)
end

# vector'vector -> scalar
_inner_product(u::Adjoint{<:Any, <:AbstractVector},
               Ac::AdjointBSplineOrRestricted, B::BSplineOrRestricted,
               v::AbstractVector) =
    dot(u', Ac*B, v)

# u or v is a matrix and the result should be an array
function _inner_product(u::Adjoint{<:Any, <:AbstractVecOrMat},
                        Ac::QuasiAdjoint{T,BB}, B::BB,
                        v::AbstractVecOrMat) where {T,BB<:BSpline{T}}
    A = parent(Ac)
    S = A === B ? B.S : (Ac*B)

    u*(S*v)
end

const FlatFuncInnerProduct{T,B} = Tuple{<:Adjoint{T,<:AbstractVecOrMat{T}},
                                        <:AdjointBasisOrRestricted{B},
                                        <:BasisOrRestricted{B},
                                        <:AbstractVecOrMat{T}}

materialize(M::Mul{<:Any, <:FlatFuncInnerProduct}) = _inner_product(M.args...)

materialize(M::Mul{<:Any, Tuple{AdjointFuncVector{<:Any,BB},
                                FuncVector{<:Any,BB}}}) where BB =
                                    _inner_product(M.args[1].args..., M.args[2].args...)

# * Norms

const FDOrFEDVRVecOrMat = Union{FDVecOrMat,
                                Mul{<:Any, <:Tuple{<:FiniteDifferencesOrRestricted, <:AbstractArray}},
                                FEDVRVecOrMat,
                                Mul{<:Any, <:Tuple{<:FEDVROrRestricted, <:AbstractArray}}}

_norm(B::FiniteDifferencesOrRestricted, c::AbstractArray, p::Real=2) =
    norm(c, p)*(step(B)^(inv(p)))

_norm(R::FEDVROrRestricted, ϕ::AbstractArray, p::Real=2) = norm(ϕ, p)

LinearAlgebra.norm(v::FDOrFEDVRVecOrMat, p::Real=2) = _norm(v.args..., p)

function LinearAlgebra.normalize!(v::FDOrFEDVRVecOrMat, p::Real=2)
    v.args[2][:] /= norm(v,p)
    v
end

LinearAlgebra.norm(v::SplineVector) = √(v'v)
function LinearAlgebra.normalize!(v::SplineVector)
    v.args[2][:] /= norm(v)
    v
end