text:
using Mocking

create_device_pool_1_patch = @patch function create_device_pool(a...; b...)
    return Dict{String, Any}("devicePool" => Dict{String, Any}())
end

create_project_1_patch = @patch function create_project(a...; b...)
    return Dict{String, Any}("project" => Dict{String, Any}("name" => "MyProject", "created" => "1472660939.152", "arn" => "arn:aws:devicefarm:us-west-2:123456789101:project:5e01a8c7-c861-4c0a-b1d5-12345EXAMPLE"))
end

create_remote_access_session_1_patch = @patch function create_remote_access_session(a...; b...)
    return Dict{String, Any}("remoteAccessSession" => Dict{String, Any}())
end

create_upload_1_patch = @patch function create_upload(a...; b...)
    return Dict{String, Any}("upload" => Dict{String, Any}("name" => "MyAppiumPythonUpload", "status" => "INITIALIZED", "created" => "1472661404.186", "type" => "APPIUM_PYTHON_TEST_PACKAGE", "arn" => "arn:aws:devicefarm:us-west-2:123456789101:upload:5e01a8c7-c861-4c0a-b1d5-5ec6e6c6dd23/b5340a65-3da7-4da6-a26e-12345EXAMPLE", "url" => "https://prod-us-west-2-uploads.s3-us-west-2.amazonaws.com/arn%3Aaws%3Adevicefarm%3Aus-west-2%3A123456789101%3Aproject%3A5e01a8c7-c861-4c0a-b1d5-12345EXAMPLE/uploads/arn%3Aaws%3Adevicefarm%3Aus-west-2%3A123456789101%3Aupload%3A5e01a8c7-c861-4c0a-b1d5-5ec6e6c6dd23/b5340a65-3da7-4da6-a26e-12345EXAMPLE/MyAppiumPythonUpload?AWSAccessKeyId=1234567891011EXAMPLE&Expires=1472747804&Signature=1234567891011EXAMPLE"))
end

delete_device_pool_1_patch = @patch function delete_device_pool(a...; b...)
    return Dict{String, Any}()
end

delete_project_1_patch = @patch function delete_project(a...; b...)
    return Dict{String, Any}()
end

delete_remote_access_session_1_patch = @patch function delete_remote_access_session(a...; b...)
    return Dict{String, Any}()
end

delete_run_1_patch = @patch function delete_run(a...; b...)
    return Dict{String, Any}()
end

delete_upload_1_patch = @patch function delete_upload(a...; b...)
    return Dict{String, Any}()
end

get_account_settings_1_patch = @patch function get_account_settings(a...; b...)
    return Dict{String, Any}("accountSettings" => Dict{String, Any}("awsAccountNumber" => "123456789101", "unmeteredDevices" => Dict{String, Any}("ANDROID" => 1, "IOS" => 2)))
end

get_device_1_patch = @patch function get_device(a...; b...)
    return Dict{String, Any}("device" => Dict{String, Any}("model" => "G2 (Sprint)", "formFactor" => "PHONE", "memory" => 16000000000, "name" => "LG G2 (Sprint)", "cpu" => Dict{String, Any}("architecture" => "armeabi-v7a", "frequency" => "MHz", "clock" => 2265.6), "platform" => "ANDROID", "arn" => "arn:aws:devicefarm:us-west-2::device:A0E6E6E1059E45918208DF75B2B7EF6C", "heapSize" => 256000000, "resolution" => Dict{String, Any}("height" => 1920, "width" => 1080), "image" => "75B2B7EF6C12345EXAMPLE", "manufacturer" => "LG", "os" => "4.2.2"))
end

get_device_pool_1_patch = @patch function get_device_pool(a...; b...)
    return Dict{String, Any}("devicePool" => Dict{String, Any}())
end

get_device_pool_compatibility_1_patch = @patch function get_device_pool_compatibility(a...; b...)
    return Dict{String, Any}("incompatibleDevices" => Any[], "compatibleDevices" => Any[])
end

get_job_1_patch = @patch function get_job(a...; b...)
    return Dict{String, Any}("job" => Dict{String, Any}())
end

get_offering_status_1_patch = @patch function get_offering_status(a...; b...)
    return Dict{String, Any}("nextPeriod" => Dict{String, Any}("D68B3C05-1BA6-4360-BC69-12345EXAMPLE" => Dict{String, Any}("offering" => Dict{String, Any}("platform" => "ANDROID", "id" => "D68B3C05-1BA6-4360-BC69-12345EXAMPLE", "type" => "RECURRING", "description" => "Android Remote Access Unmetered Device Slot"), "effectiveOn" => "1472688000", "quantity" => 1)), "current" => Dict{String, Any}("D68B3C05-1BA6-4360-BC69-12345EXAMPLE" => Dict{String, Any}("offering" => Dict{String, Any}("platform" => "ANDROID", "id" => "D68B3C05-1BA6-4360-BC69-12345EXAMPLE", "type" => "RECURRING", "description" => "Android Remote Access Unmetered Device Slot"), "quantity" => 1)))
end

get_project_1_patch = @patch function get_project(a...; b...)
    return Dict{String, Any}("project" => Dict{String, Any}("name" => "My Project", "created" => "1472660939.152", "arn" => "arn:aws:devicefarm:us-west-2:123456789101:project:5e01a8c7-c861-4c0a-b1d5-12345EXAMPLE"))
end

get_remote_access_session_1_patch = @patch function get_remote_access_session(a...; b...)
    return Dict{String, Any}("remoteAccessSession" => Dict{String, Any}())
end

get_run_1_patch = @patch function get_run(a...; b...)
    return Dict{String, Any}("run" => Dict{String, Any}("totalJobs" => 3, "counters" => Dict{String, Any}("passed" => 0, "total" => 0, "warned" => 0, "stopped" => 0, "failed" => 0, "errored" => 0, "skipped" => 0), "name" => "My Test Run", "status" => "RUNNING", "platform" => "ANDROID", "created" => "1472667509.852", "completedJobs" => 0, "arn" => "arn:aws:devicefarm:us-west-2:123456789101:run:5e01a8c7-c861-4c0a-b1d5-5ec6e6c6dd23/0fcac17b-6122-44d7-ae5a-12345EXAMPLE", "result" => "PENDING", "deviceMinutes" => Dict{String, Any}("total" => 0.0, "metered" => 0.0, "unmetered" => 0.0), "billingMethod" => "METERED", "type" => "BUILTIN_EXPLORER"))
end

get_suite_1_patch = @patch function get_suite(a...; b...)
    return Dict{String, Any}("suite" => Dict{String, Any}())
end

get_test_1_patch = @patch function get_test(a...; b...)
    return Dict{String, Any}("test" => Dict{String, Any}())
end

get_upload_1_patch = @patch function get_upload(a...; b...)
    return Dict{String, Any}("upload" => Dict{String, Any}())
end

install_to_remote_access_session_1_patch = @patch function install_to_remote_access_session(a...; b...)
    return Dict{String, Any}("appUpload" => Dict{String, Any}())
end

list_device_pools_1_patch = @patch function list_device_pools(a...; b...)
    return Dict{String, Any}("devicePools" => Any[Dict{String, Any}("name" => "Top Devices", "rules" => Any[Dict{String, Any}("operator" => "IN", "value" => "[\"arn:aws:devicefarm:us-west-2::device:123456789EXAMPLE\",\"arn:aws:devicefarm:us-west-2::device:123456789EXAMPLE\",\"arn:aws:devicefarm:us-west-2::device:123456789EXAMPLE\",\"arn:aws:devicefarm:us-west-2::device:123456789EXAMPLE\",\"arn:aws:devicefarm:us-west-2::device:123456789EXAMPLE\",\"arn:aws:devicefarm:us-west-2::device:123456789EXAMPLE\",\"arn:aws:devicefarm:us-west-2::device:123456789EXAMPLE\",\"arn:aws:devicefarm:us-west-2::device:123456789EXAMPLE\",\"arn:aws:devicefarm:us-west-2::device:123456789EXAMPLE\",\"arn:aws:devicefarm:us-west-2::device:123456789EXAMPLE\"]", "attribute" => "ARN")], "description" => "Top devices", "arn" => "arn:aws:devicefarm:us-west-2::devicepool:082d10e5-d7d7-48a5-ba5c-12345EXAMPLE"), Dict{String, Any}("name" => "My Android Device Pool", "rules" => Any[Dict{String, Any}("operator" => "IN", "value" => "[\"arn:aws:devicefarm:us-west-2::device:123456789EXAMPLE\",\"arn:aws:devicefarm:us-west-2::device:123456789EXAMPLE\",\"arn:aws:devicefarm:us-west-2::device:123456789EXAMPLE\"]", "attribute" => "ARN")], "description" => "Samsung Galaxy Android devices", "arn" => "arn:aws:devicefarm:us-west-2:123456789101:devicepool:5e01a8c7-c861-4c0a-b1d5-5ec6e6c6dd23/bf96e75a-28f6-4e61-b6a7-12345EXAMPLE")])
end

list_devices_1_patch = @patch function list_devices(a...; b...)
    return Dict{String, Any}()
end

list_offering_transactions_1_patch = @patch function list_offering_transactions(a...; b...)
    return Dict{String, Any}("offeringTransactions" => Any[Dict{String, Any}("cost" => Dict{String, Any}("amount" => 0, "currencyCode" => "USD"), "offeringStatus" => Dict{String, Any}("offering" => Dict{String, Any}("platform" => "ANDROID", "id" => "D68B3C05-1BA6-4360-BC69-12345EXAMPLE", "type" => "RECURRING", "description" => "Android Remote Access Unmetered Device Slot"), "effectiveOn" => "1472688000", "type" => "RENEW", "quantity" => 0), "transactionId" => "03728003-d1ea-4851-abd6-12345EXAMPLE", "createdOn" => "1470021420"), Dict{String, Any}("cost" => Dict{String, Any}("amount" => 250, "currencyCode" => "USD"), "offeringStatus" => Dict{String, Any}("offering" => Dict{String, Any}("platform" => "ANDROID", "id" => "D68B3C05-1BA6-4360-BC69-12345EXAMPLE", "type" => "RECURRING", "description" => "Android Remote Access Unmetered Device Slot"), "effectiveOn" => "1470021420", "type" => "PURCHASE", "quantity" => 1), "transactionId" => "56820b6e-06bd-473a-8ff8-12345EXAMPLE", "createdOn" => "1470021420"), Dict{String, Any}("cost" => Dict{String, Any}("amount" => 175, "currencyCode" => "USD"), "offeringStatus" => Dict{String, Any}("offering" => Dict{String, Any}("platform" => "ANDROID", "id" => "8980F81C-00D7-469D-8EC6-12345EXAMPLE", "type" => "RECURRING", "description" => "Android Unmetered Device Slot"), "effectiveOn" => "1465538520", "type" => "PURCHASE", "quantity" => 1), "transactionId" => "953ae2c6-d760-4a04-9597-12345EXAMPLE", "createdOn" => "1465538520"), Dict{String, Any}("cost" => Dict{String, Any}("amount" => 8.07, "currencyCode" => "USD"), "offeringStatus" => Dict{String, Any}("offering" => Dict{String, Any}("platform" => "IOS", "id" => "A53D4D73-A6F6-4B82-A0B0-12345EXAMPLE", "type" => "RECURRING", "description" => "iOS Unmetered Device Slot"), "effectiveOn" => "1459344300", "type" => "PURCHASE", "quantity" => 1), "transactionId" => "2baf9021-ae3e-47f5-ab52-12345EXAMPLE", "createdOn" => "1459344300")])
end

list_offerings_1_patch = @patch function list_offerings(a...; b...)
    return Dict{String, Any}("offerings" => Any[Dict{String, Any}("platform" => "IOS", "id" => "A53D4D73-A6F6-4B82-A0B0-12345EXAMPLE", "type" => "RECURRING", "description" => "iOS Unmetered Device Slot", "recurringCharges" => Any[Dict{String, Any}("cost" => Dict{String, Any}("amount" => 250, "currencyCode" => "USD"), "frequency" => "MONTHLY")]), Dict{String, Any}("platform" => "ANDROID", "id" => "8980F81C-00D7-469D-8EC6-12345EXAMPLE", "type" => "RECURRING", "description" => "Android Unmetered Device Slot", "recurringCharges" => Any[Dict{String, Any}("cost" => Dict{String, Any}("amount" => 250, "currencyCode" => "USD"), "frequency" => "MONTHLY")]), Dict{String, Any}("platform" => "ANDROID", "id" => "D68B3C05-1BA6-4360-BC69-12345EXAMPLE", "type" => "RECURRING", "description" => "Android Remote Access Unmetered Device Slot", "recurringCharges" => Any[Dict{String, Any}("cost" => Dict{String, Any}("amount" => 250, "currencyCode" => "USD"), "frequency" => "MONTHLY")]), Dict{String, Any}("platform" => "IOS", "id" => "552B4DAD-A6C9-45C4-94FB-12345EXAMPLE", "type" => "RECURRING", "description" => "iOS Remote Access Unmetered Device Slot", "recurringCharges" => Any[Dict{String, Any}("cost" => Dict{String, Any}("amount" => 250, "currencyCode" => "USD"), "frequency" => "MONTHLY")])])
end

list_projects_1_patch = @patch function list_projects(a...; b...)
    return Dict{String, Any}("projects" => Any[Dict{String, Any}("name" => "My Test Project", "created" => "1453163262.105", "arn" => "arn:aws:devicefarm:us-west-2:123456789101:project:7ad300ed-8183-41a7-bf94-12345EXAMPLE"), Dict{String, Any}("name" => "Hello World", "created" => "1470350112.439", "arn" => "arn:aws:devicefarm:us-west-2:123456789101:project:d6b087d9-56db-4e44-b9ec-12345EXAMPLE")])
end

list_remote_access_sessions_1_patch = @patch function list_remote_access_sessions(a...; b...)
    return Dict{String, Any}("remoteAccessSessions" => Any[])
end

list_runs_1_patch = @patch function list_runs(a...; b...)
    return Dict{String, Any}("runs" => Any[Dict{String, Any}("totalJobs" => 3, "counters" => Dict{String, Any}("passed" => 0, "total" => 0, "warned" => 0, "stopped" => 0, "failed" => 0, "errored" => 0, "skipped" => 0), "name" => "My Test Run", "status" => "RUNNING", "platform" => "ANDROID", "created" => "1472667509.852", "completedJobs" => 0, "arn" => "arn:aws:devicefarm:us-west-2:123456789101:run:5e01a8c7-c861-4c0a-b1d5-5ec6e6c6dd23/0fcac17b-6122-44d7-ae5a-12345EXAMPLE", "result" => "PENDING", "deviceMinutes" => Dict{String, Any}("total" => 0.0, "metered" => 0.0, "unmetered" => 0.0), "billingMethod" => "METERED", "type" => "BUILTIN_EXPLORER")])
end

list_samples_1_patch = @patch function list_samples(a...; b...)
    return Dict{String, Any}("samples" => Any[])
end

list_suites_1_patch = @patch function list_suites(a...; b...)
    return Dict{String, Any}("suites" => Any[])
end

list_tests_1_patch = @patch function list_tests(a...; b...)
    return Dict{String, Any}("tests" => Any[])
end

list_unique_problems_1_patch = @patch function list_unique_problems(a...; b...)
    return Dict{String, Any}("uniqueProblems" => Dict{String, Any}())
end

list_uploads_1_patch = @patch function list_uploads(a...; b...)
    return Dict{String, Any}("uploads" => Any[])
end

purchase_offering_1_patch = @patch function purchase_offering(a...; b...)
    return Dict{String, Any}("offeringTransaction" => Dict{String, Any}("cost" => Dict{String, Any}("amount" => 8.07, "currencyCode" => "USD"), "offeringStatus" => Dict{String, Any}("offering" => Dict{String, Any}("platform" => "ANDROID", "id" => "D68B3C05-1BA6-4360-BC69-12345EXAMPLE", "type" => "RECURRING", "description" => "Android Remote Access Unmetered Device Slot"), "effectiveOn" => "1472648340", "type" => "PURCHASE", "quantity" => 1), "transactionId" => "d30614ed-1b03-404c-9893-12345EXAMPLE", "createdOn" => "1472648340"))
end

renew_offering_1_patch = @patch function renew_offering(a...; b...)
    return Dict{String, Any}("offeringTransaction" => Dict{String, Any}("cost" => Dict{String, Any}("amount" => 250, "currencyCode" => "USD"), "offeringStatus" => Dict{String, Any}("offering" => Dict{String, Any}("platform" => "ANDROID", "id" => "D68B3C05-1BA6-4360-BC69-12345EXAMPLE", "type" => "RECURRING", "description" => "Android Remote Access Unmetered Device Slot"), "effectiveOn" => "1472688000", "type" => "RENEW", "quantity" => 1), "transactionId" => "e90f1405-8c35-4561-be43-12345EXAMPLE", "createdOn" => "1472648880"))
end

schedule_run_1_patch = @patch function schedule_run(a...; b...)
    return Dict{String, Any}("run" => Dict{String, Any}())
end

stop_run_1_patch = @patch function stop_run(a...; b...)
    return Dict{String, Any}("run" => Dict{String, Any}())
end

update_device_pool_1_patch = @patch function update_device_pool(a...; b...)
    return Dict{String, Any}("devicePool" => Dict{String, Any}())
end

update_project_1_patch = @patch function update_project(a...; b...)
    return Dict{String, Any}("project" => Dict{String, Any}("name" => "NewName", "created" => "1448400709.927", "arn" => "arn:aws:devicefarm:us-west-2:123456789101:project:8f75187d-101e-4625-accc-12345EXAMPLE"))
end


text:
# Descriptive Statistics

#############################
#
#   Location
#
#############################

# Geometric mean
function geomean(a::RealArray)
    s = 0.0
    n = length(a)
    for i = 1 : n
        @inbounds s += log(a[i])
    end
    return exp(s / n)
end

# Harmonic mean
function harmmean(a::RealArray)
    s = 0.0
    n = length(a)
    for i in 1 : n
        @inbounds s += inv(a[i])
    end
    return n / s
end

# Trimmed mean
function trimmean(x::RealArray, p::Real)
    n = length(x)
    n > 0 || error("x can not be empty.")
    0 <= p < 1 || error("p must be non-negative and less than 1.")
    rn = min(iround(n * p), n-1)

    sx = sort(x)
    nl = rn >> 1
    nh = (rn - nl)
    s = 0.0
    for i = (1+nl) : (n-nh)
        @inbounds s += x[i]
    end
    return s / (n - rn)
end

# compute mode, given the range of integer values
function mode{T<:Integer}(a::AbstractArray{T}, rgn::UnitRange{T})
    isempty(a) && error("mode: input array cannot be empty.")
    len = length(a)
    r0 = rgn[1]
    r1 = rgn[end]
    cnts = zeros(Int, length(rgn))
    mc = 0    # maximum count
    mv = r0   # a value corresponding to maximum count
    for i = 1:len
        @inbounds x = a[i]
        if r0 <= x <= r1
            @inbounds c = (cnts[x - r0 + 1] += 1)
            if c > mc
                mc = c
                mv = x
            end
        end
    end
    return mv
end

function modes{T<:Integer}(a::AbstractArray{T}, rgn::UnitRange{T})
    r0 = rgn[1]
    r1 = rgn[end]
    n = length(rgn)
    cnts = zeros(Int, n)
    # find the maximum count
    mc = 0
    for i = 1:length(a)
        @inbounds x = a[i]
        if r0 <= x <= r1
            @inbounds c = (cnts[x - r0 + 1] += 1)
            if c > mc
                mc = c
            end
        end
    end
    # find all values corresponding to maximum count
    ms = T[]
    for i = 1:n
        @inbounds if cnts[i] == mc
            push!(ms, rgn[i])
        end
    end
    return ms
end

# compute mode over arbitrary array
function mode{T}(a::AbstractArray{T})
    isempty(a) && error("mode: input array cannot be empty.")
    cnts = Dict{T,Int}()
    # first element
    mc = 1
    mv = a[1]
    cnts[mv] = 1
    # find the mode along with table construction
    for i = 2 : length(a)
        @inbounds x = a[i]
        if haskey(cnts, x)
            c = (cnts[x] += 1)
            if c > mc
                mc = c
                mv = x
            end
        else
            cnts[x] = 1
            # in this case: c = 1, and thus c > mc won't happen
        end
    end
    return mv
end

function modes{T}(a::AbstractArray{T})
    isempty(a) && error("modes: input array cannot be empty.")
    cnts = Dict{T,Int}()
    # first element
    mc = 1
    cnts[a[1]] = 1
    # find the mode along with table construction
    for i = 2 : length(a)
        @inbounds x = a[i]
        if haskey(cnts, x)
            c = (cnts[x] += 1)
            if c > mc
                mc = c
            end
        else
            cnts[x] = 1
            # in this case: c = 1, and thus c > mc won't happen
        end
    end
    # find values corresponding to maximum counts
    ms = T[]
    for (x, c) in cnts
        if c == mc
            push!(ms, x)
        end
    end
    return ms
end

#############################
#
#   quantile and friends
#
#############################

percentile{T<:Real}(v::AbstractArray{T}, p) = quantile(v, p * 0.01)

quantile{T<:Real}(v::AbstractArray{T}) = quantile(v, [.0, .25, .5, .75, 1.0])
nquantile{T<:Real}(v::AbstractArray{T}, n::Integer) = quantile(v, (0:n)/n)

#############################
#
#   Dispersion
#
#############################

# Variation coefficient: std / mean
variation{T<:Real}(x::AbstractArray{T}, m::Real) = stdm(x, m) / m
variation{T<:Real}(x::AbstractArray{T}) = variation(x, mean(x))

# Standard error of the mean: std(a) / sqrt(len)
sem{T<:Real}(a::AbstractArray{T}) = sqrt(var(a) / length(a))

# Median absolute deviation
mad{T<:Real}(v::AbstractArray{T}, args...;arg...) = mad!(copy(v), args...;arg...)
mad{T<:Real}(v::Range{T}, args...;arg...) = mad!([v], args...;arg...)

function mad!{T<:Real}(v::AbstractArray{T}, center::Real=median!(v); constant::Real=1.4826)
    for i in 1:length(v)
        @inbounds v[i] = abs(v[i]-center)
    end
    constant * median!(v)
end

# Interquartile range
iqr{T<:Real}(v::AbstractArray{T}) = (q = quantile(v, [.25, .75]); q[2] - q[1])

#############################
#
#   Z-scores
#
#############################

function _zscore!(Z::AbstractArray, X::AbstractArray, μ::Real, σ::Real)
    # Z and X are assumed to have the same size
    iσ = inv(σ)
    if μ == zero(μ)
        for i = 1 : length(X)
            @inbounds Z[i] = X[i] * iσ
        end
    else
        for i = 1 : length(X)
            @inbounds Z[i] = (X[i] - μ) * iσ
        end
    end
    return Z
end

@ngenerate N typeof(Z) function _zscore!{S,T,N}(Z::AbstractArray{S,N}, X::AbstractArray{T,N}, μ::AbstractArray, σ::AbstractArray)
    # Z and X are assumed to have the same size
    # μ and σ are assumed to have the same size, that is compatible with size(X)
    siz1 = size(X, 1)
    @nextract N ud d->size(μ, d)
    if size(μ, 1) == 1 && siz1 > 1
        @nloops N i d->(d>1 ? (1:size(X,d)) : (1:1)) d->(j_d = ud_d ==1 ? 1 : i_d) begin
            v = (@nref N μ j)
            c = inv(@nref N σ j)
            for i_1 = 1:siz1
                (@nref N Z i) = ((@nref N X i) - v) * c
            end
        end
    else
        @nloops N i X d->(j_d = ud_d ==1 ? 1 : i_d) begin
            (@nref N Z i) = ((@nref N X i) - (@nref N μ j)) / (@nref N σ j)
        end
    end
    return Z
end

function _zscore_chksize(X::AbstractArray, μ::AbstractArray, σ::AbstractArray)
    size(μ) == size(σ) || throw(DimensionMismatch("μ and σ should have the same size."))
    for i=1:ndims(X)
        dμ_i = size(μ,i)
        (dμ_i == 1 || dμ_i == size(X,i)) || throw(DimensionMismatch("X and μ have incompatible sizes."))
    end
end

function zscore!{ZT<:FloatingPoint,T<:Real}(Z::AbstractArray{ZT}, X::AbstractArray{T}, μ::Real, σ::Real)
    size(Z) == size(X) || throw(DimensionMismatch("Z and X must have the same size."))
    _zscore!(Z, X, μ, σ)
end

function zscore!{ZT<:FloatingPoint,T<:Real,U<:Real,S<:Real}(Z::AbstractArray{ZT}, X::AbstractArray{T},
                                                            μ::AbstractArray{U}, σ::AbstractArray{S})
    size(Z) == size(X) || throw(DimensionMismatch("Z and X must have the same size."))
    _zscore_chksize(X, μ, σ)
    _zscore!(Z, X, μ, σ)
end

zscore!{T<:FloatingPoint}(X::AbstractArray{T}, μ::Real, σ::Real) = _zscore!(X, X, μ, σ)

zscore!{T<:FloatingPoint,U<:Real,S<:Real}(X::AbstractArray{T}, μ::AbstractArray{U}, σ::AbstractArray{S}) =
    (_zscore_chksize(X, μ, σ); _zscore!(X, X, μ, σ))

function zscore{T<:Real}(X::AbstractArray{T}, μ::Real, σ::Real)
    ZT = typeof((zero(T) - zero(μ)) / one(σ))
    _zscore!(Array(ZT, size(X)), X, μ, σ)
end

function zscore{T<:Real,U<:Real,S<:Real}(X::AbstractArray{T}, μ::AbstractArray{U}, σ::AbstractArray{S})
    _zscore_chksize(X, μ, σ)
    ZT = typeof((zero(T) - zero(U)) / one(S))
    _zscore!(Array(ZT, size(X)), X, μ, σ)
end

zscore{T<:Real}(X::AbstractArray{T}) = ((μ, σ) = mean_and_std(X); zscore(X, μ, σ))
zscore{T<:Real}(X::AbstractArray{T}, dim::Int) = ((μ, σ) = mean_and_std(X, dim); zscore(X, μ, σ))

#############################
#
#   entropy and friends
#
#############################

function entropy{T<:Real}(p::AbstractArray{T})
    s = 0.
    z = zero(T)
    for i = 1:length(p)
        @inbounds pi = p[i]
        if pi > z
            s += pi * log(pi)
        end
    end
    return -s
end

function crossentropy{T<:Real}(p::AbstractArray{T}, q::AbstractArray{T})
    length(p) == length(q) || throw(DimensionMismatch("Inconsistent array length."))
    s = 0.
    z = zero(T)
    for i = 1:length(p)
        @inbounds pi = p[i]
        @inbounds qi = q[i]
        if pi > z
            s += pi * log(qi)
        end
    end
    return -s
end

function kldivergence{T<:Real}(p::AbstractArray{T}, q::AbstractArray{T})
    length(p) == length(q) || throw(DimensionMismatch("Inconsistent array length."))
    s = 0.
    z = zero(T)
    for i = 1:length(p)
        @inbounds pi = p[i]
        @inbounds qi = q[i]
        if pi > z
            s += pi * log(pi / qi)
        end
    end
    return s
end

#############################
#
#   summary
#
#############################

immutable SummaryStats{T<:FloatingPoint}
    mean::T
    min::T
    q25::T
    median::T
    q75::T
    max::T
end

function summarystats{T<:Real}(a::AbstractArray{T})
    m = mean(a)
    qs = quantile(a, [0.00, 0.25, 0.50, 0.75, 1.00])
    R = typeof(convert(FloatingPoint, zero(T)))
    SummaryStats{R}(
        convert(R, m),
        convert(R, qs[1]),
        convert(R, qs[2]),
        convert(R, qs[3]),
        convert(R, qs[4]),
        convert(R, qs[5]))
end

function Base.show(io::IO, ss::SummaryStats)
    println(io, "Summary Stats:")
    @printf(io, "Mean:         %.6f\n", ss.mean)
    @printf(io, "Minimum:      %.6f\n", ss.min)
    @printf(io, "1st Quartile: %.6f\n", ss.q25)
    @printf(io, "Median:       %.6f\n", ss.median)
    @printf(io, "3rd Quartile: %.6f\n", ss.q75)
    @printf(io, "Maximum:      %.6f\n", ss.max)
end

describe{T<:Real}(a::AbstractArray{T}) = show(summarystats(a))

ext: jl


text:
using Test
using Plots

using jInvVis
using jInv.Mesh

# tests for regular mesh
domain = [1 2 0 .1 3 4]
n      = [8 5 3]
M2D    = getRegularMesh(domain[1:4],n[1:2])
M3D    = getRegularMesh(domain,n)

println("=== test plotGrid (regular mesh) ===")
plotGrid(M2D)
# plotGrid(M3D)
# plotGrid(getCellCenteredGrid(M2D),M2D);
println("==== all passed === ")


text:
using Distributions: ValueSupport, VariateForm

struct MixedSupport <: ValueSupport end
struct MixedVariate <: VariateForm end

struct JointDistribution{A0,A,B,M} <: Distribution{MixedVariate, MixedSupport}
    model::Model{A,B,M}
    args::A0
end

(jd::JointDistribution)(nt::NamedTuple) = JointDistribution(jd.model, merge(jd.args, nt))

(m::Model)(;args...)= JointDistribution(m,(;args...))

(m::Model)(nt::NamedTuple) = JointDistribution(m,nt)

function Base.show(io::IO, d :: JointDistribution)
    m = d.model
    println(io, "Joint Distribution")
    print(io, "    Bound arguments: [")
    join(io, fieldnames(arguments(d)), ", ")
    println(io, "]")
    print(io, "    Variables: [")
    join(io, setdiff(toposort(m),arguments(m)), ", ")
    println(io, "]\n")
    println(io, convert(Expr, m))
end


text:
function potential_energy(A::AbstractSystem, p::EmpiricalPotential)
    e = 0.0
    nnlist = neighborlist(A, p.rcutoff)

    for ii in 1:length(A)
        for R in nnlist.R[ii]
            e += potential_energy(R, p)
        end
    end
    e
end


text:
module AppleAccelerateLinAlgWrapper

using libblastrampoline_jll, Libdl
using LinearAlgebra: BLAS, Adjoint, Transpose,
    UnitLowerTriangular, LowerTriangular, UnitUpperTriangular, UpperTriangular, AbstractTriangular
import LinearAlgebra

if Sys.isapple()
    default_blas_int() = Val(Int32)
else
    default_blas_int() = Val(LinearAlgebra.BlasInt)
end

const LIBBLASTRAMPOLINE_HANDLE = Ref{Ptr{Cvoid}}(C_NULL)

for f ∈ [:gemm, :getrf, :trsm]
    get = Symbol(:get_, f, :_addr)
    for (p, T) ∈ ((:s, :Float32), (:d, :Float64))
        pf = Symbol(p, f, :_)
        pf64 = Symbol(pf, Symbol("64_"))
        @eval begin
            function $get(::Val{$T}, ::Val{Int32})
                dlsym(LIBBLASTRAMPOLINE_HANDLE[], $(QuoteNode(pf)))::Ptr{Cvoid}
            end
            function $get(::Val{$T}, ::Val{Int64})
                dlsym(LIBBLASTRAMPOLINE_HANDLE[], $(QuoteNode(pf64)))::Ptr{Cvoid}
            end
        end
    end
end

function gemm!(transA::AbstractChar, transB::AbstractChar,
               alpha::Union{T, Bool},
               A::StridedMatrix{T}, B::StridedMatrix{T},
               beta::Union{T, Bool},
               C::StridedMatrix{T}, ::Val{BlasInt}) where {BlasInt <: Union{Int32,Int64},T<:Union{Float32,Float64}}
    m = size(A, transA == 'N' ? 1 : 2) % BlasInt
    ka = size(A, transA == 'N' ? 2 : 1) % BlasInt
    kb = size(B, transB == 'N' ? 1 : 2) % BlasInt
    n = size(B, transB == 'N' ? 2 : 1) % BlasInt
    if ka != kb || m != (size(C,1) % BlasInt) || n != (size(C,2) % BlasInt)
        throw(DimensionMismatch("A has size ($m,$ka), B has size ($kb,$n), C has size $(size(C))"))
    end
    fptr = get_gemm_addr(Val(T), Val(BlasInt))
#    ccall(:jl_breakpoint, Cvoid, (Any,), fptr)
    ccall(fptr, Cvoid, (Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt},
                        Ref{BlasInt}, Ref{T}, Ptr{T}, Ref{BlasInt},
                        Ptr{T}, Ref{BlasInt}, Ref{T}, Ptr{T},
                        Ref{BlasInt}, Clong, Clong),
          transA, transB, m, n,
          ka, alpha, A, max(1,stride(A,2)),
          B, max(1,stride(B,2)), beta, C,
          max(1,stride(C,2)) % BlasInt, 1, 1)
    C
end
function gemm(
    transA::AbstractChar, transB::AbstractChar, alpha::T, A::StridedMatrix{T}, B::StridedMatrix{T}, ::Val{BlasInt}
) where {BlasInt <: Union{Int32,Int64}, T <: Union{Float32,Float64}}
    gemm!(transA, transB, alpha, A, B, zero(T), similar(B, T, (size(A, transA == 'N' ? 1 : 2), size(B, transB == 'N' ? 2 : 1))), Val(BlasInt))
end

untranspose_flag(A::Adjoint, f::Bool = false) = untranspose_flag(parent(A), !f)
untranspose_flag(A::Transpose, f::Bool = false) = untranspose_flag(parent(A), !f)
untranspose_flag(A, f::Bool = false) = (A, f)

function gemm!(
    C::AbstractMatrix{T}, A::AbstractMatrix{T}, B::AbstractMatrix{T}, α = one(T), β = zero(T), ::Val{BlasInt} = default_blas_int()
) where {T<:Union{Float32,Float64},BlasInt}
    pA, fa = untranspose_flag(A)
    pB, fb = untranspose_flag(B)
    gemm!(ifelse(fa, 'T', 'N'), ifelse(fb, 'T', 'N'), α, parent(A), parent(B), β, C, Val(BlasInt))
end
function gemm(
    A::AbstractMatrix{T}, B::AbstractMatrix{T}, α = one(T), β = zero(T), ::Val{BlasInt} = default_blas_int()
) where {T<:Union{Float32,Float64},BlasInt<:Union{Int32,Int64}}
    gemm!(similar(A, (size(A,1),size(B,2))), A, B, α, β, Val(BlasInt))
end

struct LU{T,M<:StridedMatrix{T},BlasInt}
    factors::M
    ipiv::Vector{BlasInt}
    info::BlasInt
end

function getrf_ipiv(A::StridedMatrix, ::Val{BlasInt}) where {BlasInt}
    M,N = size(A)
    Vector{BlasInt}(undef, min(M,N))
end
function getrf!(
    A::StridedMatrix{T}, ipiv::StridedVector{BlasInt} = getrf_ipiv(A, default_blas_int())
) where {T <: Union{Float32,Float64}, BlasInt <: Union{Int32,Int64}}
    M, N = size(A)
    fptr = get_getrf_addr(Val(T), Val(BlasInt))
    info = Ref{BlasInt}()
    ccall(
        fptr, Cvoid, (Ref{BlasInt},Ref{BlasInt},Ptr{T},Ref{BlasInt},Ptr{BlasInt},Ref{BlasInt}),
        M % BlasInt, N % BlasInt, A, max(1,stride(A,2)) % BlasInt, ipiv, info
    )
    LU(A, ipiv, info[])
end
lu!(A::StridedMatrix, ::Val{BlasInt} = default_blas_int()) where {BlasInt} = getrf!(A, getrf_ipiv(A, Val(BlasInt)))
lu(A::StridedMatrix, ::Val{BlasInt} = default_blas_int()) where {BlasInt} = lu!(copy(A), Val(BlasInt))

uplochar(::Union{LowerTriangular,UnitLowerTriangular}) = 'L'
uplochar(::Union{UpperTriangular,UnitUpperTriangular}) = 'U'
diagchar(::Union{LowerTriangular,UpperTriangular}) = 'N'
diagchar(::Union{UnitLowerTriangular,UnitUpperTriangular}) = 'U'

# transpose wrappers should be external to the triangular wrappers of `trsm!`
function trsm!(
    B::StridedMatrix{T}, α::T, A::AbstractMatrix{T}, side::Char, ::Val{BlasInt} = default_blas_int()
) where {T<:Union{Float32,Float64},BlasInt<:Union{Int32,Int64}}
    M, N = size(A)
    pA, transa = untranspose_flag(A)
    uplo = uplochar(pA)
    diag = diagchar(pA)
    fptr = get_trsm_addr(Val(T), Val(BlasInt))
    ppA = parent(pA)
    ccall(
        fptr, Cvoid, (Ref{UInt8}, Ref{UInt8}, Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ref{T}, Ptr{Float64}, Ref{BlasInt}, Ptr{T}, Ref{BlasInt}),
        side, uplo, ifelse(transa, 'T', 'N'), diag, M % BlasInt, N % BlasInt, α, ppA, max(1,stride(ppA,2)) % BlasInt, B, max(1,stride(B,2)) % BlasInt
    )
    B
end

function rdiv!(
    A::StridedMatrix{T}, B::AbstractTriangular{T}, ::Val{BlasInt}
) where {T <: Union{Float32,Float64}, BlasInt <: Union{Int32,Int64}}
    trsm!(A, one(T), B, 'R', Val(BlasInt))
end
function ldiv!(
    A::StridedMatrix{T}, B::AbstractTriangular{T}, ::Val{BlasInt}
) where {T <: Union{Float32,Float64}, BlasInt <: Union{Int32,Int64}}
    trsm!(A, one(T), B, 'L', Val(BlasInt))
end

function _ipiv_rows!(A::LU, order::OrdinalRange, B::StridedVecOrMat)
    @inbounds for i ∈ order
        i ≠ A.ipiv[i] && LinearAlgebra._swap_rows!(B, i, A.ipiv[i])
    end
    B
end
_apply_ipiv_rows!(A::LU, B::StridedVecOrMat) = _ipiv_rows!(A, 1 : length(A.ipiv), B)
function _ipiv_cols!(A::LU, order::OrdinalRange, B::StridedVecOrMat)
    ipiv = A.ipiv
    @inbounds for i ∈ order
        i ≠ ipiv[i] && LinearAlgebra._swap_cols!(B, i, ipiv[i])
    end
    B
end
_apply_inverse_ipiv_cols!(A::LU, B::StridedVecOrMat) = _ipiv_cols!(A, length(A.ipiv) : -1 : 1, B)
function rdiv!(A::AbstractMatrix, B::LU{<:Any,<:AbstractMatrix,BlasInt}) where {BlasInt}
    rdiv!(rdiv!(A, UpperTriangular(B.factors), Val(BlasInt)), UnitLowerTriangular(B.factors), Val(BlasInt))
    _apply_inverse_ipiv_cols!(B, A) # mutates `A`
end
function ldiv!(A::LU{<:Any,<:AbstractMatrix,BlasInt}, B::AbstractMatrix) where {BlasInt}
    _apply_ipiv_rows!(A, B)
    ldiv!(ldiv!(B, UnitLowerTriangular(A.factors), Val(BlasInt)), UpperTriangular(A.factors), Val(BlasInt))
end

function rdiv!(A::AbstractMatrix, B::AbstractMatrix, ::Val{BlasInt} = default_blas_int()) where {BlasInt}
    rdiv!(A, lu!(B, Val(BlasInt)))
end
function ldiv!(A::AbstractMatrix, B::AbstractMatrix, ::Val{BlasInt} = default_blas_int()) where {BlasInt}
    ldiv!(lu!(A, Val(BlasInt)), B)
end

rdiv(A::AbstractMatrix, B, ::Val{BlasInt} = default_blas_int()) where {BlasInt} = rdiv!(copy(A), copy(B), Val(BlasInt))
ldiv(A::AbstractMatrix, B, ::Val{BlasInt} = default_blas_int()) where {BlasInt} = ldiv!(copy(A), copy(B), Val(BlasInt))

function __init__()
    LIBBLASTRAMPOLINE_HANDLE[] = libblastrampoline_jll.libblastrampoline_handle
    Sys.isapple() && BLAS.lbt_forward("/System/Library/Frameworks/Accelerate.framework/Versions/A/Accelerate")
end

end # module


text: const SEPARATOR = "\n\n" * "#"^80 * "\n\n"


text:
julia> print(j2(41,3,3))
Prisoner killing order: 2, 5, 8, 11, 14, 17, 20, 23, 26, 29, 32, 35, 38, 0, 4, 9, 13, 18, 22, 27, 31, 36, 40, 6, 12, 19, 25, 33, 39, 7, 16, 28, 37, 10, 24, 1, 21, 3.
Survivors: [15, 30, 34]

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
### A Pluto.jl notebook ###
# v0.12.20

using Markdown
using InteractiveUtils

# ╔═╡ e7c05ce8-58ef-11eb-2da6-bfffd4491883
using Pkg, DrWatson

# ╔═╡ 7fb874b8-58f0-11eb-3e17-91468135856f
begin
	#@quickactivate "StatisticalRethinkingStan"
	using StanSample, StanOptimize
	using StatisticalRethinking
end

# ╔═╡ f7e9345a-58ef-11eb-347d-1b1c98f3663b
md" ## Clip-07-12s.jl"

# ╔═╡ 3c15dcca-58f3-11eb-046c-33c66612701d
md" ##### Entropy."

# ╔═╡ f7e8f63e-58ef-11eb-094b-85d2f930cbfc
H(p) = -sum(p .* log.(p))

# ╔═╡ 4a2cd688-58f3-11eb-071a-7d5191b2acc4
md" ##### Cross Entropy."

# ╔═╡ 313dbcfa-58f3-11eb-0b61-73d83a62cb76
H(p, q) = -sum(p .* log.(q))

# ╔═╡ 514bd3c4-58f3-11eb-27eb-055303c639b1
md" ##### Kullback-Leibler divergence."

# ╔═╡ 37a70450-58f0-11eb-00b1-c526f9615ff4
D(p, q) = sum(p .* log.(p ./ q))

# ╔═╡ 5e898a3a-58f4-11eb-1c3d-ef781cf3e9f0
md" ### Snippet 7.12"

# ╔═╡ 48b01fc0-58f0-11eb-10e1-8969cfc5022f
begin
	p = [0.3, 0.7]
	q = [0.25, 0.75]
	earth = [0.7, 0.3]
	mars = [0.01, 0.99]
end

# ╔═╡ 65e714ae-58f0-11eb-0f86-bb6dad500c12
H(p)

# ╔═╡ 90e27340-58f2-11eb-367e-cb2aae23b616
H([0.01, 0.99])

# ╔═╡ c0ee72a4-6668-11eb-1078-a5f09dfbecff
H([0.7, 0.15, 0.15])

# ╔═╡ 69d6127c-58f0-11eb-29c9-0ba1d12c61c2
D(p, q)

# ╔═╡ b6ffed52-58f0-11eb-1ab2-7902bd8a2b52
begin
	qrange = 0.001:0.01:1.0
	res = []
	for qstep in qrange
		qs = [qstep, 1-qstep]
		append!(res, [D(p, qs)])
	end
	plot(qrange, res, xlab="q[1]", ylab="Divergence q from p", leg=false)
	vline!([0.3])
end

# ╔═╡ edf40570-58f3-11eb-0650-03caff8b1c72
H(p, q)

# ╔═╡ 42415416-58f4-11eb-2cc3-c951f0416798
D(p, q)

# ╔═╡ 49cab272-58f4-11eb-02a8-375cce076851
H(p, q) - H(p)

# ╔═╡ 2c131646-58f6-11eb-12f2-f3b8dfb82dcb
md" ##### Divergence from earth -> mars."

# ╔═╡ 09550666-5994-11eb-3c4c-11461379f117
md"
!!! note
    Reverse arguments?
"

# ╔═╡ 755016a2-58f5-11eb-02c6-235a5d44e9c7
D(mars, earth)

# ╔═╡ 40d5189a-58f6-11eb-08c3-95d784d16c87
md" ##### Divergence from mars -> earth."

# ╔═╡ 6bc4b400-58f5-11eb-1de2-c9c8b44ae7ae
D(earth, mars)

# ╔═╡ 0ffd4f9a-58f0-11eb-1f66-236f6b48ab3f
md" ## End of clip-07-12s.jl"

# ╔═╡ Cell order:
# ╟─f7e9345a-58ef-11eb-347d-1b1c98f3663b
# ╠═e7c05ce8-58ef-11eb-2da6-bfffd4491883
# ╠═7fb874b8-58f0-11eb-3e17-91468135856f
# ╟─3c15dcca-58f3-11eb-046c-33c66612701d
# ╠═f7e8f63e-58ef-11eb-094b-85d2f930cbfc
# ╟─4a2cd688-58f3-11eb-071a-7d5191b2acc4
# ╠═313dbcfa-58f3-11eb-0b61-73d83a62cb76
# ╟─514bd3c4-58f3-11eb-27eb-055303c639b1
# ╠═37a70450-58f0-11eb-00b1-c526f9615ff4
# ╟─5e898a3a-58f4-11eb-1c3d-ef781cf3e9f0
# ╠═48b01fc0-58f0-11eb-10e1-8969cfc5022f
# ╠═65e714ae-58f0-11eb-0f86-bb6dad500c12
# ╠═90e27340-58f2-11eb-367e-cb2aae23b616
# ╠═c0ee72a4-6668-11eb-1078-a5f09dfbecff
# ╠═69d6127c-58f0-11eb-29c9-0ba1d12c61c2
# ╠═b6ffed52-58f0-11eb-1ab2-7902bd8a2b52
# ╠═edf40570-58f3-11eb-0650-03caff8b1c72
# ╠═42415416-58f4-11eb-2cc3-c951f0416798
# ╠═49cab272-58f4-11eb-02a8-375cce076851
# ╟─2c131646-58f6-11eb-12f2-f3b8dfb82dcb
# ╟─09550666-5994-11eb-3c4c-11461379f117
# ╠═755016a2-58f5-11eb-02c6-235a5d44e9c7
# ╟─40d5189a-58f6-11eb-08c3-95d784d16c87
# ╠═6bc4b400-58f5-11eb-1de2-c9c8b44ae7ae
# ╟─0ffd4f9a-58f0-11eb-1f66-236f6b48ab3f


text:
import Pkg

MIRROR = "https://zenodo.org/record/3633228/files/"
TARFILES = ["ModPolCls.tar", "ModPolEta.tar", "ModPolAtk.tar", "ModPolWeb.tar"]

cd(Pkg.dir("ClassPolynomials", "deps"))

for f in TARFILES
    download("$MIRROR/$f", f)
    run(`tar xf $f`)
    rm(f)
end


text:
import DecisionTree

function DecisionTreeModel(
        single_label_name::Symbol;
        name::AbstractString = "",
        nsubfeatures::Integer = 2,
        ntrees::Integer = 20,
        isclassificationmodel::Bool = false,
        isregressionmodel::Bool = false,
        levels::AbstractVector = [],
        )
    hyperparameters = Dict()
    hyperparameters[:nsubfeatures] = nsubfeatures
    hyperparameters[:ntrees] = ntrees
    hyperparameters = fix_type(hyperparameters)
    underlyingrandomforest = FitNotYetRunUnderlyingObject()
    result = DecisionTreeModel(
        name,
        isclassificationmodel,
        isregressionmodel,
        single_label_name,
        levels,
        hyperparameters,
        underlyingrandomforest,
        )
    return result
end

"""
"""
function get_underlying(
        x::DecisionTreeModel;
        saving::Bool = false,
        loading::Bool = false,
        )
    result = x.underlyingrandomforest
    return result
end

"""
"""
function fit!(
        estimator::DecisionTreeModel,
        featuresarray::AbstractArray,
        labelsarray::AbstractArray,
        )
    @info(string("Starting to train DecisionTree model."))
    randomforest = try
        DecisionTree.build_forest(
            labelsarray,
            featuresarray,
            estimator.hyperparameters[:nsubfeatures],
            estimator.hyperparameters[:ntrees],
            )
    catch e
        @warn(
            string(
                "While training DecisionTree model, ignored error: ",
                e,
                )
            )
        FitFailedUnderlyingObject()
    end
    @info(string("Finished training DecisionTree model."))
    estimator.underlyingrandomforest = randomforest
    return estimator
end

"""
"""
function predict(
        estimator::DecisionTreeModel,
        featuresarray::AbstractArray,
        )
    if estimator.isclassificationmodel && !estimator.isregressionmodel
        probabilitiesassoc = predict_proba(
            estimator,
            featuresarray,
            )
        predictionsvector = single_labelprobabilitiestopredictions(
            probabilitiesassoc
            )
        return predictionsvector
    elseif !estimator.isclassificationmodel && estimator.isregressionmodel
        if is_nothing(estimator.underlyingrandomforest)
            predicted_values = fill(Float64(0), size(featuresarray,1))
        else
            predicted_values = DecisionTree.apply_forest(
                estimator.underlyingrandomforest,
                featuresarray,
                )
        end
        return predicted_values
    else
        error(
            "Could not figure out if model is classification or regression"
            )
    end
end

"""
"""
function predict(
        estimator::DecisionTreeModel,
        featuresarray::AbstractArray,
        positive_class::AbstractString,
        threshold::AbstractFloat,
        )
    if estimator.isclassificationmodel && !estimator.isregressionmodel
        probabilitiesassoc = predict_proba(
            estimator,
            featuresarray,
            )
        predictionsvector = single_labelprobabilitiestopredictions(
            probabilitiesassoc,
            positive_class,
            threshold,
            )
        return predictionsvector
    else
        error(
            "Can only use the `threshold` argument with classification models"
            )
    end
end

"""
"""
function predict_proba(
        estimator::DecisionTreeModel,
        featuresarray::AbstractArray,
        )
    if estimator.isclassificationmodel && !estimator.isregressionmodel
        if is_nothing(estimator.underlyingrandomforest)
            predictedprobabilities = fill(
                Float64(0),
                size(featuresarray, 1),
                length(estimator.levels),
                )
            predictedprobabilities[:, 1] = 1
        else
            predictedprobabilities = DecisionTree.apply_forest_proba(
                estimator.underlyingrandomforest,
                featuresarray,
                estimator.levels,
                )
        end
        result = Dict()
        for i = 1:length(estimator.levels)
            result[estimator.levels[i]] = predictedprobabilities[:, i]
        end
        result = fix_type(result)
        return result
    elseif !estimator.isclassificationmodel && estimator.isregressionmodel
        error("predict_proba is not defined for regression models")
    else
        error(
            "Could not figure out if model is classification or regression"
            )
    end
end

"""
"""
function single_labelmulticlassdfrandomforestclassifier_DecisionTree(
        feature_names::AbstractVector,
        single_label_name::Symbol,
        single_label_levels::AbstractVector;
        name::AbstractString = "",
        nsubfeatures::Integer = 2,
        ntrees::Integer = 10,
        feature_contrasts::Union{Nothing, AbstractFeatureContrasts} = nothing,
        )
    dftransformer = MutableDataFrame2DecisionTreeTransformer(
        feature_names,
        single_label_name;
        levels = single_label_levels,
        )
    randomforestestimator = DecisionTreeModel(
        single_label_name;
        name = name,
        nsubfeatures = nsubfeatures,
        ntrees = ntrees,
        isclassificationmodel = true,
        isregressionmodel = false,
        levels = single_label_levels,
        )
    probapackager = ImmutablePackageSingleLabelPredictProbaTransformer(
        single_label_name,
        )
    predpackager = ImmutablePackageSingleLabelPredictionTransformer(
        single_label_name,
        )
    finalpipeline = dftransformer |>
                    randomforestestimator |>
                    probapackager |>
                    predpackager
    finalpipeline.name = name
    if !is_nothing(feature_contrasts)
        set_feature_contrasts!(finalpipeline, feature_contrasts)
    end
    return finalpipeline
end

"""
"""
function single_labelmulticlassdataframerandomforestclassifier(
        feature_names::AbstractVector,
        single_label_name::Symbol,
        single_label_levels::AbstractVector;
        name::AbstractString = "",
        package::Symbol = :none,
        nsubfeatures::Integer = 2,
        ntrees::Integer = 10,
        feature_contrasts::Union{Nothing, AbstractFeatureContrasts} = nothing,
        )
    if package == :DecisionTree
        result =
            single_labelmulticlassdfrandomforestclassifier_DecisionTree(
                feature_names,
                single_label_name,
                single_label_levels;
                name = name,
                nsubfeatures = nsubfeatures,
                ntrees = ntrees,
                feature_contrasts = feature_contrasts
                )
        return result
    else
        error("$(package) is not a valid value for package")
    end
end

"""
"""
function single_labeldataframerandomforestregression_DecisionTree(
        feature_names::AbstractVector,
        single_label_name::Symbol;
        name::AbstractString = "",
        nsubfeatures::Integer = 2,
        ntrees::Integer = 10,
        feature_contrasts::Union{Nothing, AbstractFeatureContrasts} = nothing,
        )
    dftransformer = MutableDataFrame2DecisionTreeTransformer(
        feature_names,
        single_label_name,
        )
    randomforestestimator = DecisionTreeModel(
        single_label_name;
        name = name,
        nsubfeatures = nsubfeatures,
        ntrees = ntrees,
        isclassificationmodel = false,
        isregressionmodel = true,
        )
    predpackager = ImmutablePackageSingleLabelPredictionTransformer(
        single_label_name,
        )
    finalpipeline = dftransformer |>
                    randomforestestimator |>
                    predpackager
    finalpipeline.name = name
    if !is_nothing(feature_contrasts)
        set_feature_contrasts!(finalpipeline, feature_contrasts)
    end
    return finalpipeline
end

"""
"""
function single_labeldataframerandomforestregression(
        feature_names::AbstractVector,
        single_label_name::Symbol;
        name::AbstractString = "",
        package::Symbol = :none,
        nsubfeatures::Integer = 2,
        ntrees::Integer = 10,
        feature_contrasts::Union{Nothing, AbstractFeatureContrasts} = nothing,
        )
    if package == :DecisionTree
        result = single_labeldataframerandomforestregression_DecisionTree(
            feature_names,
            single_label_name;
            name = name,
            nsubfeatures = nsubfeatures,
            ntrees = ntrees,
            feature_contrasts = feature_contrasts
        )
        return result
    else
        error("$(package) is not a valid value for package")
    end
end


text:
import BitBasis: BitStr, BitStr64

abstract type AbstractArrayReg{D,T,AT} <: AbstractRegister{D} end

struct NoBatch end
_asint(x::Int) = x
_asint(::NoBatch) = 1

"""
    ArrayReg{D,T,MT<:AbstractMatrix{T}} <: AbstractArrayRegister{D}
    ArrayReg{D}(raw)
    ArrayReg(raw::AbstractVecOrMat; nlevel=2)
    ArrayReg(r::ArrayReg)

Simulated full amplitude register type, it uses an array to represent
corresponding one or a batch of quantum states. `T`
is the numerical type for each amplitude, it is `ComplexF64` by default.

!!! warning

    `ArrayReg` constructor will not normalize the quantum state. If you need a
    normalized quantum state remember to use `normalize!(register)` on the register or
    normalize the input raw array with `normalize` or [`batched_normalize!`](@ref).
"""
mutable struct ArrayReg{D,T,MT<:AbstractMatrix{T}} <: AbstractArrayReg{D,T,MT}
    state::MT
    function ArrayReg{D,T,MT}(state::MT) where {D,T, MT<:AbstractMatrix{T}}
        _check_reg_input(state, D, 1)
        return new{D,T,MT}(state)
    end
end

ArrayReg(raw; nlevel=2) = ArrayReg{nlevel}(raw)
ArrayReg{D}(raw::AbstractVector) where D = ArrayReg{D}(reshape(raw, :, 1))
function ArrayReg{D}(raw::MT) where {D,T,MT<:AbstractMatrix{T}}
    return ArrayReg{D,T,MT}(raw)
end

ArrayReg(r::ArrayReg{D}) where {D} = ArrayReg{D}(copy(r.state))
ArrayReg(r::ArrayReg{D,T,<:Transpose}) where {D,T} =
        ArrayReg{D}(Transpose(copy(r.state.parent)))

Base.copy(r::ArrayReg) = ArrayReg(r)
Base.similar(r::ArrayReg{D}) where D = ArrayReg{D}(similar(state(r)))
Base.similar(r::ArrayReg{D}, state::AbstractMatrix) where D = ArrayReg{D}(state)

"""
    BatchedArrayReg{D,T,MT<:AbstractMatrix{T}} <: AbstractArrayReg{D}
    BatchedArrayReg(raw, nbatch; nlevel=2)
    BatchedArrayReg{D}(raw, nbatch)

Simulated batched full amplitude register type, it uses an array to represent
corresponding one or a batch of quantum states. `T`
is the numerical type for each amplitude, it is `ComplexF64` by default.

!!! warning

    `BatchedArrayReg` constructor will not normalize the quantum state. If you need a
    normalized quantum state remember to use `normalize!(register)` on the register or
    normalize the input raw array with `normalize` or [`batched_normalize!`](@ref).
"""
mutable struct BatchedArrayReg{D,T,MT<:AbstractMatrix{T}} <: AbstractArrayReg{D,T,MT}
    state::MT
    nbatch::Int
    function BatchedArrayReg{D,T,MT}(state::MT, nbatch::Int) where {D,T, MT<:AbstractMatrix{T}}
        _check_reg_input(state, D, nbatch)
        return new{D,T,MT}(state, nbatch)
    end
end
BatchedArrayReg(raw::AbstractMatrix, nbatch::Int=size(raw, 2); nlevel=2) = BatchedArrayReg{nlevel}(raw, nbatch)
function BatchedArrayReg{D}(raw::MT, nbatch::Int) where {D,T,MT<:AbstractMatrix{T}}
    return BatchedArrayReg{D,T,MT}(raw, nbatch)
end

BatchedArrayReg(r::BatchedArrayReg{D}) where {D} = BatchedArrayReg{D}(copy(r.state), r.nbatch)
BatchedArrayReg(r::BatchedArrayReg{D,T,<:Transpose}) where {D,T} =
        BatchedArrayReg{D}(Transpose(copy(r.state.parent)), r.nbatch)

Base.copy(r::BatchedArrayReg) = BatchedArrayReg(r)
Base.similar(r::BatchedArrayReg{D}) where D = BatchedArrayReg{D}(similar(state(r)), r.nbatch)
Base.similar(r::BatchedArrayReg{D}, state::AbstractMatrix) where D = BatchedArrayReg{D}(state, r.nbatch)
YaoBase.viewbatch(r::ArrayReg, ind::Int) = ind == 1 ? r : error("Index `$ind` out of bounds, should be `1`.")

# convert
function ArrayReg(r::BatchedArrayReg{D}) where D
    if nbatch(r) != 1
        error("can not convert a `BatchedArrayReg` with `nbatch != 1` to a `ArrayReg`")
    else
        return ArrayReg{D}(r.state)
    end
end
BatchedArrayReg(r::ArrayReg{D}...) where D = BatchedArrayReg{D}(hcat(state.(r)...), length(r))

function _check_reg_input(raw::AbstractMatrix{T}, D::Integer, B::Integer) where T
    _warn_type(raw)
    if D <= 0
        error("invalid number of level: $D")
    end
    if B < 0
        error("invalid batch size: $B")
    end
    ispow(size(raw, 1), D) ||
        throw(DimensionMismatch("Expect first dimension size to be power of $D"))
    if !(ispow(size(raw, 2) ÷ B, D) && size(raw, 2) % B == 0)
        throw(
            DimensionMismatch(
                "Expect second dimension size to be an integer multiple of batch size $B",
            ),
        )
    end
end

function _warn_type(raw::AbstractMatrix{T}) where T
    T <: Complex || @warn "Input matrix element type is not `Complex`, got `$(eltype(raw))`"
end

"""
    arrayreg(state; nbatch::Union{Integer,NoBatch}=NoBatch(), nlevel::Integer=2)

Create an array register, if nbatch is a integer, it will return a `BatchedArrayReg`.
"""
function arrayreg(state; nbatch::Union{Integer,NoBatch}=NoBatch(), nlevel::Integer=2)
    if nbatch isa NoBatch
        return ArrayReg{nlevel}(state)
    else
        return BatchedArrayReg{nlevel}(state, nbatch)
    end
end

Adapt.@adapt_structure ArrayReg
Adapt.@adapt_structure BatchedArrayReg

const AdjointArrayReg{D,T,MT} = AdjointRegister{D,<:AbstractArrayReg{D,T,MT}}
const ArrayRegOrAdjointArrayReg{D,T,MT} =
    Union{AbstractArrayReg{D,T,MT},AdjointArrayReg{D,T,MT}}

"""
    datatype(register) -> Int

Returns the numerical data type used by register.

!!! note

    `datatype` is not the same with `eltype`, since `AbstractRegister` family
    is not exactly the same with `AbstractArray`, it is an iterator of several
    registers.
"""
datatype(r::AbstractArrayReg{D,T}) where {D,T} = T

"""
    arrayreg([T=ComplexF64], bit_str; nbatch=NoBatch())

Construct an array register from bit string literal.
For bit string literal please read [`@bit_str`](@ref).

# Examples

```jldoctest; setup=:(using YaoArrayRegister)
julia> arrayreg(bit"1010")
ArrayReg{2, ComplexF64, Array...}
    active qudits: 4/4
    nlevel: 2

julia> arrayreg(ComplexF32, bit"1010")
ArrayReg{2, ComplexF32, Array...}
    active qudits: 4/4
    nlevel: 2
```
"""
arrayreg(bitstr::BitStr; nbatch::Union{Int,NoBatch}=NoBatch()) = arrayreg(ComplexF64, bitstr; nbatch=nbatch)
arrayreg(::Type{T}, bitstr::BitStr; nbatch::Union{Int,NoBatch}=NoBatch()) where {T} = arrayreg(onehot(T, bitstr, _asint(nbatch)); nbatch=nbatch, nlevel=2)

transpose_storage(reg::AbstractArrayReg{D,T,<:Transpose}) where {D,T} = arrayreg(copy(reg.state); nbatch=nbatch(reg), nlevel=D)
transpose_storage(reg::AbstractArrayReg) =
    arrayreg(transpose(copy(transpose(reg.state))); nbatch=nbatch(reg), nlevel=nlevel(reg))

# NOTE: ket bra is not copyable
function Base.copyto!(dst::AbstractArrayReg, src::AbstractArrayReg)
    copyto!(state(dst), state(src))
    return dst
end

function Base.copyto!(dst::AdjointArrayReg, src::AdjointArrayReg)
    copyto!(state(dst), state(src))
    return dst
end

# register interface
YaoBase.nqudits(r::AbstractArrayReg{2}) = log2i(length(r.state) ÷ _asint(nbatch(r)))
YaoBase.nqudits(r::AbstractArrayReg{D}) where {D} = logdi(length(r.state) ÷ _asint(nbatch(r)), D)
YaoBase.nactive(r::AbstractArrayReg{D}) where {D} = logdi(size(r.state, 1), D)
YaoBase.viewbatch(r::BatchedArrayReg{D}, ind::Int) where D = @inbounds ArrayReg{D}(view(rank3(r), :, :, ind))

function YaoBase.addbits!(r::AbstractArrayReg{D}, n::Int) where {D}
    raw = state(r)
    M, N = size(raw)
    r.state = similar(r.state, M * (D ^ n), N)
    fill!(r.state, 0)
    r.state[1:M, :] = raw
    return r
end

function YaoBase.insert_qudits!(reg::AbstractArrayReg{D}, loc::Int; nqudits::Int = 1) where D
    na = nactive(reg)
    focus!(reg, 1:loc-1)
    reg2 = join(zero_state(nqudits; nbatch = nbatch(reg), nlevel=D), reg) |> relax! |> focus!((1:na+nqudits)...)
    reg.state = reg2.state
    reg
end

function YaoBase.probs(r::ArrayReg)
    if size(r.state, 2) == 1
        return vec(r.state .|> abs2)
    else
        return dropdims(sum(r.state .|> abs2, dims = 2), dims = 2)
    end
end

function YaoBase.probs(r::BatchedArrayReg)
    if size(r.state, 2) == nbatch(r)
        return r.state .|> abs2
    else
        probs = r |> rank3 .|> abs2
        return dropdims(sum(probs, dims = 2), dims = 2)
    end
end

function YaoBase.reorder!(r::AbstractArrayReg, orders)
    @inbounds for i = 1:size(r.state, 2)
        r.state[:, i] = reorder(r.state[:, i], orders)
    end
    return r
end

function YaoBase.collapseto!(r::AbstractArrayReg, bit_config::Integer)
    st = normalize!(r.state[Int(bit_config)+1, :])
    fill!(r.state, 0)
    r.state[Int(bit_config)+1, :] .= st
    return r
end

"""
    fidelity(r1::AbstractArrayReg, r2::AbstractArrayReg)

Calcuate the fidelity between `r1` and `r2`, if `r1` or `r2` is not pure state
(`nactive(r) != nqudits(r)`), the fidelity is calcuated by purification. See also
[`pure_state_fidelity`](@ref), [`purification_fidelity`](@ref).

    fidelity'(pair_or_reg1, pair_or_reg2) -> (g1, g2)

Obtain the gradient with respect to registers and circuit parameters.
For pair input `ψ=>circuit`, the returned gradient is a pair of `gψ=>gparams`,
with `gψ` the gradient of input state and `gparams` the gradients of circuit parameters.
For register input, the return value is a register.
"""
function YaoBase.fidelity(r1::BatchedArrayReg{D}, r2::BatchedArrayReg{D}) where {D}
    B1, B2 = nbatch(r1), nbatch(r2)
    B1 == B2 || throw(DimensionMismatch("Register batch not match!"))
    B = nbatch(r1)

    state1 = rank3(r1)
    state2 = rank3(r2)
    size(state1) == size(state2) || throw(DimensionMismatch("Register size not match!"))
    if size(state1, 2) == 1
        res = map(b -> pure_state_fidelity(state1[:, 1, b], state2[:, 1, b]), 1:B)
    else
        res = map(b -> purification_fidelity(state1[:, :, b], state2[:, :, b]), 1:B)
    end
    return res
end

function YaoBase.fidelity(r1::BatchedArrayReg, r2::ArrayReg)
    B = nbatch(r1)
    state1 = rank3(r1)
    state2 = rank3(r2)
    nqudits(r1) == nqudits(r2) || throw(DimensionMismatch("Register size not match!"))
    if size(state1, 2) == 1
        res = map(b -> pure_state_fidelity(state1[:, 1, b], state2[:, 1, 1]), 1:B)
    else
        res = map(b -> purification_fidelity(state1[:, :, b], state2[:, :, 1]), 1:B)
    end
    return res
end

YaoBase.fidelity(r1::ArrayReg, r2::BatchedArrayReg) = YaoBase.fidelity(r2, r1)

function YaoBase.fidelity(r1::ArrayReg, r2::ArrayReg)
    state1 = state(r1)
    state2 = state(r2)
    nqudits(r1) == nqudits(r2) || throw(DimensionMismatch("Register size not match!"))

    if size(state1, 2) == 1
        return pure_state_fidelity(state1[:, 1], state2[:, 1])
    else
        return purification_fidelity(state1, state2)
    end
end

YaoBase.tracedist(r1::ArrayReg, r2::ArrayReg) = tracedist(ρ(r1), ρ(r2))
YaoBase.tracedist(r1::BatchedArrayReg, r2::BatchedArrayReg) = tracedist.(r1, r2)

# properties
"""
    state(register::AbstractArrayReg) -> raw array

Returns the raw array storage of `register`. See also [`statevec`](@ref).
"""
state(r::AbstractArrayReg) = r.state
state(r::AdjointArrayReg) = adjoint(state(parent(r)))

"""
    statevec(r::ArrayReg) -> array

Return a state matrix/vector by droping the last dimension of size 1. See also [`state`](@ref).

!!! warning
    `statevec` is not type stable. It may cause performance slow down.
"""
statevec(r::ArrayRegOrAdjointArrayReg) = matvec(state(r))

"""
    relaxedvec(r::AbstractArrayReg) -> AbstractArray

Return a vector representation of state, with all qudits activated.
See also [`state`](@ref), [`statevec`](@ref).
"""
relaxedvec(r::ArrayReg) = vec(state(r))
relaxedvec(r::BatchedArrayReg) = reshape(state(r), :, nbatch(r))

"""
    hypercubic(r::ArrayReg) -> AbstractArray

Return the hypercubic form (high dimensional tensor) of this register, only active qudits are considered.
See also [`rank3`](@ref).
"""
BitBasis.hypercubic(r::ArrayRegOrAdjointArrayReg) =
    reshape(state(r), ntuple(i -> 2, Val(nactive(r)))..., :)

"""
    rank3(r::ArrayReg)

Return the rank 3 tensor representation of state,
the 3 dimensions are (activated space, remaining space, batch dimension).
See also [`rank3`](@ref).
"""
function rank3(r::ArrayRegOrAdjointArrayReg)
    s = state(r)
    reshape(s, size(s, 1), :, _asint(nbatch(r)))
end

"""
    join(regs...)

concat a list of registers `regs` to a larger register, each register should
have the same batch size. See also [`repeat`](@ref).
"""
Base.join(rs::AbstractArrayReg...) = _join(join_datatype(rs...), rs...)
Base.join(r::AbstractArrayReg) = r
function _join(::Type{T}, rs0::AbstractArrayReg{D}, rs::AbstractArrayReg{D}...) where {T,D}
    state = batched_kron(rank3(rs0), rank3.(rs)...)
    return arrayreg(reshape(state, size(state, 1), :), nbatch=nbatch(rs[1]), nlevel=D)
end

join_datatype(r::AbstractArrayReg{D,T}, rs::AbstractArrayReg{D,T}...) where {D,T} = join_datatype(T, r, rs...)
join_datatype(::Type{T}, r::AbstractArrayReg{D,T1}, rs::AbstractArrayReg...) where {T,T1,D} =
    join_datatype(promote_type(T, T1), rs...)
join_datatype(::Type{T}) where {T} = T

# initialization methods
"""
    product_state([T=ComplexF64], bit_str; nbatch=NoBatch())

Create an [`ArrayReg`](@ref) with bit string literal
defined with [`@bit_str`](@ref). See also [`zero_state`](@ref),
[`rand_state`](@ref), [`uniform_state`](@ref).

# Examples

```jldoctest; setup=:(using YaoArrayRegister)
julia> product_state(bit"100"; nbatch=2)
BatchedArrayReg{2, ComplexF64, Transpose...}
    active qudits: 3/3
    nlevel: 2
    nbatch: 2

julia> r1 = product_state(ComplexF32, bit"100"; nbatch=2)
BatchedArrayReg{2, ComplexF32, Transpose...}
    active qudits: 3/3
    nlevel: 2
    nbatch: 2

julia> r2 = product_state(ComplexF32, [0, 0, 1]; nbatch=2)
BatchedArrayReg{2, ComplexF32, Transpose...}
    active qudits: 3/3
    nlevel: 2
    nbatch: 2

julia> r1 ≈ r2   # because we read bit strings from right to left, vectors from left to right.
true
```
"""
product_state(bit_str::BitStr; nbatch::Union{NoBatch,Int} = NoBatch()) =
    product_state(ComplexF64, bit_str; nbatch = nbatch)

product_state(bit_str::AbstractVector; nbatch::Union{NoBatch,Int} = NoBatch()) =
    product_state(ComplexF64, bit_str; nbatch = nbatch)

"""
    product_state([T=ComplexF64], total::Int, bit_config::Integer; nbatch=1, no_transpose_storage=false)

Create an [`ArrayReg`](@ref) with bit configuration `bit_config`, total number of bits `total`.
See also [`zero_state`](@ref), [`rand_state`](@ref), [`uniform_state`](@ref).

# Examples

```jldoctest; setup=:(using YaoArrayRegister)
julia> product_state(4, 3; nbatch=2)
BatchedArrayReg{2, ComplexF64, Transpose...}
    active qudits: 4/4
    nlevel: 2
    nbatch: 2

julia> product_state(4, 0b1001; nbatch=2)
BatchedArrayReg{2, ComplexF64, Transpose...}
    active qudits: 4/4
    nlevel: 2
    nbatch: 2

julia> product_state(ComplexF32, 4, 0b101)
ArrayReg{2, ComplexF32, Array...}
    active qudits: 4/4
    nlevel: 2
```

!!! warning

    This interface will not check whether the number of required digits
    for the bit configuration matches the total number of bits.
"""
product_state(total::Int, bit_config::Integer; kwargs...) =
    product_state(ComplexF64, total, bit_config; kwargs...)

product_state(::Type{T}, bit_str::BitStr{N}; kwargs...) where {T,N} =
    product_state(T, N, buffer(bit_str); kwargs...)

product_state(::Type{T}, bit_configs::AbstractVector; kwargs...) where {T} =
    product_state(T, bit_literal(bit_configs...); kwargs...)

function product_state(
    ::Type{T},
    total::Int,
    bit_config::Integer;
    nbatch::Union{Int,NoBatch} = NoBatch(),
    nlevel::Int=2,
    no_transpose_storage::Bool = false,
) where {T}
    if nbatch isa NoBatch || no_transpose_storage
        raw = onehot(T, total, bit_config, _asint(nbatch), nlevel)
    else
        raw = zeros(T, _asint(nbatch), nlevel ^ total)
        raw[:, Int(bit_config)+1] .= 1
        raw = transpose(raw)
    end
    return arrayreg(raw; nbatch=nbatch, nlevel=nlevel)
end

function BitBasis.onehot(::Type{T}, nbits::Int, x::Integer, nbatch::Int, nlevel::Int) where {T}
    v = zeros(T, nlevel ^ nbits, nbatch)
    v[x+1, :] .= 1
    return v
end

"""
    zero_state([T=ComplexF64], n::Int; nbatch::Int=1)

Create an [`AbstractArrayReg`](@ref) with total number of bits `n`.
See also [`product_state`](@ref), [`rand_state`](@ref), [`uniform_state`](@ref).

# Examples

```jldoctest; setup=:(using YaoArrayRegister)
julia> zero_state(4)
ArrayReg{2, ComplexF64, Array...}
    active qudits: 4/4
    nlevel: 2

julia> zero_state(ComplexF32, 4)
ArrayReg{2, ComplexF32, Array...}
    active qudits: 4/4
    nlevel: 2

julia> zero_state(ComplexF32, 4; nbatch=3)
BatchedArrayReg{2, ComplexF32, Transpose...}
    active qudits: 4/4
    nlevel: 2
    nbatch: 3
```
"""
zero_state(n::Int; kwargs...) = zero_state(ComplexF64, n; kwargs...)
zero_state(::Type{T}, n::Int; kwargs...) where {T} = product_state(T, n, 0; kwargs...)

"""
    rand_state([T=ComplexF64], n::Int; nbatch=1, no_transpose_storage=false)

Create a random [`AbstractArrayReg`](@ref) with total number of qudits `n`.

# Examples

```jldoctest; setup=:(using YaoArrayRegister)
julia> rand_state(4)
ArrayReg{2, ComplexF64, Array...}
    active qudits: 4/4
    nlevel: 2

julia> rand_state(ComplexF64, 4)
ArrayReg{2, ComplexF64, Array...}
    active qudits: 4/4
    nlevel: 2

julia> rand_state(ComplexF64, 4; nbatch=2)
BatchedArrayReg{2, ComplexF64, Transpose...}
    active qudits: 4/4
    nlevel: 2
    nbatch: 2
```
"""
rand_state(n::Int; kwargs...) = rand_state(ComplexF64, n; kwargs...)

function rand_state(
    ::Type{T},
    n::Int;
    nbatch::Union{Int,NoBatch} = NoBatch(),
    no_transpose_storage::Bool = false,
    nlevel = 2,
) where {T}
    raw =
        nbatch isa NoBatch || no_transpose_storage ? randn(T, nlevel ^ n, _asint(nbatch)) :
        transpose(randn(T, _asint(nbatch), nlevel ^ n))
    return normalize!(arrayreg(raw; nbatch=nbatch, nlevel=nlevel))
end

"""
    uniform_state([T=ComplexF64], n; nbatch=1, no_transpose_storage=false)

Create a uniform state: ``\\frac{1}{2^n} \\sum_k |k⟩``. This state
can also be created by applying [`H`](@ref) (Hadmard gate) on ``|00⋯00⟩`` state.

# Example

```jldoctest; setup=:(using YaoArrayRegister)
julia> uniform_state(4; nbatch=2)
BatchedArrayReg{2, ComplexF64, Transpose...}
    active qudits: 4/4
    nlevel: 2
    nbatch: 2

julia> uniform_state(ComplexF32, 4; nbatch=2)
BatchedArrayReg{2, ComplexF32, Transpose...}
    active qudits: 4/4
    nlevel: 2
    nbatch: 2
```
"""
uniform_state(n::Int; kwargs...) = uniform_state(ComplexF64, n; kwargs...)

function uniform_state(::Type{T}, n::Int;
    nbatch::Union{Int,NoBatch} = NoBatch(),
    nlevel::Int = 2,
    no_transpose_storage::Bool = false,
) where {T}
    if (nbatch isa NoBatch) || no_transpose_storage
        raw = ones(T, nlevel ^ n, _asint(nbatch))
    else
        raw = transpose(ones(T, _asint(nbatch), nlevel ^ n))
    end
    return normalize!(arrayreg(raw; nbatch=nbatch, nlevel=nlevel))
end

"""
    oneto(n::Int) -> f(register)
    oneto(r::AbstractArrayReg, n::Int=nqudits(r))

Returns an register with `1:n` qudits activated.
"""
oneto(n::Int) = r -> oneto(r, n)
oneto(r::AbstractArrayReg{D}, n::Int = nqudits(r)) where {D} =
    arrayreg(reshape(copy(r.state), D ^ n, :), nbatch=nbatch(r), nlevel=D)
oneto(r::AbstractArrayReg{D,T,<:Transpose}, n::Int = nqudits(r)) where {D,T} =
    transpose_storage(arrayreg(reshape(r.state, D ^ n, :); nbatch=nbatch(r), nlevel=D))

"""
    repeat(register, n)

Create an [`ArrayReg`](@ref) by copying the original `register` for `n` times on
batch dimension.

# Example

```jldoctest; setup=:(using YaoArrayRegister)
julia> repeat(arrayreg(bit"101"; nbatch=3), 4)
BatchedArrayReg{2, ComplexF64, Array...}
    active qudits: 3/3
    nlevel: 2
    nbatch: 12
```
"""
Base.repeat(r::AbstractArrayReg{D}, n::Int) where D =
    BatchedArrayReg{D}(hcat((state(r) for k = 1:n)...), n * _asint(nbatch(r)))

# NOTE: overload this to make printing more compact
#       but do not alter the way how type parameters print
function Base.show(io::IO, reg::ArrayReg{D,T,MT}) where {D,T,MT}
    print(io, "ArrayReg{$D, $T, $(nameof(MT))...}")
    print(io, "\n    active qudits: ", nactive(reg), "/", nqudits(reg))
    print(io, "\n    nlevel: ", nlevel(reg))
end

function Base.show(io::IO, reg::BatchedArrayReg{D,T,MT}) where {D,T,MT}
    print(io, "BatchedArrayReg{$D, $T, $(nameof(MT))...}")
    print(io, "\n    active qudits: ", nactive(reg), "/", nqudits(reg))
    print(io, "\n    nlevel: ", nlevel(reg))
    print(io, "\n    nbatch: ", nbatch(reg))
end

"""
    mutual_information(reg::AbstractArrayReg, part1, part2)

Compute the mutual information between locations `part1` and locations `part2` in a quantum state `reg`.
"""
function mutual_information(reg::AbstractArrayReg, part1, part2)
    von_neumann_entropy(reg, part1) .+ von_neumann_entropy(reg, part2) .- von_neumann_entropy(reg, part1 ∪ part2)
end

"""
    von_neumann_entropy(reg::AbstractArrayReg, part)
    von_neumann_entropy(ρ::DensityMatrix)

The entanglement entropy between `part` and the rest part in quantum state `reg`.
If the input is a density matrix, it returns the entropy of a mixed state.
"""
von_neumann_entropy(reg::BatchedArrayReg, part) = von_neumann_entropy.(reg, Ref(part))
von_neumann_entropy(reg::ArrayReg, part) = von_neumann_entropy(density_matrix(reg, part))

function Base.iterate(it::Union{BatchedArrayReg{D}, AdjointRegister{D, <:BatchedArrayReg{D}}} where D, state = 1)
    if state > nbatch(it)
        return nothing
    else
        return viewbatch(it, state), state + 1
    end
end

Base.length(r::BatchedArrayReg) = r.nbatch
Base.length(r::AdjointRegister{D, <:BatchedArrayReg{D}}) where {D} = length(parent(r))

"""
    nbatch(register) -> Int

Returns the number of batches.
"""
nbatch(r::BatchedArrayReg) = r.nbatch
nbatch(r::ArrayReg) = NoBatch()
nbatch(r::AdjointArrayReg) = nbatch(parent(r))


text:
module FWI_64x200

    using LinearAlgebra, PyPlot, JLD, Statistics, Random, Dates
    using JUDI.TimeModeling, SegyIO, JUDI.SLIM_optim, ImageFiltering, IterativeSolvers

    export fwi

    function fwi(m, m0, niterations)

        m = adjoint(m)
        m0 = adjoint(m0)

        n = size(m)
        dm = vec(m - m0)

        # Setup info and model structure
        d = (25., 25.)
        o = (0., 0.)
        nsrc = 21	# number of sources

        ## Set up receiver geometry in meters
        nxrec = 201
        xrec = range(0f0, stop=5000f0, length=nxrec)
        yrec = 0f0
        zrec = range(0f0, stop=0f0, length=nxrec)

        # receiver sampling and recording time
        timeR = 2000f0   # receiver recording time [ms]
        dtR = 1f0    # receiver sampling interval [ms]

        # Set up receiver structure
        recGeometry = Geometry(xrec, yrec, zrec; dt=dtR, t=timeR, nsrc=nsrc)

        ## Set up source geometry (cell array with source locations for each shot)
        xsrc = convertToCell(range(0f0, stop=5000f0, length=nsrc))
        ysrc = convertToCell(range(0f0, stop=0f0, length=nsrc))
        zsrc = convertToCell(range(0f0, stop=0f0, length=nsrc))

        # source sampling and number of time steps
        timeS = timeR  # ms
        dtS = dtR   # ms

        # Set up source structure
        srcGeometry = Geometry(xsrc, ysrc, zsrc; dt=dtS, t=timeS)

        # setup wavelet
        f0 = 0.006f0     # KHz
        wavelet = ricker_wavelet(timeS, dtS, f0)
        q = judiVector(srcGeometry, wavelet);

        # Set up model and info structure for linear operators
        model = Model(n, d, o, m; nb=200)
        model0 = Model(n, d, o, m0; nb=200)
        ntComp = get_computational_nt(srcGeometry, recGeometry, model)
        info = Info(prod(n), nsrc, ntComp)

        opt = Options(space_order=16) # Use linearized inverse scattering imaging condition for the Jacobian

        # Setup operators
        Pr = judiProjection(info, recGeometry)
        F = judiModeling(info, model; options=opt)
        F0 = judiModeling(info, model0; options=opt)
        Ps = judiProjection(info, srcGeometry)
        J = judiJacobian(Pr * F0 * adjoint(Ps), q)

        # Nonlinear modeling
        d_obs = Pr * F * adjoint(Ps) * q # Practical observation data generated on m
        # d_syn = Pr * F0* adjoint(Ps) * q # Synthetic observation data generated on m0
        # qad = Ps * adjoint(F) * adjoint(Pr) * d_obs

        # # Linearized modeling
        # # J.options.file_name = "linearized_shot"
        # dD = J * dm
        # rtm = adjoint(J) * dD # data from linearized modeling (Born modeling)
        # grad = adjoint(J) * (d_syn - d_obs) # data from Nonlinear modeling, which gives the RTM in practice

        # rtm = adjoint(reshape(rtm, n))
        # grad = adjoint(reshape(grad, n))

        # Right-hand preconditioners (model topmute)
        idx_wb = find_water_bottom(reshape(dm, n)) # find the index of the water bottom (you could also manually add this)
        Tm = judiTopmute(model0.n, idx_wb, 2)  # Mute water column
        S = judiDepthScaling(model0) # design a depth scaling operator
        Mr = S * Tm

        # Bound constraints
        vmin = ones(Float32, model0.n) .+ 0.3f0
        vmax = ones(Float32, model0.n) .+ 5.5f0

        mmin = vec((1f0 ./ vmax).^2)	# convert to slowness squared [s^2/km^2]
        mmax = vec((1f0 ./ vmin).^2)

        ####################################################################################################
        ## FWI using Stochastic gradient decent method
        # Optimization parameters
        # batchsize = 10
        # fhistory_SGD = zeros(Float32, niterations)

        # # Projection operator for bound constraints
        # proj(x) = reshape(median([vec(mmin) vec(x) vec(mmax)]; dims=2), model0.n)

        # # Main loop
        # for j = 1:niterations

        #     # get fwi objective function value and gradient
        #     i = randperm(d_obs.nsrc)[1:batchsize]
        #     fval, gradient = fwi_objective(model0, q[i], d_obs[i]; options=opt)
        #     println("FWI iteration no: ", j, "; function value: ", fval)
        #     fhistory_SGD[j] = fval

        #     # linesearch
        #     step = backtracking_linesearch(model0, q[i], d_obs[i], fval, gradient, proj; alpha=1f0)

        #     # Update model and bound projection
        #     model0.m = proj(model0.m + reshape(step, model0.n))
        # end
        # obj = fhistory_SGD

        ####################################################################################################
        ## FWI using Gauss-Newton method
        # Optimization parameters
        maxiter_GN = 1
        fhistory_GN = zeros(Float32, niterations)

        # Projection operator for bound constraints
        proj(x) = reshape(median([vec(mmin) vec(x) vec(mmax)]; dims=2), model0.n)

        # Gauss-Newton method
        for j=1:niterations

            d_syn = Pr * F0* adjoint(Ps) * q # Synthetic observation data generated on m0
            fval = .5f0*norm(d_syn - d_obs)^2

            println("FWI iteration no: ", j, "; function value: ", fval, "; at ", now() )
            fhistory_GN[j] = fval

            # GN update steps
            J1 = J * Mr
            step =  zeros(Float32, n[1]*n[2])
            lsqr!(step, J1, d_syn - d_obs; maxiter=maxiter_GN, verbose=false)

            # update model and bound constraints
            model0.m = proj(model0.m - reshape(step, model0.n))  # alpha=1
        end
        obj = fhistory_GN

        ####################################################################################################
        fwi = adjoint(model0.m) # inversion result in slowness squared [s^2/km^2]
        # fwi = sqrt.(1f0 ./ adjoint(model0.m)) # inversion result in km/s

        return fwi, obj

    end
end


text:
export AlazarAbortAsyncRead, AlazarAbortCapture, AlazarBeforeAsyncRead
export AlazarBoardsInSystemBySystemID, AlazarBusy, AlazarConfigureAuxIO
export AlazarConfigureLSB, AlazarConfigureRecordAverage, AlazarForceTrigger
export AlazarForceTriggerEnable, AlazarGetBoardBySystemID, AlazarGetBoardKind
export AlazarGetChannelInfo, AlazarInputControl, AlazarNumOfSystems
export AlazarGetParameter, AlazarGetParameterUL
export AlazarPostAsyncBuffer, AlazarRead, AlazarReadEx, AlazarResetTimeStamp
export AlazarSetBWLimit, AlazarSetCaptureClock, AlazarSetExternalClockLevel
export AlazarSetExternalTrigger, AlazarSetLED, AlazarSetParameter
export AlazarSetParameterUL, AlazarSetRecordCount, AlazarSetRecordSize
export AlazarSetTriggerDelay, AlazarSetTriggerOperation, AlazarSetTriggerTimeOut
export AlazarSetTriggerTimeoutTicks, AlazarSleepDevice
export AlazarStartCapture, AlazarTriggered, AlazarWaitAsyncBufferComplete

AlazarAbortAsyncRead(handle::U32) =
    ccall((:AlazarAbortAsyncRead,ats), U32, (U32,), handle)

AlazarAbortCapture(handle::U32) =
    ccall((:AlazarAbortCapture,ats), U32, (U32,), handle)

AlazarBeforeAsyncRead(handle::U32, channels, transferOffset, samplesPerRecord,
        recordsPerBuffer, recordsPerAcquisition, flags) =
    ccall((:AlazarBeforeAsyncRead,ats), U32, (U32, U32, Clong, U32, U32, U32, U32),
        handle, channels, transferOffset, samplesPerRecord, recordsPerBuffer,
        recordsPerAcquisition, flags)

AlazarBoardsInSystemBySystemID(sid::Integer) =
    ccall((:AlazarBoardsInSystemBySystemID,ats),Culong,(Culong,),sid)

AlazarBusy(handle::U32) = ccall((:AlazarBusy,ats), U32, (U32,), handle)

AlazarConfigureAuxIO(handle::U32, mode, parameter) =
    ccall((:AlazarConfigureAuxIO,ats), U32, (U32, U32, U32), handle, mode, parameter)

AlazarConfigureLSB(handle::U32, valueLSB0, valueLSB1) =
    ccall((:AlazarConfigureLSB,ats), U32, (U32, U32, U32), handle, valueLSB0, valueLSB1)

AlazarConfigureRecordAverage(handle::U32, mode, samplesPerRecord,
        recordsPerAverage, options) =
    ccall((:AlazarConfigureRecordAverage,ats), U32, (U32, U32, U32, U32, U32),
        handle, mode, samplesPerRecord, recordsPerAverage, options)

AlazarForceTrigger(handle::U32) =
    ccall((:AlazarForceTrigger,ats),U32,(U32,),handle)

AlazarForceTriggerEnable(handle::U32) =
    ccall((:AlazarForceTriggerEnable,ats),U32,(U32,),handle)

AlazarGetBoardKind(handle::U32) =
    ccall((:AlazarGetBoardKind,ats),U32,(U32,),handle)

AlazarGetBoardBySystemID(systemID, boardID) =
    ccall((:AlazarGetBoardBySystemID,ats),U32,(U32,U32),systemID,boardID)

AlazarGetChannelInfo(handle::U32, memorySize_samples, bitsPerSample) =
    ccall((:AlazarGetChannelInfo,ats), U32, (U32, Ptr{U32}, Ptr{U8}),
        handle, memorySize_samples, bitsPerSample)

AlazarGetParameter(handle::U32, channel, parameter, value) =
    ccall((:AlazarGetParameter,ats), U32, (U32, U8, U32, Ptr{Clong}),
        handle, channel, parameter, value)

AlazarGetParameterUL(handle::U32, channel, parameter, value) =
    ccall((:AlazarGetParameterUL,ats), U32, (U32, U8, U32, Ptr{U32}),
        handle, channel, parameter, value)

AlazarInputControl(handle::U32, channel, coupling, inputRange, impedance) =
    ccall((:AlazarInputControl,ats), U32, (U32, U8, U32, U32, U32),
        handle, channel, coupling, inputRange, impedance)

AlazarNumOfSystems() = ccall((:AlazarNumOfSystems,ats),U32,())

AlazarPostAsyncBuffer(handle::U32, buffer, bufferLength) =
    ccall((:AlazarPostAsyncBuffer,ats), U32, (U32, Ptr{Void}, U32),
        handle, buffer, bufferLength)

AlazarRead(handle::U32, channelId, buffer, elementSize, record,
        transferOffset, transferLength) =
    ccall((:AlazarRead,ats), U32, (U32, U32, Ptr{Void}, Cint, Clong, Clong, U32),
        handle, channelId, buffer, elementSize, record,
        transferOffset, transferLength)

AlazarReadEx(handle::U32, channelId, buffer, elementSize, record,
        transferOffset, transferLength) =
    ccall((:AlazarReadEx,ats), U32, (U32, U32, Ptr{Void}, Cint, Clong, Clonglong, U32),
    handle, channelId, buffer, elementSize, record, transferOffset, transferLength)

AlazarResetTimeStamp(handle::U32, option) =
    ccall((:AlazarResetTimeStamp,ats),U32,(U32,U32),handle, option)

AlazarSetBWLimit(handle::U32, channel, enable) =
    ccall((:AlazarSetBWLimit,ats), U32, (U32, U32, U32), handle, channel, enable)

AlazarSetCaptureClock(handle::U32, source, rate, edge, decimation) =
    ccall((:AlazarSetCaptureClock,ats), U32, (U32, U32, U32, U32, U32),
        handle, source, rate, edge, decimation) #int(source), int(rate), int(edge)

AlazarSetExternalClockLevel(handle::U32, level_percent) =
    ccall((:AlazarSetExternalClockLevel,ats), U32, (U32, Cfloat),
        handle, level_percent)

AlazarSetExternalTrigger(handle::U32, coupling, range) =
    ccall((:AlazarSetExternalTrigger,ats), U32, (U32,U32,U32),
        handle, coupling, range)

AlazarSetLED(handle::U32, ledState) =
    ccall((:AlazarSetLED,ats), U32, (U32, U32), handle, ledState)

AlazarSetParameter(handle::U32, channelId, parameterId, value) =
    ccall((:AlazarSetParameter,ats), U32, (U32, U8, U32, Clong),
        handle, channelId, parameterId, value)

AlazarSetParameterUL(handle::U32, channelId, parameterId, value) =
    ccall((:AlazarSetParameterUL,ats), U32, (U32, U8, U32, U32),
        handle, channelId, parameterId, value)

AlazarSetRecordCount(handle::U32, count) =
    ccode((:AlazarSetRecordCount,ats), U32, (U32, U32), handle, count)

AlazarSetRecordSize(handle::U32, preTriggerSamples, postTriggerSamples) =
    ccall((:AlazarSetRecordSize,ats), U32, (U32, U32, U32),
        handle, preTriggerSamples, postTriggerSamples)

AlazarSetTriggerDelay(handle::U32, delay_samples) =
    ccall((:AlazarSetTriggerDelay,ats),U32,(U32,U32),handle, delay_samples)

AlazarSetTriggerOperation(handle::U32, operation,
        engine1, source1, slope1, level1, engine2, source2, slope2, level2) =
    ccall((:AlazarSetTriggerOperation,ats), U32,
        (U32, U32, U32, U32, U32, U32, U32, U32, U32, U32),
        handle, operation, engine1, source1, slope1, level1,
        engine2, source2, slope2, level2)

AlazarSetTriggerTimeOut(handle::U32, timeout_clocks) =
    ccall((:AlazarSetTriggerTimeOut,ats), U32, (U32, U32), handle, timeout_clocks)

AlazarSleepDevice(handle::U32, sleepState) =
    ccall((:AlazarSleepDevice,ats), U32, (U32, U32), handle, sleepState)

AlazarStartCapture(handle::U32) = ccall((:AlazarStartCapture,ats), U32, (U32,), handle)

AlazarTriggered(handle::U32) = ccall((:AlazarTriggered,ats), U32, (U32,), handle)

AlazarWaitAsyncBufferComplete(handle::U32, buffer, timeout_ms) =
    ccall((:AlazarWaitAsyncBufferComplete,ats), U32, (U32, Ptr{Void}, U32),
        handle, buffer, timeout_ms)


text:
using MINLPTests, JuMP, Ipopt, AmplNLWriter, Test

const OPTIMIZER = MINLPTests.JuMP.with_optimizer(
    AmplNLWriter.Optimizer, Ipopt.amplexe, ["print_level=0"]
)

@testset "MINLPTests" begin
    ###
    ### src/nlp tests.
    ###

    MINLPTests.test_nlp(OPTIMIZER, exclude = [
        "005_011",  # Uses the function `\`
        "006_010",  # User-defined function
        "007_010"   # Infeasible model
    ], objective_tol = 1e-5, primal_tol = 1e-5, dual_tol = NaN)

    @testset "nlp_007_010" begin
        MINLPTests.nlp_007_010(OPTIMIZER, 1e-5, NaN, NaN;
            termination_target = MOI.INFEASIBLE,
            primal_target = MOI.NO_SOLUTION)
    end

    ###
    ### src/nlp-cvx tests.
    ###

    MINLPTests.test_nlp_cvx(OPTIMIZER, exclude = [
        "109_010"  # Ipopt fails to converge
    ])
end


text:
# Simple Linear Regression
# Originally written by Keegan Go for lsqpy
# Translated into LinearLeastSquares.jl by Karanveer Mohan and David Zeng

using LinearLeastSquares
using Gadfly
# Set the random seed to get consistent data
srand(1)

# Number of examples to use
n = 100

# Specify the true value of the variable
true_coeffs = [2; -2; 0.5]

# Generate data
x_data = rand(n) * 5
x_data_expanded = hcat([x_data .^ i for i in 1 : 3]...)
y_data = x_data_expanded * true_coeffs + 0.5 * rand(n)

slope = Variable()
offset = Variable()
line = offset + x_data * slope
residuals = line - y_data
fit_error = sum_squares(residuals)
optval = minimize!(fit_error)

# plot the data and the line
t = [0; 5; 0.1]
p = plot(
  layer(x=x_data, y=y_data, Geom.point),
  layer(x=t, y=evaluate(slope) * t + evaluate(offset), Geom.line),
  Theme(panel_fill=color("white"))
)
# draw(PNG("linear_regression.png", 16cm, 12cm), p)