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

struct LocalNetworkGatewaysApi <: SwaggerApi
    client::Swagger.Client
end

"""

Creates or updates a local network gateway in the specified resource group.
Param: resourceGroupName::String (required)
Param: localNetworkGatewayName::String (required)
Param: parameters::LocalNetworkGateway (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: LocalNetworkGateway
"""
function _swaggerinternal_localNetworkGatewaysCreateOrUpdate(_api::LocalNetworkGatewaysApi, resourceGroupName::String, localNetworkGatewayName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    Swagger.validate_param("localNetworkGatewayName", "localNetworkGatewaysCreateOrUpdate", :minLength, localNetworkGatewayName, 1)

    _ctx = Swagger.Ctx(_api.client, "PUT", LocalNetworkGateway, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/localNetworkGateways/{localNetworkGatewayName}", ["azure_auth"], parameters)
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "localNetworkGatewayName", localNetworkGatewayName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function localNetworkGatewaysCreateOrUpdate(_api::LocalNetworkGatewaysApi, resourceGroupName::String, localNetworkGatewayName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_localNetworkGatewaysCreateOrUpdate(_api, resourceGroupName, localNetworkGatewayName, parameters, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function localNetworkGatewaysCreateOrUpdate(_api::LocalNetworkGatewaysApi, response_stream::Channel, resourceGroupName::String, localNetworkGatewayName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_localNetworkGatewaysCreateOrUpdate(_api, resourceGroupName, localNetworkGatewayName, parameters, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Deletes the specified local network gateway.
Param: resourceGroupName::String (required)
Param: localNetworkGatewayName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: Nothing
"""
function _swaggerinternal_localNetworkGatewaysDelete(_api::LocalNetworkGatewaysApi, resourceGroupName::String, localNetworkGatewayName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    Swagger.validate_param("localNetworkGatewayName", "localNetworkGatewaysDelete", :minLength, localNetworkGatewayName, 1)

    _ctx = Swagger.Ctx(_api.client, "DELETE", Nothing, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/localNetworkGateways/{localNetworkGatewayName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "localNetworkGatewayName", localNetworkGatewayName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function localNetworkGatewaysDelete(_api::LocalNetworkGatewaysApi, resourceGroupName::String, localNetworkGatewayName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_localNetworkGatewaysDelete(_api, resourceGroupName, localNetworkGatewayName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function localNetworkGatewaysDelete(_api::LocalNetworkGatewaysApi, response_stream::Channel, resourceGroupName::String, localNetworkGatewayName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_localNetworkGatewaysDelete(_api, resourceGroupName, localNetworkGatewayName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Gets the specified local network gateway in a resource group.
Param: resourceGroupName::String (required)
Param: localNetworkGatewayName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: LocalNetworkGateway
"""
function _swaggerinternal_localNetworkGatewaysGet(_api::LocalNetworkGatewaysApi, resourceGroupName::String, localNetworkGatewayName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    Swagger.validate_param("localNetworkGatewayName", "localNetworkGatewaysGet", :minLength, localNetworkGatewayName, 1)

    _ctx = Swagger.Ctx(_api.client, "GET", LocalNetworkGateway, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/localNetworkGateways/{localNetworkGatewayName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "localNetworkGatewayName", localNetworkGatewayName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function localNetworkGatewaysGet(_api::LocalNetworkGatewaysApi, resourceGroupName::String, localNetworkGatewayName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_localNetworkGatewaysGet(_api, resourceGroupName, localNetworkGatewayName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function localNetworkGatewaysGet(_api::LocalNetworkGatewaysApi, response_stream::Channel, resourceGroupName::String, localNetworkGatewayName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_localNetworkGatewaysGet(_api, resourceGroupName, localNetworkGatewayName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Gets all the local network gateways in a resource group.
Param: resourceGroupName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: LocalNetworkGatewayListResult
"""
function _swaggerinternal_localNetworkGatewaysList(_api::LocalNetworkGatewaysApi, resourceGroupName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", LocalNetworkGatewayListResult, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/localNetworkGateways", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function localNetworkGatewaysList(_api::LocalNetworkGatewaysApi, resourceGroupName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_localNetworkGatewaysList(_api, resourceGroupName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function localNetworkGatewaysList(_api::LocalNetworkGatewaysApi, response_stream::Channel, resourceGroupName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_localNetworkGatewaysList(_api, resourceGroupName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Updates a local network gateway tags.
Param: resourceGroupName::String (required)
Param: localNetworkGatewayName::String (required)
Param: parameters::TagsObject (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: LocalNetworkGateway
"""
function _swaggerinternal_localNetworkGatewaysUpdateTags(_api::LocalNetworkGatewaysApi, resourceGroupName::String, localNetworkGatewayName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    Swagger.validate_param("localNetworkGatewayName", "localNetworkGatewaysUpdateTags", :minLength, localNetworkGatewayName, 1)

    _ctx = Swagger.Ctx(_api.client, "PATCH", LocalNetworkGateway, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/localNetworkGateways/{localNetworkGatewayName}", ["azure_auth"], parameters)
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "localNetworkGatewayName", localNetworkGatewayName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function localNetworkGatewaysUpdateTags(_api::LocalNetworkGatewaysApi, resourceGroupName::String, localNetworkGatewayName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_localNetworkGatewaysUpdateTags(_api, resourceGroupName, localNetworkGatewayName, parameters, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function localNetworkGatewaysUpdateTags(_api::LocalNetworkGatewaysApi, response_stream::Channel, resourceGroupName::String, localNetworkGatewayName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_localNetworkGatewaysUpdateTags(_api, resourceGroupName, localNetworkGatewayName, parameters, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export localNetworkGatewaysCreateOrUpdate, localNetworkGatewaysDelete, localNetworkGatewaysGet, localNetworkGatewaysList, localNetworkGatewaysUpdateTags


text:
# Autogenerated wrapper script for SoapyLMS7_jll for x86_64-apple-darwin
export LimeUtil, libLMS7Support, libLimeSuite

using soapysdr_jll
using libusb_jll
JLLWrappers.@generate_wrapper_header("SoapyLMS7")
JLLWrappers.@declare_executable_product(LimeUtil)
JLLWrappers.@declare_library_product(libLMS7Support, "libLMS7Support.dylib")
JLLWrappers.@declare_library_product(libLimeSuite, "@rpath/libLimeSuite.20.10-1.dylib")
function __init__()
    JLLWrappers.@generate_init_header(soapysdr_jll, libusb_jll)
    JLLWrappers.@init_executable_product(
        LimeUtil,
        "bin/LimeUtil",
    )

    JLLWrappers.@init_library_product(
        libLMS7Support,
        "lib/SoapySDR/modules0.8/libLMS7Support.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libLimeSuite,
        "lib/libLimeSuite.20.10.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()

ext: jl


text: Dict{String, Any}("answer" => Dict{String, Any}("value" => "\\x64", "type" => "string"))


text:
using NearestNeighbors
using ProgressMeter

push!(LOAD_PATH, @__DIR__)
using PCREO

function unify(bucket, database)
    @assert iszero(last(bucket)[1] - first(bucket)[1])
    a_name = bucket[1][2]
    b_name = bucket[2][2]
    a = database[a_name]
    b = database[b_name]
    if size(a.points) != size(b.points)
        println("NOTE: $a_name and $b_name have identical energy" *
                " but different dimensionality.")
        return nothing
    end
    identical_facets = (sort!(length.(a.facets)) == sort!(length.(b.facets)))
    candidate_isometries = Matrix{Float64}[]
    b_tree = KDTree(b.points)
    dist = minimum(matching_distance(mat * a.points, b_tree)
                   for mat in candidate_isometries)
    @assert 0.0 <= dist
    if dist < 1.0e-10
        println("Unifying $a_name and $b_name.")
        if !identical_facets
            println("WARNING: $a_name and $b_name are isometric but have" *
                    " different facets.")
        end
        a_entries = readdir(joinpath(PCREO_DATABASE_DIRECTORY, a_name))
        b_entries = readdir(joinpath(PCREO_DATABASE_DIRECTORY, b_name))
        if length(a_entries) >= length(b_entries)
            for b_entry in b_entries
                mv(joinpath(PCREO_DATABASE_DIRECTORY, b_name, b_entry),
                   joinpath(PCREO_DATABASE_DIRECTORY, a_name, b_entry))
            end
            rm(joinpath(PCREO_DATABASE_DIRECTORY, b_name))
        else
            for a_entry in a_entries
                mv(joinpath(PCREO_DATABASE_DIRECTORY, a_name, a_entry),
                   joinpath(PCREO_DATABASE_DIRECTORY, b_name, a_entry))
            end
            rm(joinpath(PCREO_DATABASE_DIRECTORY, a_name))
        end
    else
        println("WARNING: $a_name and $b_name have identical energy" *
                " but are not isometric ($dist).")
    end
end

function main()
    while true
        dirnames = readdir(PCREO_DATABASE_DIRECTORY)
        database = Dict{String,PCREORecord}()
        @showprogress for dirname in dirnames
            database[dirname] = PCREORecord(dirname)
        end
        energy_buckets = bucket_by_first(
            sort!([(database[dirname].energy, dirname)
                   for dirname in dirnames]),
            1.0e-10)
        for bucket in energy_buckets
            if length(bucket) > 1
                try
                    unify(bucket, database)
                catch e
                    if e isa AssertionError
                        println("ERROR: ", e)
                    else
                        rethrow(e)
                    end
                end
            end
        end
    end
end

main()


text:
facts("Statistical models") do
  n_weeks = 52
  n_scenarios = 5
  year = 1999
  envflow = .02m^3/s
  maxflow = 5.m^3/s

  name1 = "RivN1"
  name2 = "RivN2"
  name3 = "RivD3"
  name4 = "RivD4"

  nr1 = gen_river(name1, n_weeks, year, n_scenarios)
  nr2 = gen_river(name2, n_weeks, year, n_scenarios)
  dr1 = gen_river(name3, n_weeks, year, n_scenarios, envflow, maxflow)
  dr2 = gen_river(name4, n_weeks, year, n_scenarios, envflow, maxflow)

  TSesNR1 = getScenarios(nr1)

  nrs = [nr1, nr2]
  nrs_names = [name1, name2]
  drs = [dr1, dr2]
  drs_names = [name3, name4]
  rs = [nr1, nr2, dr1, dr2]
  rs_names = [name1, name2, name3, name4]

  r_name = "Reservoir: as simple as possible"
  r_min_capacity = 2.e8m^3
  r_max_capacity = 8.e8m^3

  const a = 1.
  const a_discharge = a * 1.m^3/s
  const a_power = a * 1.W
  const a_height = a * 1.m
  const b = 2.
  const b_discharge = b * 1.m^3/s
  const b_power = b * 1.W
  const c = 3.

  context("Utilities: split an array into random subsets") do
    @fact ReservoirManagement.splitIntoSets_sizes(8, [1.]) --> [8]
    @fact ReservoirManagement.splitIntoSets_sizes(8, [.8, .2]) --> [6, 2]
    @fact_throws splitIntoSets_sizes(8, [1.1]) # Relative sizes do not sum up to 1.

    tses = gen_tses(52, 2000, 8)
    tses_100 = splitIntoSets(tses, [1.])
    @fact length(tses_100) --> 1
    @fact length(tses_100[1]) --> 8
    tses_80_20 = splitIntoSets(tses, [.8, .2])
    @fact length(tses_80_20) --> 2
    @fact length(tses_80_20[1]) --> 6
    @fact length(tses_80_20[2]) --> 2
  end

  context("t-Student: model from weekly data") do
    cis = ReservoirManagement.computeCI_TStudent(nr1, .9, ci_period=Week(1))

    @fact length(cis) --> 3
    for i in 1:3 # mean, low, high
      @fact getPeriod(nr1) --> meta(cis[i])
      @fact length(cis[i]) --> 52
      @fact sum(values(cis[i]) .< 0.) --> 0
    end

    @fact sum(values(cis[2]) .> values(cis[1])) --> 0 # low < mean
    @fact sum(values(cis[2]) .> values(cis[3])) --> 0 # low < high
    @fact sum(values(cis[1]) .> values(cis[3])) --> 0 # mean < high
  end
  #
  # context("Adam2015: model from weekly data") do
  #   cis = computeCI_Adam2015(nr1, .9)
  #
  #   @fact length(cis) --> 3
  #   for i in 1:3 # mean, low, high
  #     @fact getPeriod(nr1) --> meta(cis[i])
  #     @fact length(cis[i]) --> 52
  #     @fact sum(values(cis[i]) .< 0.) --> 0
  #   end
  #
  #   # TODO: not working yet!
  #   # @fact sum(values(cis[2]) .> values(cis[1])) --> 0 # low < mean
  #   # @fact sum(values(cis[2]) .> values(cis[3])) --> 0 # low < high
  #   # @fact sum(values(cis[1]) .> values(cis[3])) --> 0 # mean < high
  # end

  context("Adam2015: scenarios from model (low level interface)") do
    scs = ReservoirManagement.generateScenarios_Adam2015(ReservoirManagement.fit_Adam2015(TSesNR1), Week(1), 2000)
    @fact length(scs) --> 2
    dates = scs[1]
    values = scs[2]
    @fact length(dates) --> length(values)
    @fact length(dates) --> 52
  end

  # context("Adam2015: scenarios from model (high level interface, river)") do
  #   scs = collect(generateScenarios(nr1, 5, Week(1), :Adam2015, Year(1)))
  #   @fact length(scs) --> 5
  #   for s in 1:5
  #     @fact getName(scs[s]) --> getName(nr1)
  #     @fact countScenarios(scs[s]) --> 1
  #     @fact getPeriod(scs[s]) --> Week(1)
  #     @fact countTimeSteps(scs[s]) --> 52
  #     return
  #   end
  # end
end


text:
"""
    getMaxDD(arr)
calculates the maximum Drawdown of profit and loss series
(not the aggregated profit and loss series)
be aware: data must be sorted by date upfront!
"""
function getMaxDD(profitLossArray::Array{Float64,1})::Float64
    peak=0.0
    maxDD=0.0
    equity = 0.0
    aggregatedSeries = accumulatedProfitAndLossSeries=accumulate(+, profitLossArray);

    for n in eachindex(aggregatedSeries)
        equity = aggregatedSeries[n]
        #Check if new peak is available
        if equity>peak
            peak = equity
        else
            if abs(peak-equity) > maxDD
                maxDD = abs(peak - equity)
            end
        end
    end
    return maxDD
end

"""
    getMaxDDRate_relativeToLastPeak(arr)
calculates maximum Drawdown rate (not percentage!)
relative to the last Peak before DD (not relative to initial Capital)
parameter arr is the profit and loss series (not the aggregated profit and loss series)
be aware: data must be sorted by date upfront!
be aware: returns inf if there is no DD *only positiv trades*
          returns -inf if there are no positive trades
"""
function getMaxDDRate_relativeToLastPeak(arr::Array{Float64,1}, )::Float64
    equity = 0.0
    peak=0.0
    last_peak_before_maxDD=0.0
    last_peak_before_maxDD_position = 0
    maxDD=0.0                                                                   # initialize maxDD with Type FLOAT
    maxDDRate=0.0
    aggregatedSeries = getAccumulatedProfitLossSeries(arr);                     # create aggregatedSeries

    for n in eachindex(aggregatedSeries)
        equity = aggregatedSeries[n]
        #Check if new peak is available
        if equity>peak
            peak = equity
        elseif abs(peak-equity) > maxDD
            maxDD = abs(peak - equity)
            last_peak_before_maxDD_position = n
        end
    end

    #if no dd don´t search for last peak!
    if maxDD > 0
        if maximum(aggregatedSeries[1:last_peak_before_maxDD_position]) > 0
            last_peak_before_maxDD = maximum(aggregatedSeries[1:last_peak_before_maxDD_position])
        else
            last_peak_before_maxDD = 0.0
        end
    else
        last_peak_before_maxDD = 0.0
    end

    println("LastPeak:", last_peak_before_maxDD)

    if maxDD == 0.0 && last_peak_before_maxDD== 0.0                             #Assume there are only positiv trades
        return Inf
    elseif maxDD > 0 && last_peak_before_maxDD == 0.0                           #Assume there are only negativ trades
        return -Inf
    else
        return maxDD/last_peak_before_maxDD
    end
end

"""
    getMaxDD_rate
sets the maximum Drawdown in relation to invested capital
"""
function getMaxDD_Rate(profitLossSeriesArray::Array{Float64,1}, initialCapital::Float64):: Float64
    maxDD = getMaxDD(profitLossSeriesArray)
    return (maxDD/initialCapital)
end

"""
    getMaxDD_percentage
sets the maximum Drawdown in relation to invested capital
"""
function getMaxDDPercentage(profitLossSeriesArray::Array{Float64,1}, initialCapital::Float64):: Float64
    maxDD = getMaxDD(profitLossSeriesArray)
    return (maxDD/initialCapital)*100
end

max_stars_repo_name: RainerNeu/tradeStatistics

max_issues_repo_name: RainerNeu/tradeStatistics

max_forks_repo_name: RainerNeu/tradeStatistics


text:
function calFlxCorrection!(
    ocn :: Ocean;
    τ   :: Float64 = 10 * 86400.0,
    cfgs...
)
    do_convadjust = cfgs[:do_convadjust]

    Δt = cfgs[:Δt]
    r = Δt / τ
    rr = r / (1.0 + r)

    calTsSsMixed!(ocn)

    # Euler backward method
    @. ΔT_mixed = rr * (ocn.in_fld.Tclim - ocn.Ts_mixed)
    @. ΔS_mixed = rr * (ocn.in_fld.Sclim - ocn.Ss_mixed)

    @loop_hor ocn i j let

        T_ML = ocn.T_ML[i, j]
        S_ML = ocn.S_ML[i, j]
        FLDO = ocn.FLDO[i, j]
        h_ML = ocn.h_ML[i, j]

        if FLDO == -1  # Whole ocean layer. In the future I want to restrict FLDO so that code can be cleaner.

            ΔT_ML = 0.0
            for k=1:ocn.Nz[i, j]
                ΔT_ML += ΔT_mixed[k, i, j] * ocn.hs[k, i, j]
            end
            ΔT_ML /= h_ML

            ΔS_ML = 0.0
            for k=1:ocn.Nz[i, j]
                ΔS_ML += ΔS_mixed[k, i, j] * ocn.hs[k, i, j]
            end
            ΔS_ML /= h_ML

            T_ML += ΔT_ML
            S_ML += ΔS_ML
            ocn.T_ML[i, j] = T_ML
            ocn.S_ML[i, j] = S_ML

            ocn.Ts[:, i, j] .= T_ML
            ocn.Ss[:, i, j] .= S_ML

        else
            # With in mixed-layer
            if FLDO == 1
                ocn.T_ML += ΔT_mixed[1, i, j]
                ocn.S_ML += ΔS_mixed[1, i, j]
            elseif FLDO > 1

                ΔT_ML = 0.0
                for k=1:FLDO-1
                    ΔT_ML += ΔT_mixed[k, i, j] * ocn.hs[k, i, j]
                end
                ΔT_ML += ΔT_mixed[FLDO, k, j] * (ocn.zs[FLDO, i, j] + h_ML)
                ΔT_ML /= h_ML

                ΔS_ML = 0.0
                for k=1:FLDO-1
                    ΔS_ML += ΔS_mixed[k, i, j] * ocn.hs[k, i, j]
                end
                ΔS_ML += ΔS_mixed[FLDO, k, j] * (ocn.zs[FLDO, i, j] + h_ML)
                ΔS_ML /= h_ML

                T_ML += ΔT_ML
                S_ML += ΔS_ML
                ocn.T_ML[i, j] = T_ML
                ocn.S_ML[i, j] = S_ML

                ocn.Ts[1:FLDO-1, i, j] .= T_ML
                ocn.Ss[1:FLDO-1, i, j] .= S_ML

            end

            # Deeper ocean
            if FLDO != -1
                @. ocn.Ts[FLDO:end, i, j] += ΔT_mixed[FLDO:end, i, j]
                @. ocn.Ss[FLDO:end, i, j] += ΔS_mixed[FLDO:end, i, j]
            end

        end

        OC_updateB!(ocn, i, j)

    end

    @. ocn.qflx_T_correction = ΔT_mixed * ocn.hs * ρc_sw   / Δt   # + => warming
    @. ocn.qflx_S_correction = ΔS_mixed * ocn.hs           / Δt   # + => saltier

    if do_convadjust
        @loop_hor ocn i j let
            OC_doConvectiveAdjustment!(ocn, i, j;)
        end
    end

end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
module ArrayFire

export AFArray
using Cxx
#using Base.Meta

import Cxx: CppEnum

addHeaderDir("/usr/local/include", kind=C_System)
Libdl.dlopen("/usr/local/lib/libafcpu.so", Libdl.RTLD_GLOBAL)
cxxinclude("arrayfire.h")

setDevice(i) = icxx"af::setDevice($i);"

const f32 = CppEnum{:af_dtype}(0)
const c32 = CppEnum{:af_dtype}(1)
const f64 = CppEnum{:af_dtype}(2)
const c64 = CppEnum{:af_dtype}(3)
const b8  = CppEnum{:af_dtype}(4)
const s32 = CppEnum{:af_dtype}(5)
const u32 = CppEnum{:af_dtype}(6)
const u8  = CppEnum{:af_dtype}(7)
const s64 = CppEnum{:af_dtype}(8)
const u64 = CppEnum{:af_dtype}(9)

aftype(::Type{Float32}) = f32
aftype(::Type{Complex{Float32}}) = c32
aftype(::Type{Float64}) = f64
aftype(::Type{Complex{Float64}}) = c64
aftype(::Type{Bool}) = b8
aftype(::Type{Int32}) = s32
aftype(::Type{UInt32}) = u32
aftype(::Type{UInt8}) = u8
aftype(::Type{Int64}) = s64
aftype(::Type{UInt64}) = u64

include("array.jl")

end


text:
using Test
using DFTK

include("testcases.jl")

@testset "Compare different SCF algorithms" begin
    Ecut = 3
    n_bands = 6
    fft_size = [9, 9, 9]
    tol = 1e-7

    Si = ElementPsp(silicon.atnum, psp=load_psp(silicon.psp))
    model = model_DFT(silicon.lattice, [Si => silicon.positions], [:lda_xc_teter93])
    basis = PlaneWaveBasis(model, Ecut, silicon.kcoords, silicon.ksymops; fft_size=fft_size)

    # Run nlsolve without guess
    ρ0 = RealFourierArray(basis)  # RFA of zeros
    scfres = self_consistent_field(basis; ρ=ρ0, tol=tol)
    ρ_nl = scfres.ρ.fourier

    # Run DM
    println("\nTesting direct minimization")
    dmres = direct_minimization(basis; g_tol=tol)
    ρ_dm = dmres.ρ.fourier
    @test maximum(abs.(ρ_dm - ρ_nl)) < sqrt(tol) / 10

    # Run other SCFs with SAD guess
    ρ0 = guess_density(basis, [Si => silicon.positions])
    for solver in (scf_nlsolve_solver, scf_damping_solver, scf_anderson_solver,
                   scf_CROP_solver)
        println("\nTesting $solver")
        scfres = self_consistent_field(basis; ρ=ρ0, solver=solver(), tol=tol)
        ρ_alg = scfres.ρ.fourier
        @test maximum(abs.(ρ_alg - ρ_nl)) < sqrt(tol) / 10
    end

    # Run other mixing with default solver (the others are too slow...)
    for mixing in (KerkerMixing(), SimpleMixing(), SimpleMixing(.5))
        println("\n Testing $mixing")
        scfres = self_consistent_field(basis; ρ=ρ0, mixing=mixing, tol=tol)
        ρ_alg = scfres.ρ.fourier
        @test maximum(abs.(ρ_alg - ρ_nl)) < sqrt(tol) / 10
    end
end


text:
using Test
using StaticArrays
using PolygonArea
using PolygonArea: Point, rotate, translate

@testset "Points" begin
    p = Point(1.0, 0.0)
    @test typeof(p) == Point{Float64}
    @test rotate(p, π/2) ≈ Point(0.0, 1.0)
    @test rotate(p, π/2, center=Point(2.0, 0.0)) ≈ Point(2.0, -1.0)
    @test rotate(p, π/4, center=(1.0, 0.0)) == Point(1.0, 0.0)
    @test translate(p, SVector(1.0, 0.0)) ≈ Point(2.0, 0.0)

    p = Point(1//2, 3//2)
    @test typeof(p) ==  Point{Rational{Int}}
    @test translate(p, SVector(1//1, 0//1)) == Point(3//2, 3//2)
end


text:
"""
Dictionary types which keys are fixed at creation time
"""
abstract AbstractFixedSizeDict{Keys}

"""
Dictionary types which keys and values are fixed at creation time
"""
immutable FixedKeyValueDict{Keys, Values<:Tuple} <: AbstractFixedSizeDict{Keys}
    keys::Keys
    values::Values
end
"""
Dictionary types which keys are fixed at creation time
"""
immutable FixedKeyDict{Keys, Values<:AbstractVector} <: AbstractFixedSizeDict{Keys}
    keys::Keys
    values::Values
end

"""
Generic constructor that dispatches XDict(pair, pair, ...) to XDict((pair, pair))
"""
function call{T<:AbstractFixedSizeDict}(::Type{T}, args...)
    T(args)
end

# TODO set operations

"""
Constructor for a list of keys together with a list of values which should have the same length
"""
function FixedKeyValueDict{N}(keys::NTuple{N, Symbol}, values::NTuple{N, Any})
    # TODO enforce unique keys?!
    keys = map(k->Val{k}(), keys)
    FixedKeyValueDict(keys, values)
end

"""
Constructor for a list of pairs of key => value.
Arbitrary data structures of length 2 can be used
"""
function FixedKeyValueDict(key_values)
    kv = tuple(key_values...) # needs to be a tuple, since both key and value needs to be a tuple
    keys = map(k->Val{first(k)}(), kv)
    values = map(last, kv)
    FixedKeyValueDict(keys, values)
end

"""
Constructor for a list of keys together with a list of values which should have the same length
"""
function FixedKeyDict{N}(keys::NTuple{N, Symbol}, values::AbstractVector)
    keys = map(k->Val{k}(), keys)
    FixedKeyDict(keys, values)
end

"""
Constructor for a list of pairs of key => value.
Arbitrary data structures of length 2 can be used
"""
function FixedKeyDict(key_values)
    kv = collect(key_values) # needs collect for some datastructures like Dict
    keys = ntuple(length(kv)) do i # ntuple, to always get a tuple
        Val{first(kv[i])}()
    end
    values = [v for (k,v) in kv] # for-comprehension to always get a Vector
    FixedKeyDict(keys, values)
end

keys{T<:AbstractFixedSizeDict}(fsd::T) = fsd.keys
# this is a generated function, since T.parameters results in the code to be slow otherwise
@generated function keys{T<:AbstractFixedSizeDict}(::Type{T})
    ks = map(x-> x(), T.parameters[1].parameters)
    :(tuple($(ks...)))
end
# use @generated to do field lookup at compile time instead of runtime
@generated function key2index{Keys, Key}(x::AbstractFixedSizeDict{Keys}, k::Val{Key})
    index = findfirst(Keys.parameters, k)
    index == 0 && throw(KeyError("key $Key not found in $x"))
    :($index)
end
# also a generated function, since we can inline the result for every type, which
# doesn't happen otherwise
@generated function haskey{Keys, Key}(sd::AbstractFixedSizeDict{Keys}, ::Key)
    :($(Key in Keys.parameters))
end

values(x::AbstractFixedSizeDict) = x.values
length(x::AbstractFixedSizeDict) = length(x.values)

function Base.start(x::AbstractFixedSizeDict)
    (1, keys(x)) # we pass around the keys so that we don't have to get them multiple times
end
function Base.next(x::AbstractFixedSizeDict, state)
    index, ks = state
    (ks[index], x.values[index]), (index+1, ks)
end
function Base.done(x::AbstractFixedSizeDict, state)
    state[1] > length(x)
end

# k should be typed, but that gives weird behaviour on 0.5
function getindex{Keys, Key}(sd::AbstractFixedSizeDict{Keys}, k::Val{Key})
    @inbounds return sd.values[key2index(sd, k)]
end

function setindex!{Keys, Key}(sd::FixedKeyDict{Keys}, value, k::Val{Key})
    @inbounds return sd.values[key2index(sd, k)] = value
end

function get{Keys}(f::Function, sd::AbstractFixedSizeDict{Keys}, k)
    if haskey(sd, k)
        return sd[k]
    end
    f()
end
function get{Keys}(sd::AbstractFixedSizeDict{Keys}, k, default)
    if haskey(sd, k)
        return sd[k]
    end
    default
end

"""
generates the expression to acces a field of a dict via a Val{Symbol}
"""
function getfield_expr(dict, key)
    :($dict[Val{$key}()])
end

"""
Allows getfield like access to the keys of a FixedDict
"""
macro get(expr)
    if expr.head == :(.)
        return esc(getfield_expr(expr.args...))
    elseif expr.head == :(=)
        fieldexpr, val = expr.args
        getfieldex = getfield_expr(fieldexpr.args...)
        return esc(Expr(:(=), getfieldex, val))
    else
        quote
            errored_expr = $(string(expr))
            throw(ArgumentError(
                "Expression of the form $errored_expr not allowed. Try a.key, or a.key=value"
            ))
        end
    end
end

# we need to redefine equality, since dicts with the same set of keys and values should
# match regardless of order
function ==(a::AbstractFixedSizeDict, b::AbstractFixedSizeDict)
    length(a) == length(b) || return false # if not same length, can't be equal
    for (key, value) in a
        if haskey(b, key)
            value == b[key] || return false # if not same value, they no equal!
        else
            return false
        end
    end
    true
end


text:
@testset "likelihood.jl" begin
    println("             > likelihood.jl")
    @testset "consistency" begin
        using BinomialSynapses: likelihood_indices

        m_out = 1024
        m_in  = 1024
        Ns = CuArray(rand(1:128, m_out))
        ps = CUDA.rand(m_out)
        qs = CUDA.rand(m_out)
        σs = CUDA.rand(m_out)
        τs = CUDA.rand(m_out)
        model = BinomialModel(Ns, ps, qs, σs, τs);

        ns = CuArray(rand(1:128, m_out, m_in))
        ks = CUDA.zeros(Int, m_out, m_in);

        ks .= ns .% 64

        u1      = likelihood(ks, model, 0.3f0)
        u2, idx = likelihood_indices(ks, model, 0.3f0)

        @test u1 ≈ u2
        @test maximum(idx) <= m_in
        @test minimum(idx) >= 1

        idx = Array(idx)
        for r in eachrow(idx)
            @test issorted(r)
        end

        if RUN_BENCHMARKS
            println("")
            println("Benchmarking function likelihood!: should take about 300μs")
            display(@benchmark CUDA.@sync likelihood($ks, $model, 0.3f0))
            println("")
            println("")
            println("Benchmarking function likelihood_resample!: should take about 4ms")
            state = BinomialState(ns, ks)
            obs   = BinomialObservation(0.3f0, 0.1f0)
            display(@benchmark CUDA.@sync likelihood_resample!($state, $model, $obs))
            println("")
        end
    end
    @testset "correctness of values" begin
        observation = 3.0
        M_out = 6
        M_in = 4
        Ns = 10 .* CUDA.ones(Int, M_out)
        ps = CUDA.rand(M_out)
        qs = cu([1f0, 1f0, 3f0, 3f0, 1f0, 3f0])
        σs = Float32(0.1) .* CUDA.ones(M_out)
        τs = CUDA.rand(M_out)
        model = BinomialModel(Ns, ps, qs, σs, τs);

        ks = cu([3 3 3 3; 1 1 1 1; 1 1 1 1; 3 3 3 3; 3 1 1 1; 1 3 3 3])

        u, idx = likelihood_indices(ks, model, observation)
        u = Array(u)
        idx = Array(idx)

        # Combinations of q and k that correspond to the observation should have
        # a high likelihood
        @test u[1] > 3.9894
        @test u[3] > 3.9894

        # Combinations of q and k that do not match the observation should have
        # a low likelihood
        @test u[2] < 1.5f-10
        @test u[4] < 1.5f-10

        # Resampling should pick the particles that match q and the observation
        @test all(idx[5:6,:].==1)
    end
end


text:
function realization(ngrid,H,Br,m,y,bx,by,pd,q)

# computes the velocity and acceleration on a lat,long grid
# at a series of snapshots.
#
# input
#
# ngrid -  number of grid points in colatitude
# H    - layer thickness (km)
# Np   - dimensionless buoyancy frequency (N/Omega)
# Br   - radial magnetic field  (mT)
# m    - angular order m
# y    - meridional coordinate
# bx,by - x and y components of magnetic perturbation
# pd    - period  (years)
# q     - quality factor
#
# output
#
# tvec      - array to times to evaluate velocity
# lat,long  -  array of latitude and longitude
# vy,vx     -  theta and phi components of velocity at CMB
# ay,ax    -  theta and phi components of acceleration

# array of times
dt = 0.5;
tvec = (0.0 : dt : abs(pd));

# latitude and longitude grid
#ngrid = 100;
dlat = pi / ngrid;
dlong = pi / ngrid;
lat = (-pi/2+dlat/2 : dlat : pi/2-dlat/2);
long = (dlong/2 : dlong : 2 * pi-dlong/2);
colat = (dlat/2 : dlat : pi-dlat/2);  # reverse order

# allocate memory for velocity and acceleration
vx = Array{Float64,3}(undef,length(tvec),ngrid,2*ngrid);
vy = Array{Float64,3}(undef,length(tvec),ngrid,2*ngrid);
ax = Array{Float64,3}(undef,length(tvec),ngrid,2*ngrid);
ay = Array{Float64,3}(undef,length(tvec),ngrid,2*ngrid);

# loop over times
for k = 1 : length(tvec)

   tstep = tvec[k];
   vxt,vyt,axt,ayt = velocity(H,Br,m,y,lat,long,bx,by,pd,q,tstep);
   vx[k,:,:]=vxt;
   vy[k,:,:]=vyt;
   ax[k,:,:]=axt;
   ay[k,:,:]=ayt;

end

# collect into arrays
t = collect(tvec);
clt = collect(colat);
lg = collect(long);

return t,clt,lg,vx,vy,ax,ay

end

max_stars_repo_name: bbuffett/MagRossby

max_issues_repo_name: bbuffett/MagRossby

max_forks_repo_name: bbuffett/MagRossby


text:
# primitive fitting API

"""
An abstract type that supertypes all the fitted shapes.
"""
abstract type FittedShape end

"""
Return the default parameters as a `NamedTuple` for the passed `FittedShape`.

# Implementation
Example definition for `MyShape<:FittedShape`:
`defaultshapeparameters(::Type{MyShape})=(myshape=(ϵ=1, α=deg2rad(5),),)`.
"""
function defaultshapeparameters end

"""
Fit a primitive shape to point-normal pairs
(passed as an array of points and an array of normals).
The type of the primitive and a parameter named tuple is also passed.
The RANSACCloud is also passed if needed for more advanced methods.
Return `nothing`, if can't fit the shape to the points.

# Implementation
Signature: `fit(::Type{MyShape}, p, n, pc, params)`.

It should return `nothing` if the fit wasn't succesfull, `MyShape` otherwise.
If the fit results multiple shapes, it can return an array of `MyShape`s as well.
"""
function fit end

"""
Compute the score of a candidate and return its score and the corresponding points.

# Implementation
Signature: `scorecandidate(pc, candidate::MyShape, subsetID, params)`.

`subsetID` is the index of the subset that is used to estimate the score.
You must count the compatible points (that are enabled) in the given subset.
You can get the points by: `pc.vertices[pc.subsets[subsetID]]` and normals similarly.
Use the [`estimatescore`](@ref) function to estimate a score,
then return a  `(score, inpoints)`,
where `inpoints` is the indexes of the points, that are counted.
"""
function scorecandidate end

"""
Refit a primitive to the whole point cloud, say search for all the compatible points
in the point cloud.

# Implementation
Signature: `refit(s::T, pc, params) where {T<:MyShape}`.

Search all the compatible (and enabled) points in the whole point cloud.
After finding the indexes, return an [`ExtractedShape`](@ref).
"""
function refit end

"""
Return a string that tells the "human-readable" type
(plane, spehere, etc.) of a `FittedShape`.

# Implementation
`strt(s::MyShape) = "myshape"`
"""
function strt end

# functions to work with primitives

"""
    ExtractedShape{S<:FittedShape}

Store an extraced primitive (`FittedShape`)
and the points that belong to the shape as `Vector{Int}`.

# Implementation

When constructing in a `refit` function, just call the constructor with:
`ExtractedShape(shape, compatible_points)`.
"""
struct ExtractedShape{S<:FittedShape}
    shape::S
    inpoints::Vector{Int}
end

Base.show(io::IO, x::ExtractedShape{A}) where {A} =
    print(io, "Cand: (", x.shape, "), $(length(x.inpoints)) ps")

"""
    struct IterationCandidates

Store shapes, their score and their points in a struct of arrays.
"""
struct IterationCandidates
    shapes::Vector{FittedShape}
    scores::Vector{ConfidenceInterval}
    inpoints::Vector{Vector{Int}}
end

function IterationCandidates()
    return IterationCandidates(FittedShape[], ConfidenceInterval[], Vector{Vector{Int}}())
end

Base.length(x::IterationCandidates) = size(x.shapes, 1)

Base.show(io::IO, x::IterationCandidates) =
    print(io, "IterationCandidates: $(size(x.shapes,1)) candidates")

"""
    recordscore!(ic::IterationCandidates, shape, score, inpoints)

Record the score and inpoints of a shape to an `IterationCandidates`.
"""
function recordscore!(ic::IterationCandidates, shape, score, inpoints)
    push!(ic.shapes, shape)
    push!(ic.scores, score)
    push!(ic.inpoints, inpoints)
    return ic
end

"""
    deleteat!(ic::IterationCandidates, args...)

Forward `deleteat!` to field arrays.
"""
function deleteat!(ic::IterationCandidates, arg)
    deleteat!(ic.shapes, arg)
    deleteat!(ic.scores, arg)
    deleteat!(ic.inpoints, arg)
    return ic
end

"""
    findhighestscore(A::IterationCandidates)

Find the largest expected value in an `IterationCandidates`.

Indicate if there's an overlap.
"""
function findhighestscore(A::IterationCandidates)
    (length(A) > 0) || return (index = 0, overlap = false)
    ind = 1
    scores = A.scores
    highest = E(scores[1])
    for i in eachindex(scores)
        esc = E(scores[i])
        if esc > highest
            highest = esc
            ind = i
        end
    end

    for i in eachindex(scores)
        i == ind && continue
        isoverlap(scores[i], scores[ind]) && return (index = ind, overlap = true)
    end
    return (index = ind, overlap = false)
end

"""
    forcefitshapes!(points, normals, parameters, candidates, level_array, octree_lev, pc)

Call `fit` based on the given parameters and add valid candidates to
"""
function forcefitshapes!(points, normals, parameters, candidates, level_array, octree_lev, pc)
    @extract parameters.iteration : shape_types
    for s in shape_types
        fitted = fit(s, points, normals, pc, parameters)
        fitted === nothing && continue
        push2candidatesandlevels!(candidates, fitted, level_array, octree_lev)
    end
    return nothing
end

"""
    scorecandidates!(pc, iterationcandidates, candidates, subsetID, params, octree_levels)

Call [`scorecandidate`](@ref) for all candidates.
Reset candidate, and octree level lists (`candidates` and `octree_levels` respectively).
"""
function scorecandidates!(pc, iterationcandidates, candidates, subsetID, params, octree_levels)
    for i in eachindex(candidates)
        sc, ip = scorecandidate(pc, candidates[i], subsetID, params)
        pc.levelscore[octree_levels[i]] += E(sc)
        recordscore!(iterationcandidates, candidates[i], sc, ip)
    end
    empty!(candidates)
    empty!(octree_levels)
    return nothing
end

"""
    invalidate_indexes!(pc, indexlist)

Invalidate indexes of the to-be-extracted shape's points, based on the `indexlist`.
"""
function invalidate_indexes!(pc, indexlist)
    for a in indexlist
        pc.isenabled[a] = false
    end
    return nothing
end

"""
    removeinvalidshapes!(pc, candidates::IterationCandidates)

Remove shapes from `candidates`, which are invalid, say have points that are not enabled.
"""
function removeinvalidshapes!(pc, candidates::IterationCandidates)
    toremove = Int[]
    for i in eachindex(candidates.inpoints)
        for a in candidates.inpoints[i]
            if ! pc.isenabled[a]
                push!(toremove, i)
                break
            end
        end
    end
    deleteat!(candidates, toremove)
    return nothing
end

function samplepointcloud!(indexv, pointcloud, octree, params)
    @extract params : iter_params=iteration
    @extract iter_params : drawN minsubsetN prob_det τ

    r1 = rand(1:pointcloud.size)

    #TODO: helyettesíteni valami okosabbal,
    # pl mindig az első enabled - ha már nagyon sok ki van véve,
    # akkor az gyorsabb lesz
    while ! pointcloud.isenabled[r1]
        r1 = rand(1:pointcloud.size)
    end
    # search for r1 in octree
    current_leaf = findleaf(octree, pointcloud.vertices[r1])
    # get the depth
    max_depth = current_leaf.data.depth
    # get the best level index == best depth
    curr_level_i = argmax(pointcloud.levelweight[1:max_depth])
    # get the best cell, based on the index
    curr_cell = getnthcell(current_leaf, curr_level_i)
    # an indexer array for random indexing
    cell_ind = curr_cell.data.incellpoints
    # frome the above, those that are enabled
    bool_cell_ind = @view pointcloud.isenabled[cell_ind]
    enabled_inds = cell_ind[bool_cell_ind]
    # if there's less enabled vertice than needed, skip the rest
    size(enabled_inds, 1) < drawN && return (false, 0)
    indexv[1] = r1

    # made up heuristic
    if size(enabled_inds, 1) < 20*drawN
        # if there's few randoms, then just choose the first ones
        # random index
        sel = 0
        # sd index
        seli = 2
        while true
            sel += 1
            # don't choose the first one
            enabled_inds[sel] == indexv[1] && continue
            indexv[seli] = enabled_inds[sel]
            seli += 1
            seli == drawN+1 && return (false, 1)
        end
        route = 1
    else
        # if there's enough random, randomly select
        for idk in 2:drawN
            nexti = rand(1:size(enabled_inds, 1))
            if indexv[1] == enabled_inds[nexti]
                # try oncemore
                nexti = rand(1:size(enabled_inds, 1))
            end
            #TODO: check if the same
            #TODO: do something with it
            indexv[idk] = enabled_inds[nexti]
        end
        route = 2
    end

    if ! allisdifferent(indexv)
        @logmsg IterLow1 "Selected indexes have same element: $indexv; route $route was taken."
        return (false, 1)
    end
    return (true, curr_level_i)
end

function samplepointcloud2!(indexv, pointcloud, octree, params)
    @extract params : iter_params=iteration
    @extract iter_params : drawN minsubsetN prob_det τ
    @assert drawN > 1 "At least 2 point is expected in this sampling function!"

    r1 = rand(1:pointcloud.size)

    #TODO: helyettesíteni valami okosabbal,
    # pl mindig az első enabled - ha már nagyon sok ki van véve,
    # akkor az gyorsabb lesz
    while ! pointcloud.isenabled[r1]
        r1 = rand(1:pointcloud.size)
    end
    # search for r1 in octree
    current_leaf = findleaf(octree, pointcloud.vertices[r1])
    # get the depth
    max_depth = current_leaf.data.depth
    # get the best level index == best depth
    curr_level_i = argmax(pointcloud.levelweight[1:max_depth])
    # get the best cell, based on the index
    curr_cell = getnthcell(current_leaf, curr_level_i)
    # an indexer array for random indexing
    cell_ind = curr_cell.data.incellpoints
    # frome the above, those that are enabled
    bool_cell_ind = @view pointcloud.isenabled[cell_ind]
    count(bool_cell_ind) < drawN && return (false, 0)

    enabled_inds = cell_ind[bool_cell_ind]
    # if there's less enabled vertice than needed, skip the rest
    indexv[1] = r1

    # if there's enough random, randomly select
    for idk in 2:drawN
        nexti = rand(1:size(enabled_inds, 1))
        if indexv[1] == enabled_inds[nexti]
            # try oncemore
            nexti = rand(1:size(enabled_inds, 1))
        end
        #TODO: check if the same
        #TODO: do something with it
        indexv[idk] = enabled_inds[nexti]
    end

    if ! allisdifferent(indexv)
        @logmsg IterLow1 "Selected indexes have same element: $indexv."
        return (false, 1)
    end
    return (true, curr_level_i)
end

function samplepointcloud3!(indexv, pointcloud, octree, params)
    @extract params : iter_params=iteration
    @extract iter_params : drawN minsubsetN prob_det τ
    @assert drawN > 1 "At least 2 point is expected in this sampling function!"

    r1 = rand(1:pointcloud.size)

    #TODO: helyettesíteni valami okosabbal,
    # pl mindig az első enabled - ha már nagyon sok ki van véve,
    # akkor az gyorsabb lesz
    while ! pointcloud.isenabled[r1]
        r1 = rand(1:pointcloud.size)
    end
    # search for r1 in octree
    current_leaf = findleaf(octree, pointcloud.vertices[r1])
    # get the depth
    max_depth = current_leaf.data.depth
    # get the best level index == best depth
    curr_level_i = argmax(pointcloud.levelweight[1:max_depth])
    # get the best cell, based on the index
    curr_cell = getnthcell(current_leaf, curr_level_i)
    # an indexer array for random indexing
    cell_ind = curr_cell.data.incellpoints
    # frome the above, those that are enabled
    bool_cell_ind = @view pointcloud.isenabled[cell_ind]
    count(bool_cell_ind) < drawN && return (false, 0)

    enabled_inds = cell_ind[bool_cell_ind]
    # if there's less enabled vertice than needed, skip the rest
    indexv[1] = r1

    # if there's enough random, randomly select
    indexv[2:end] = rand(enabled_inds, drawN-1)

    if ! allisdifferent(indexv)
        @logmsg IterLow1 "Selected indexes have same element: $indexv."
        return (false, 1)
    end
    return (true, curr_level_i)
end

function samplepointcloud4!(indexv, pointcloud, params)
    @extract params : iter_params=iteration
    @extract iter_params : drawN minsubsetN prob_det τ
    @assert drawN > 1 "At least 2 point is expected in this sampling function!"

    r1 = rand(1:pointcloud.size)

    #TODO: helyettesíteni valami okosabbal,
    # pl mindig az első enabled - ha már nagyon sok ki van véve,
    # akkor az gyorsabb lesz
    while ! pointcloud.isenabled[r1]
        r1 = rand(1:pointcloud.size)
    end
    # search for r1 in octree
    current_leaf = findleaf(pointcloud.octree, pointcloud.vertices[r1])
    # get the depth
    max_depth = current_leaf.data.depth
    # get the best level index == best depth
    curr_level_i = argmax(pointcloud.levelweight[1:max_depth])
    # get the best cell, based on the index
    curr_cell = getnthcell(current_leaf, curr_level_i)
    # an indexer array for random indexing
    cell_ind = curr_cell.data.incellpoints
    # frome the above, those that are enabled
    enabled_inds = @view cell_ind[pointcloud.isenabled[cell_ind]]
    size(enabled_inds, 1) < drawN && return (false, 0)

    # if there's less enabled vertice than needed, skip the rest
    indexv[1] = r1

    # if there's enough random, randomly select
    for idk in 2:drawN
        nexti = rand(1:size(enabled_inds, 1))
        if indexv[1] == enabled_inds[nexti]
            # try oncemore
            nexti = rand(1:size(enabled_inds, 1))
        end
        #TODO: check if the same
        #TODO: do something with it
        indexv[idk] = enabled_inds[nexti]
    end

    if ! allisdifferent(indexv)
        @logmsg IterLow1 "Selected indexes have same element: $indexv."
        return (false, 1)
    end
    return (true, curr_level_i)
end


text:
module test_octo_pretty

using Test
using Octo.Pretty

nts = [(id=1, name="John"), (id=2, name="Tom")]
@test Pretty.table(nts) == """
|   id | name   |
| ---- | ------ |
|    1 | John   |
|    2 | Tom    |
2 rows."""

Pretty.set(colsize=30)
nts = [(id=1, case="글씨가 넘흐 길어서 짤릴 때"),
       (id=2, case="글씨에 점  하나가"),
       (id=3, case="글씨가 거의 비슷해"),
       (id=4, case="안짤리는 넘")]
@test Pretty.table(nts) == """
|   id | case              |
| ---- | ----------------- |
|    1 | 글씨가 넘흐 길... |
|    2 | 글씨에 점  하나가 |
|    3 | 글씨가 거의 비... |
|    4 | 안짤리는 넘       |
4 rows."""

nts = Vector{NamedTuple{(:a,),Tuple{Int}}}()
@test Pretty.table(nts) == """
| a   |
| --- |
empty row."""

Pretty.set(false)
@test sprint(show, nts) in ("NamedTuple{(:a,), Tuple{Int64}}[]", "NamedTuple{(:a,),Tuple{Int64}}[]")
Pretty.set(true)

@test Pretty._regularize_text("가1", 1) == "."
@test Pretty._regularize_text("가1", 2) == "가"
@test Pretty._regularize_text("가1", 3) == "가1"

nt = NamedTuple{(:id, :name, :salary),Tuple{Union{Missing, Int32},Union{Missing, String},Union{Missing, Float32}}}((1, "John", 10000.5f0))
@test Pretty.table(nt) == """
|   id | name   |    salary |
| ---- | ------ | --------- |
|    1 | John   |   10000.5 |"""

end # module test_octo_pretty


text:
# derivatives can either be applied in fourier space by multipliying by im*k or
# in map space by finite differences. this tracks which option to use for a
# given field
abstract type ∂modes end
struct fourier∂ <: ∂modes end
struct map∂ <: ∂modes end
promote_rule(::Type{map∂}, ::Type{fourier∂}) = fourier∂

# Flat{Nside,θpix,∂mode} is a flat sky pixelization with `Nside` pixels per side for (Ny, Nx)
# and pixels of width `θpix` arcmins, where derivatives are done according to ∂mode
abstract type Flat{Nside,θpix,∂mode<:∂modes,D} <: Pix end

# abstract type for the kind of data that can be stored in FlatFields
const AbstractRank2or3Array{T} = Union{AbstractArray{T,2},AbstractArray{T,3}}

# for convenience
Flat(;Nside, θpix=θpix₀, ∂mode=fourier∂, D=1) = Flat{Nside,θpix,∂mode,D}
Nside(::Type{P}) where {N,P<:Flat{N}} = N

# default angular resolution used by a number of convenience constructors
θpix₀ = 1

@doc """
The number of threads used by FFTW for CPU FFTs (default is the environment
variable `FFTW_NUM_THREADS`, or if that is not specified its
`Sys.CPU_THREADS÷2`). This must be set before creating any `FlatField` objects.
"""
FFTW_NUM_THREADS = nothing
@init global FFTW_NUM_THREADS = parse(Int,get(ENV,"FFTW_NUM_THREADS","$(Sys.CPU_THREADS÷2)"))

@doc """
Time-limit for FFT planning on CPU (default: 5 seconds). This must be set before
creating any `FlatField` objects.
"""
FFTW_TIMELIMIT = 5

function FlatInfo(T, Arr, θpix, Nside, D)

    FFTW.set_num_threads(FFTW_NUM_THREADS)

    Ny, Nx = Nside .* (1,1)
    Δx   = T(deg2rad(θpix/60))
    FFT  = plan_rfft(Arr{T}(undef,Ny,Nx,(D==1 ? () : (D,))...), (1,2); (Arr <: Array ? (timelimit=FFTW_TIMELIMIT,) : ())...)
    Δℓx  = T(2π/(Nx*Δx))
    Δℓy  = T(2π/(Ny*Δx))
    nyq  = T(2π/(2Δx))
    Ωpix = T(Δx^2)
    ky   = ifftshift(-Ny÷2:(Ny-1)÷2) .* Δℓy
    kx   = ifftshift(-Nx÷2:(Nx-1)÷2) .* Δℓx
    kmag = @. sqrt(kx'^2 + ky^2)
    ϕ    = @. angle(kx' + im*ky)[1:Ny÷2+1,:]
    sin2ϕ, cos2ϕ = @. sin(2ϕ), cos(2ϕ)
    if iseven(Ny)
        sin2ϕ[end, end:-1:(Nx÷2+2)] .= sin2ϕ[end, 2:Nx÷2]
    end

    @namedtuple(T, θpix, Nx, Ny, Nside, Δx, Δℓx, Δℓy, nyq, Ωpix, kx, ky, kmag, sin2ϕ=Arr(sin2ϕ), cos2ϕ=Arr(cos2ϕ), FFT)

end

function Cℓ_to_2D(::Type{P}, ::Type{T}, Cℓ) where {T,N,P<:Flat{N}}
    Complex{T}.(nan2zero.(Cℓ.(fieldinfo(P,T).kmag[1:fieldinfo(P,T).Ny÷2+1,:])))
end

@doc doc"""
    pixwin(θpix, ℓ)

Returns the pixel window function for square flat-sky pixels of width `θpix` (in
arcmin) evaluated at some `ℓ`s. This is the scaling of k-modes, the scaling of
the power spectrum will be pixwin^2.
"""
pixwin(θpix, ℓ) = @. sinc(ℓ*deg2rad(θpix/60)/2π)


text:
using Measures
using FactCheck

const imm = Length(:mm, 1)
const rmm = Length(:mm, 1//1)

const xx = Length(:xx, 1.0)

facts("Bounding boxes") do
    context("Valid bounding box constructors") do
        @fact BoundingBox((0mm, 1mm), 1mm, 2mm) -->
            BoundingBox((0mm, 1mm), 1mm, 2mm)
        @fact BoundingBox((2mm, 1mm, 3mm), 1mm, 2mm, 4mm) -->
            BoundingBox((2mm, 1mm, 3mm), 1mm, 2mm, 4mm)
    end

    context("Invalid bounding box constructors") do
        @fact_throws MethodError BoundingBox((0mm,), 1mm, 2mm)
        @fact_throws MethodError BoundingBox((0mm, 3mm, 8mm), 1mm, 4mm)
        @fact_throws MethodError BoundingBox((0mm, 3mm), 1mm, 2mm, 4mm)
        @fact_throws MethodError BoundingBox((0mm, 6mm, 3mm, 8mm), 1mm, 2mm, 4mm)
    end
end

facts("Measure constructors") do
    context("Length constructor") do
        @fact Length(:mm, 1) --> Length{:mm, Int}(1)
        @fact Length(:mm, 1.1) --> Length{:mm, Float64}(1.1)
        @fact 1mm |> typeof --> Length{:mm, Float64}
        @fact convert(Length{:mm, Float64}, 1imm) |> typeof --> Length{:mm, Float64}
    end

end

function check_binop(∘) #\circ
    @fact (1imm ∘ 2imm).value === 1 ∘ 2 --> true
    @fact (1mm ∘ 2imm).value === 1.0 ∘ 2 --> true
    @fact (2mm ∘ 1imm).value === 2.0 ∘ 1 --> true
    @fact typeof(2mm ∘ 3xx) <: Length --> false
    @fact typeof(2mm ∘ 3xx) <: Measures.BinaryOp --> true
end

function check_unaryop(>:)
    @fact >:1imm.value === >:1 --> true
    @fact >:1.0imm.value === >:1 --> false
    @fact >:1.0imm.value === >:1.0 --> true
end

facts("operators") do
    context("type promotion in +") do
        @fact typeof(1imm + 2imm) --> Length{:mm, Int}
        @fact typeof(1imm + 2mm) --> Length{:mm, Float64}
        @fact typeof(1mm + 3xx) --> Measures.Add{Length{:mm, Float64}, Length{:xx, Float64}}
    end

    context("equality") do
        @fact 1imm --> 1mm
        @fact 1imm --> not(1xx)
        @fact 1imm + 1xx --> not(1mm + 1xx) # perhaps make this true?
    end

    context("hash") do
        @fact hash(1imm) --> hash(1mm)
        @fact hash(1imm) --> not(hash(1xx))
        @fact isequal(1imm, 1mm) --> true
        @fact isequal(1imm, 1xx) --> false
    end

    context("+") do
        check_unaryop(+)
        check_binop(+)
    end
    context("-") do
        check_unaryop(-)
        check_binop(-)
    end
    context("max") do
        check_binop(max)
    end
    context("min") do
        check_binop(max)
    end
end

FactCheck.exitstatus()


text:
###############################################################################
#
# Parse string to unicode or latex strings
#
###############################################################################
"""
    parse_symbol(letter::String, unicode::Bool=false)

Return the unicode or latex string, given
- `id` Identifier for the greek letter, e.g., "alpha"
- `unicode` Optional. If true, return unicode; otherwise, return upgreek string
"""
function parse_symbol(id::String, unicode::Bool=false)
    # greek letter
    if id=="alpha" || id=="α"
        letter = (unicode ? "α" : "\$\\upalpha\$");
    elseif id=="Alpha" || id=="ALPHA" || id=="Α"
        letter = (unicode ? "α" : "\$\\upAlpha\$");
    elseif id=="beta" || id=="β"
        letter = (unicode ? "β" : "\$\\upbeta\$");
    elseif id=="Beta" || id=="BETA" || id=="Β"
        letter = (unicode ? "Β" : "\$\\upBeta\$");
    elseif id=="gamma" || id=="γ"
        letter = (unicode ? "γ" : "\$\\upgamma\$");
    elseif id=="Gamma" || id=="GAMMA" || id=="Γ"
        letter = (unicode ? "Γ" : "\$\\upGamma\$");
    elseif id=="delta" || id=="δ"
        letter = (unicode ? "δ" : "\$\\updelta\$");
    elseif id=="Delta" || id=="DELTA" || id=="Δ"
        letter = (unicode ? "Δ" : "\$\\upDelta\$");
    elseif id=="epsilon" || id=="ϵ"
        letter = (unicode ? "ϵ" : "\$\\upepsilon\$");
    elseif id=="Epsilon" || id=="EPSILON" || id=="Ε"
        letter = (unicode ? "Ε" : "\$\\upEpsilon\$");
    elseif id=="zeta" || id=="ζ"
        letter = (unicode ? "ζ" : "\$\\upzeta\$");
    elseif id=="Zeta" || id=="ZETA" || id=="Ζ"
        letter = (unicode ? "Ζ" : "\$\\upZeta\$");
    elseif id=="eta" || id=="η"
        letter = (unicode ? "η" : "\$\\upeta\$");
    elseif id=="Eta" || id=="ETA" || id=="Η"
        letter = (unicode ? "Η" : "\$\\upEta\$");
    elseif id=="theta" || id=="θ"
        letter = (unicode ? "θ" : "\$\\uptheta\$");
    elseif id=="Theta" || id=="THETA" || id=="Θ"
        letter = (unicode ? "Θ" : "\$\\upTheta\$");
    elseif id=="iota" || id=="ι"
        letter = (unicode ? "ι" : "\$\\upiota\$");
    elseif id=="Iota" || id=="IOTA" || id=="Ι"
        letter = (unicode ? "Ι" : "\$\\upIota\$");
    elseif id=="kappa" || id=="κ"
        letter = (unicode ? "κ" : "\$\\upkappa\$");
    elseif id=="Kappa" || id=="KAPPA" || id=="Κ"
        letter = (unicode ? "Κ" : "\$\\upKappa\$");
    elseif id=="lambda" || id=="λ"
        letter = (unicode ? "λ" : "\$\\uplambda\$");
    elseif id=="Lambda" || id=="LAMBDA" || id=="Λ"
        letter = (unicode ? "Λ" : "\$\\upLambda\$");
    elseif id=="mu" || id=="μ"
        letter = (unicode ? "μ" : "\$\\upmu\$");
    elseif id=="Mu" || id=="MU" || id=="M"
        letter = (unicode ? "M" : "\$\\upMu\$");
    elseif id=="nu" || id=="ν"
        letter = (unicode ? "ν" : "\$\\upnu\$");
    elseif id=="Nu" || id=="NU" || id=="N"
        letter = (unicode ? "N" : "\$\\upNu\$");
    elseif id=="xi" || id=="ξ"
        letter = (unicode ? "ξ" : "\$\\upxi\$");
    elseif id=="Xi" || id=="XI" || id=="Ξ"
        letter = (unicode ? "Ξ" : "\$\\upXi\$");
    elseif id=="pi" || id=="π"
        letter = (unicode ? "π" : "\$\\uppi\$");
    elseif id=="Pi" || id=="PI" || id=="Π"
        letter = (unicode ? "Π" : "\$\\upPi\$");
    elseif id=="rho" || id=="ρ"
        letter = (unicode ? "ρ" : "\$\\uprho\$");
    elseif id=="Rho" || id=="RHO" || id=="Ρ"
        letter = (unicode ? "Ρ" : "\$\\upRho\$");
    elseif id=="sigma" || id=="σ"
        letter = (unicode ? "σ" : "\$\\upsigma\$");
    elseif id=="Sigma" || id=="SIGMA" || id=="Σ"
        letter = (unicode ? "Σ" : "\$\\upSigma\$");
    elseif id=="tau" || id=="τ"
        letter = (unicode ? "τ" : "\$\\uptau\$");
    elseif id=="Tau" || id=="TAU" || id=="Τ"
        letter = (unicode ? "Τ" : "\$\\upTau\$");
    elseif id=="upsilon" || id=="υ"
        letter = (unicode ? "υ" : "\$\\upupsilon\$");
    elseif id=="Upsilon" || id=="UPSILON" || id=="Υ"
        letter = (unicode ? "Υ" : "\$\\upUpsilon\$");
    elseif id=="phi" || id=="ϕ"
        letter = (unicode ? "ϕ" : "\$\\upphi\$");
    elseif id=="Phi" || id=="PHI" || id=="Φ"
        letter = (unicode ? "Φ" : "\$\\upPhi\$");
    elseif id=="chi" || id=="χ"
        letter = (unicode ? "χ" : "\$\\upchi\$");
    elseif id=="Chi" || id=="CHI" || id=="Χ"
        letter = (unicode ? "Χ" : "\$\\upChi\$");
    elseif id=="psi" || id=="ψ"
        letter = (unicode ? "ψ" : "\$\\uppsi\$");
    elseif id=="Psi" || id=="PSI" || id=="Ψ"
        letter = (unicode ? "Ψ" : "\$\\upPsi\$");
    elseif id=="omega" || id=="ω"
        letter = (unicode ? "ω" : "\$\\upomega\$");
    elseif id=="Omega" || id=="OMEGA" || id=="Ω"
        letter = (unicode ? "Ω" : "\$\\upOmega\$");

    # unit symbol
    elseif lowercase(id) in ["deg", "degree", "°"]
        letter = (unicode ? "°" : "\$^{\\circ}\$");
    elseif lowercase(id) in ["degc", "degreeC", "°C"]
        letter = (unicode ? "°C" : "\$^{\\circ}\$C");

    # else, return id
    else
        letter = id;
    end

    return letter
end


text:
# DESCRIPTION
# Trim initialization.
#
# ------------------------------------------------------------------------------

# initialize state vectors and initial guess
# xf = [u v w p q r phi theta psi x y z]'
xf0=zeros(NFSTATES)
# xr = [MBC_flap, dMBC_flap, MBC_lag, dMBC_lag, MBC_inflow, azimuth, dyn twist]'
xr0=zeros(NRSTATES)
# xtr = TR_inflow
xtr0=zeros(NTRSTATES)
# initial guess for controls [%]
u0=[50.0; 50.0; 50.0; 50.0]
# set non-zero inflow states
xr0[17]=0.05
# initialize dynamic twist blade loads value [lb]
xr0[21]=WEIGHT/NB
# tail rotor inflow
xtr0=0.05
# initial guess for trim solution
xf0[1]=VXTRIM
xf0[2]=VYTRIM
xf0[3]=VZTRIM
xf0[9]=PSITRIM
xf0[10]=XNTRIM
xf0[11]=YETRIM
xf0[12]=-ALTTRIM
x0=[xf0; xr0; xtr0]
# set trim targets
xfdot0=zeros(NFSTATES)
xfdot0[9]=PSIDTRIM
xfdot0[10:12]=[VXTRIM*cos(PSITRIM)-VYTRIM*cos(PSITRIM); VXTRIM*sin(PSITRIM)+
               VYTRIM*cos(PSITRIM); VZTRIM]
xrdot0=zeros(NRSTATES)
xtrdot0=zeros(NTRSTATES)
xdot0=[xfdot0; xrdot0; xtrdot0]
xdot_targ=xdot0[TRIMTARG]
# Option to trim in zero sideselip or zero bank angle
if (VXTRIM>=60*1.688)
    if PSIDTRIM==0
        # Trim heading for zero bank angle above 60 kts
        x0[7]=0.0;
        TRIMVARS[7:8]=[8; 9]
    else
        x0[9]=atan2(VYTRIM,VXTRIM)
        TRIMVARS[7:8]=[7; 8]
        # initial guess for phi (using turn coord)
        x0[7]=atan(PSIDTRIM*sqrt(VXTRIM^2+VYTRIM^2+VZTRIM^2)/G)
        # initial guess for r
        x0[6]=G/sqrt(VXTRIM^2+VYTRIM^2+VZTRIM^2)*sin(x0[7])
        # initial guess for q
        x0[5]=x0[6]*tan(x0[7])
    end
else
    # otherwise set heading to flight path
    x0[9]=atan2(VYTRIM,VXTRIM)
    TRIMVARS[7:8]=[7; 8]
end
# time [s]
t=0.0

max_stars_repo_name: umbertosaetti/J-GenHel

max_issues_repo_name: umbertosaetti/J-GenHel

max_forks_repo_name: umbertosaetti/J-GenHel


text:
using ComplexDiff
using Test

@testset "ComplexDiff.jl" begin
    include("derivative.jl")
    include("jacobiantest.jl")
end


text:
# using Base.Test
# using MIToS.Information
# using MIToS.MSA

let false_clusters = SequenceClusters(zeros(20),zeros(20),rand(20)),
  seq = res"ARNDCQEGHILKMFPSTWYV-"

  print("""

Tests for ResidueContingencyTables
==================================
""")

  print("""

Tests for ResidueCount
----------------------
""")

  @test size(ResidueCount{Int, 1, false}().table) == (20,)
  @test size(ResidueCount{Int, 1, true}().table) == (21,)

  @test size(ResidueCount{Int, 2, false}().table) == (20,20)
  @test size(ResidueCount{Int, 2, true}().table) == (21,21)

  @test size(ResidueCount{Int, 3, false}().table) == (20,20,20)
  @test size(ResidueCount{Int, 3, true}().table) == (21,21,21)

  for ndim = 1:4
    @test size(ResidueCount{Int, ndim, true}().marginals) == (21, ndim)
    @test size(ResidueCount{Int, ndim, false}().marginals) == (20, ndim)
  end

  print("""

Iterate for N & UseGap:
""")
  for ndim = 1:4, usegap = Bool[true, false]
    println("# Iteration = N: $(ndim) UseGap: $(usegap)")
    nres = usegap ? 21 : 20

    println("Test nresidues for ResidueCount{Float64, $(ndim), $(usegap)}")
    N = zeros(ResidueCount{Float64, ndim, usegap})
    @test nresidues(N) == nres

    println("Test zeros for ResidueCount{Float64, $(ndim), $(usegap)}")
    @test sum(N.table) == zero(BigFloat)

    println("Test iteration interface for ResidueCount{Float64, $(ndim), $(usegap)}")
    @test collect(N) == zeros(Int, nres^ndim)

    println("Test size for ResidueCount{Float64, $(ndim), $(usegap)}")
    @test size(N) == size(N.table)

    println("Test indexing with i for ResidueCount{Float64, $(ndim), $(usegap)}")
    @test N[1] == zero(BigFloat)
    @test N[Int(Residue('R'))] == zero(BigFloat)

    println("Test setindex! with i for ResidueCount{Float64, $(ndim), $(usegap)}")
    N[3] = 1
    @test N[3] == one(BigFloat)

    println("Test update! for ResidueCount{Float64, $(ndim), $(usegap)}")
    fill!(N.table, 1)
    update!(N)
    @test N.total == length(N.table)
    @test N.marginals[1] == nres ^ ( ndim - 1 )

    println("Test apply_pseudocount! with Laplace Smoothing for ResidueCount{Float64, $(ndim), $(usegap)}")
    P = apply_pseudocount!(zeros(ResidueCount{Float64, ndim, usegap}), AdditiveSmoothing(1.0))
    @test P[1,1] == 1.0
    @test P.total == length(P.table)
    @test P.marginals[1] == nres ^ ( ndim - 1 )

    println("Test fill! with AdditiveSmoothing(0.05) for ResidueCount{Float64, $(ndim), $(usegap)}")
    fill!(P, AdditiveSmoothing(0.05))
    @test P[1,1] == 0.05
    @test P.total == length(P.table) * 0.05
    @test P.marginals[1] == 0.05 * ( nres ^ ( ndim - 1 ) )
  end
  print("""
Iteration end

""")

  print("""

Iterate for UseGap in [true, false]:
""")
  for usegap = Bool[true, false]
    println("# Iteration = UseGap: $(usegap)")
    nres = usegap ? 21 : 20

    N = zeros(ResidueCount{Float64, 2, usegap})

    println("Test indexing with ij for ResidueCount{Float64, 2, $(usegap)}")
    @test N[1,1] == zero(BigFloat)
    @test N[Int(Residue('R')), Int(Residue('R'))] == zero(BigFloat)

    println("Test setindex! with ij for ResidueCount{Float64, 2, $(usegap)}")
    N[3,3] = 1
    @test N[3,3] == one(BigFloat)

    println("Test count! for ResidueCount{Float64, 2, $(usegap)} and ResidueCount{Float64, 1, $(usegap)}")
    C = zeros(ResidueCount{Float64, 2, usegap})
    count!(C, seq, seq)
    @test C.table == eye(nres)
    @test (count!(zeros(ResidueCount{Float64, 1, usegap}), seq)).total == nres

    println("Test count! for ResidueCount{Float64, 5, $(usegap)} and  ResidueCount{Float64, 5, $(usegap)}}")
    @test (count!(zeros(ResidueCount{Float64, 5, usegap}), seq, seq, seq, seq, seq)).total == nres

    println("Test count! for ResidueCount{Float64, 4, $(usegap)} and  ResidueCount{Float64, 4, $(usegap)}}")
    C4 = count!(zeros(ResidueCount{Float64, 4, usegap}), seq, seq, seq, seq)
    @test C4[20,20,20,20] == 1.0
    @test C4[20,20,20,19] == 0.0
    # Also test sum
    @test sum(C4) == nres
  end
  print("""
Iteration end

""")

  print("""
Test similar
""")
  for ndim = 1:4
    @test size(similar(ResidueCount{Int, ndim, false}()).marginals) == (20, ndim)
    @test size(similar(ResidueCount{Int, ndim, true}()).marginals) == (21, ndim)
  end

  @test size(similar(ResidueCount{Int, 1, false}(),4).marginals) == (20, 4)
  @test isa(similar(ResidueCount{Int, 1, false}(),BigFloat).total, BigFloat)

  print("""
Test count! with SequenceClusters
""")
  @test_throws BoundsError count!(zeros(ResidueCount{Float64, 1, true}), false_clusters, seq)
  @test count!(zeros(ResidueCount{Float64, 1, false}), false_clusters, seq).table == getweight(false_clusters)
  @test count!(zeros(ResidueCount{Float64, 2, false}), false_clusters, seq, seq).marginals[:,1] == getweight(false_clusters)
  @test count!(zeros(ResidueCount{Float64, 3, false}), false_clusters, seq, seq, seq).marginals[:,1] == getweight(false_clusters)

  print("""
Test count
""")
  @test count(seq, weight=false_clusters).table == getweight(false_clusters)
  @test count(seq, seq, weight=false_clusters).marginals[:,1] == getweight(false_clusters)
  @test count(seq, seq, seq, weight=false_clusters).marginals[:,1] == getweight(false_clusters)

  @test count(seq, usegap=true).total == 21
  @test count(seq, seq, usegap=true).total == 21
  @test count(seq, seq, seq, usegap=true).total == 21

  @test count(AdditiveSmoothing(1.0), seq, usegap=true).total == 21 + 21
  @test count(AdditiveSmoothing(1.0), seq, seq, usegap=true).total == 21 + (21*21)

  @test count(AdditiveSmoothing(1), seq, usegap=true) == count(AdditiveSmoothing(1.0), seq, usegap=true)
  @test count(AdditiveSmoothing(1), seq, seq, usegap=true) == count(AdditiveSmoothing(1.0), seq, seq, usegap=true)

  print("""

Tests for ResidueProbability
----------------------
""")

  @test size(ResidueProbability{BigFloat, 1, false}().table) == (20,)
  @test size(ResidueProbability{BigFloat, 1, true}().table) == (21,)

  @test size(ResidueProbability{BigFloat, 2, false}().table) == (20,20)
  @test size(ResidueProbability{BigFloat, 2, true}().table) == (21,21)

  @test size(ResidueProbability{BigFloat, 3, false}().table) == (20,20,20)
  @test size(ResidueProbability{BigFloat, 3, true}().table) == (21,21,21)

  for ndim = 1:4
    @test size(ResidueProbability{BigFloat, ndim, true}().marginals) == (21, ndim)
    @test size(ResidueProbability{BigFloat, ndim, false}().marginals) == (20, ndim)
  end

  print("""

Iterate for N & UseGap:
""")
  for ndim = 1:4, usegap = Bool[true, false]
    println("# Iteration = N: $(ndim) UseGap: $(usegap)")
    nres = usegap ? 21 : 20

    println("Test zeros for ResidueProbability{BigFloat, $(ndim), $(usegap)}")
    N = zeros(ResidueProbability{BigFloat, ndim, usegap})
    @test sum(N.table) == zero(BigFloat)

    println("Test iteration interface for ResidueProbability{BigFloat, $(ndim), $(usegap)}")
    @test collect(N) == zeros(BigFloat, nres^ndim)

    println("Test size for ResidueProbability{BigFloat, $(ndim), $(usegap)}")
    @test size(N) == size(N.table)

    println("Test indexing with i for ResidueProbability{BigFloat, $(ndim), $(usegap)}")
    @test N[1] == zero(BigFloat)
    @test N[Int(Residue('R'))] == zero(BigFloat)

    println("Test setindex! with i for ResidueProbability{BigFloat, $(ndim), $(usegap)}")
    N[3] = 1.0
    @test N[3] == one(BigFloat)

    println("Test update! for ResidueProbability{BigFloat, $(ndim), $(usegap)}")
    N = zeros(ResidueProbability{BigFloat, ndim, usegap})
    N[1] = 0.5
    update!(N)
    @test N.marginals[1,1] == big"0.5"

    println("Test normalize! for ResidueProbability{BigFloat, $(ndim), $(usegap)}")
    normalize!(N)
    @test N[1] == one(BigFloat)
    @test N.marginals[1,1] == one(BigFloat)
    @test sum(N) == one(BigFloat)

    println("Test fill! with ResidueCount{Int, $(ndim), $(usegap)} for ResidueProbability{BigFloat, $(ndim), $(usegap)}")
    C = fill!(zeros(ResidueCount{Int, ndim, usegap}), AdditiveSmoothing(1))
    P = fill!(N, C)
    @test P[1] == one(BigFloat) / length(C.table)
    @test P.marginals[1,1] == one(BigFloat) / nres
    @test_approx_eq_eps sum(P) one(BigFloat) 1e-72
  end

  print("""
Iteration end

""")

  print("""
Test Float64
""")

  let N = zeros(ResidueProbability{Float64, 3, true}), C = fill!(zeros(ResidueCount{Int, 3, true}), AdditiveSmoothing(1))
    P = fill!(N, C)
    @test_approx_eq sum(P) one(Float64)
    @test isa(sum(P), Float64)
  end

  print("""

Iterate for UseGap in [true, false]:
""")
  for usegap = Bool[true, false]
    println("# Iteration = UseGap: $(usegap)")
    nres = usegap ? 21 : 20

    N = zeros(ResidueProbability{BigFloat, 2, usegap})

    println("Test indexing with ij for ResidueProbability{BigFloat, 2, $(usegap)}")
    @test N[1,1] == zero(BigFloat)
    @test N[Int(Residue('R')), Int(Residue('R'))] == zero(BigFloat)

    println("Test setindex! with ij for ResidueProbability{BigFloat, 2, $(usegap)}")
    N[3,3] = 1.0
    @test N[3,3] == one(BigFloat)
  end
  print("""
Iteration end

""")

  for ndim = 1:4
    @test size(similar(ResidueProbability{BigFloat, ndim, false}()).marginals) == (20, ndim)
    @test size(similar(ResidueProbability{BigFloat, ndim, true}()).marginals) == (21, ndim)
  end

  @test size(similar(ResidueProbability{BigFloat, 1, false}(),4).marginals) == (20, 4)

  print("""
Test pseudofrequencies
""")
  Pab = fill!(zeros(ResidueProbability{BigFloat, 2, false}), count(seq, seq))

  Gab = blosum_pseudofrequencies!(zeros(ResidueProbability{BigFloat, 2, false}), Pab)
  @test_approx_eq sum(Gab) 1.0

  @test_approx_eq sum( apply_pseudofrequencies!(copy(Pab), Gab, 1, 1) ) 1.0

  Pab_with_pseudofrequency = apply_pseudofrequencies!(copy(Pab), Gab, 10, 0)
  @test Pab_with_pseudofrequency == Pab

  print("""
Test probabilities
""")

  @test probabilities(Float64, seq, weight=false_clusters).table == getweight(false_clusters)/sum(getweight(false_clusters))
  @test probabilities(Float64, seq, seq, weight=false_clusters).marginals[:,1] == getweight(false_clusters)/sum(getweight(false_clusters))
  @test probabilities(Float64, seq, seq, seq, weight=false_clusters).marginals[:,1] == getweight(false_clusters)/sum(getweight(false_clusters))

  @test_approx_eq probabilities(Float64, seq)[1] 0.05
  @test_approx_eq probabilities(Float32, seq, seq)[2,2] Float32(0.05)
  @test_approx_eq probabilities(BigFloat, seq, seq, seq)[3,3,3] one(BigFloat)/20

  @test_approx_eq probabilities(Float64, AdditiveSmoothing(1.0), seq, usegap=true)[1] 1/21.0
  Pab = probabilities(Float64, AdditiveSmoothing(1.0), seq, seq, usegap=true)
  @test_approx_eq Pab[1,1] 2.0/(21 + 21*21)
  @test_approx_eq Pab[1,2] 1.0/(21 + 21*21)
  @test_approx_eq sum(Pab) 1.0

  @test_approx_eq probabilities(Float64, 1, 0, seq, seq)[2,2] 0.05

  print("""
Test delete_dimensions
""")
  Pxyz = probabilities(Float64, seq, seq, seq);
  @test delete_dimensions(Pxyz, 3) == probabilities(Float64, seq, seq)
  @test delete_dimensions(Pxyz, 3, 2) == probabilities(Float64, seq)
  Pxy = delete_dimensions(Pxyz, 3);
  @test delete_dimensions!(Pxy, Pxyz, 1) == probabilities(Float64, seq, seq)
  @test_approx_eq sum(Pxy) 1.0
  Nxyz = count(seq, seq, seq);
  Nxy = delete_dimensions(Nxyz, 3);
  @test Nxy == count(seq, seq)
  @test delete_dimensions(Nxyz, 3, 2) == count(seq)
  @test delete_dimensions!(Nxy, Nxyz, 1) == count(seq, seq)
  @test sum(Nxy) == Nxyz.total

end


text:
"""
Any StaticInjection struct that wants to support dynamic injectors must implement this
method to set the value.

The method is only for internal uses.
"""
function set_dynamic_injector!(
    static_injector::T,
    dynamic_injector::U,
) where {T <: Union{Nothing, StaticInjection}, U <: Union{Nothing, DynamicInjection}}
    error("set_dynamic_injector! is not implemented for $T")
end

"""
All DynamicInjection structs must implement this method to set the value.

The method is only for internal uses.
"""
function set_static_injector!(
    dynamic_injector::T,
    static_injector::U,
) where {T <: Union{Nothing, DynamicInjection}, U <: Union{Nothing, StaticInjection}}
    error("set_static_injector! is not implemented for $T")
end