text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: rbin
using AWS.Compat
using AWS.UUIDs

"""
    create_rule(resource_type, retention_period)
    create_rule(resource_type, retention_period, params::Dict{String,<:Any})

Creates a Recycle Bin retention rule. For more information, see  Create Recycle Bin
retention rules in the Amazon EC2 User Guide.

# Arguments
- `resource_type`: The resource type to be retained by the retention rule. Currently, only
  Amazon EBS snapshots are supported.
- `retention_period`: Information about the retention period for which the retention rule
  is to retain resources.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"Description"`: A brief description for the retention rule.
- `"ResourceTags"`: Information about the resource tags to use to identify resources that
  are to be retained by the retention rule. The retention rule retains only deleted snapshots
  that have one or more of the specified tag key and value pairs. If a snapshot is deleted,
  but it does not have any of the specified tag key and value pairs, it is immediately
  deleted without being retained by the retention rule. You can add the same tag key and
  value pair to a maximum or five retention rules.
- `"Tags"`: Information about the tags to assign to the retention rule.
"""
function create_rule(
    ResourceType, RetentionPeriod; aws_config::AbstractAWSConfig=global_aws_config()
)
    return rbin(
        "POST",
        "/rules",
        Dict{String,Any}(
            "ResourceType" => ResourceType, "RetentionPeriod" => RetentionPeriod
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function create_rule(
    ResourceType,
    RetentionPeriod,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return rbin(
        "POST",
        "/rules",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}(
                    "ResourceType" => ResourceType, "RetentionPeriod" => RetentionPeriod
                ),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    delete_rule(identifier)
    delete_rule(identifier, params::Dict{String,<:Any})

Deletes a Recycle Bin retention rule. For more information, see  Delete Recycle Bin
retention rules in the Amazon EC2 User Guide.

# Arguments
- `identifier`: The unique ID of the retention rule to delete.

"""
function delete_rule(identifier; aws_config::AbstractAWSConfig=global_aws_config())
    return rbin(
        "DELETE",
        "/rules/$(identifier)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function delete_rule(
    identifier,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return rbin(
        "DELETE",
        "/rules/$(identifier)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    get_rule(identifier)
    get_rule(identifier, params::Dict{String,<:Any})

Gets information about a Recycle Bin retention rule.

# Arguments
- `identifier`: The unique ID of the retention rule.

"""
function get_rule(identifier; aws_config::AbstractAWSConfig=global_aws_config())
    return rbin(
        "GET",
        "/rules/$(identifier)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function get_rule(
    identifier,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return rbin(
        "GET",
        "/rules/$(identifier)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    list_rules(resource_type)
    list_rules(resource_type, params::Dict{String,<:Any})

Lists the Recycle Bin retention rules in the Region.

# Arguments
- `resource_type`: The resource type retained by the retention rule. Only retention rules
  that retain the specified resource type are listed.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"MaxResults"`: The maximum number of results to return for the request in a single page.
  The remaining results can be seen by sending another request with the returned nextToken
  value. This value can be between 5 and 500. If maxResults is given a larger value than 500,
  you receive an error.
- `"NextToken"`: The token to use to retrieve the next page of results.
- `"ResourceTags"`: The tags used to identify resources that are to be retained by the
  retention rule.
"""
function list_rules(ResourceType; aws_config::AbstractAWSConfig=global_aws_config())
    return rbin(
        "POST",
        "/list-rules",
        Dict{String,Any}("ResourceType" => ResourceType);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function list_rules(
    ResourceType,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return rbin(
        "POST",
        "/list-rules",
        Dict{String,Any}(
            mergewith(_merge, Dict{String,Any}("ResourceType" => ResourceType), params)
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    list_tags_for_resource(resource_arn)
    list_tags_for_resource(resource_arn, params::Dict{String,<:Any})

Lists the tags assigned a specific resource.

# Arguments
- `resource_arn`: The Amazon Resource Name (ARN) of the resource for which to list the tags.

"""
function list_tags_for_resource(
    resourceArn; aws_config::AbstractAWSConfig=global_aws_config()
)
    return rbin(
        "GET",
        "/tags/$(resourceArn)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function list_tags_for_resource(
    resourceArn,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return rbin(
        "GET",
        "/tags/$(resourceArn)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    tag_resource(tags, resource_arn)
    tag_resource(tags, resource_arn, params::Dict{String,<:Any})

Assigns tags to the specified resource.

# Arguments
- `tags`: Information about the tags to assign to the resource.
- `resource_arn`: The Amazon Resource Name (ARN) of the resource to which to assign the
  tags.

"""
function tag_resource(Tags, resourceArn; aws_config::AbstractAWSConfig=global_aws_config())
    return rbin(
        "POST",
        "/tags/$(resourceArn)",
        Dict{String,Any}("Tags" => Tags);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function tag_resource(
    Tags,
    resourceArn,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return rbin(
        "POST",
        "/tags/$(resourceArn)",
        Dict{String,Any}(mergewith(_merge, Dict{String,Any}("Tags" => Tags), params));
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    untag_resource(resource_arn, tag_keys)
    untag_resource(resource_arn, tag_keys, params::Dict{String,<:Any})

Unassigns a tag from a resource.

# Arguments
- `resource_arn`: The Amazon Resource Name (ARN) of the resource from which to unassign the
  tags.
- `tag_keys`: Information about the tags to unassign from the resource.

"""
function untag_resource(
    resourceArn, tagKeys; aws_config::AbstractAWSConfig=global_aws_config()
)
    return rbin(
        "DELETE",
        "/tags/$(resourceArn)",
        Dict{String,Any}("tagKeys" => tagKeys);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function untag_resource(
    resourceArn,
    tagKeys,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return rbin(
        "DELETE",
        "/tags/$(resourceArn)",
        Dict{String,Any}(mergewith(_merge, Dict{String,Any}("tagKeys" => tagKeys), params));
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    update_rule(identifier)
    update_rule(identifier, params::Dict{String,<:Any})

Updates an existing Recycle Bin retention rule. For more information, see  Update Recycle
Bin retention rules in the Amazon EC2 User Guide.

# Arguments
- `identifier`: The unique ID of the retention rule to update.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"Description"`: The retention rule description.
- `"ResourceTags"`: Information about the resource tags to use to identify resources that
  are to be retained by the retention rule. The retention rule retains only deleted snapshots
  that have one or more of the specified tag key and value pairs. If a snapshot is deleted,
  but it does not have any of the specified tag key and value pairs, it is immediately
  deleted without being retained by the retention rule.  You can add the same tag key and
  value pair to a maximum or five retention rules.
- `"ResourceType"`: The resource type to be retained by the retention rule. Currently, only
  Amazon EBS snapshots are supported.
- `"RetentionPeriod"`: Information about the retention period for which the retention rule
  is to retain resources.
"""
function update_rule(identifier; aws_config::AbstractAWSConfig=global_aws_config())
    return rbin(
        "PATCH",
        "/rules/$(identifier)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function update_rule(
    identifier,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return rbin(
        "PATCH",
        "/rules/$(identifier)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end


text:
using Test
using DynamicalBilliards
using DynamicalBilliards.Testing

function test_lyapunov_spectrum(p, bd, t = 1e6, error_level = 1e-2)
    # calculate spectrum
    exps = lyapunovspectrum!(p, bd, t)
    sumpair = exps[1] + exps[4]

    # these properties should be true for all billiards
    @testset "λ₁ + λ₄ ≈ 0" begin
        @test abs(sumpair) < error_level
    end

    @testset "λ₂ ≈ 0, λ₃ ≈ 0" begin
        @test abs(exps[2]) < error_level
        @test abs(exps[3]) < error_level
    end

    exps = lyapunovspectrum!(p, bd, 10000)
    sumpair = exps[1] + exps[4]

    # these properties should be true for all billiards
    @testset "λ₁ + λ₄ ≈ 0" begin
        @test abs(sumpair) < error_level
    end

    @testset "λ₂ ≈ 0, λ₃ ≈ 0" begin
        @test abs(exps[2]) < error_level
        @test abs(exps[3]) < error_level
    end
end

billiards_testset("Properties of Lyapunov spectrum",
                  test_lyapunov_spectrum; caller = simple_tests)

function test_lyapunov_magnetic_limit(args...)
    bd = billiard_bunimovich(0.1)

    t = 1e5
    error_level = 5e-2

    for j ∈ 1:10
        pmag = randominside(bd, 1e-3)
        plin = Particle(pmag.pos, pmag.vel, SVector{2,Float64}(0,0))

        exp_mag = lyapunovspectrum!(pmag, bd, t)
        exp_lin = lyapunovspectrum!(plin, bd, t)

        @test exp_mag[1] - exp_lin[1] < error_level
    end
end

billiards_testset("Lyapunov spectrum for small B",
                  identity; caller=test_lyapunov_magnetic_limit)

function test_lyapunov_values(args...)
    t = 20000.0
    radius = 1.0

    spaces = [2.0, 2.5, 3.0, 3.5, 4.0]
    # based on Gaspard et al (see docs) & DynamicalBilliards v2.5
    expected_values = [3.6, 1.4, 0.8, 0.6, 0.5]
    error_level = 0.2

    for (i, space) in enumerate(spaces)
        bd = billiard_polygon(6, space/(sqrt(3)); setting = "periodic")
        disc = Disk([0., 0.], radius)
        billiard = Billiard(bd.obstacles..., disc)
        p = randominside(billiard)
        λ = lyapunovspectrum(p, billiard, t)[1]
        @test abs(λ - expected_values[i]) < error_level
    end
end

billiards_testset("Lyapunov numerical values",
                  identity, caller=test_lyapunov_values)

using LinearAlgebra
function test_perturbationgrowth(p, bd)

    tmax = 100.0
    error_level = 1e-1

    t, R, o = perturbationgrowth(p, bd, tmax)
    λ = lyapunovspectrum(p, bd, tmax)
    Δ = perturbationevolution(R)

    norms = log.(norm.(Δ))
    actual = norms[end]; i = length(norms) - 1
    while isinf(actual)
        actual = norms[i]
        i -= 1
    end

    λ_estimate = actual/t[i]
    @test abs(λ[1] - λ_estimate) < error_level

    nmax = 500
    error_level = 1e-1

    t, R, o = perturbationgrowth(p, bd, nmax)
    λ = lyapunovspectrum(p, bd, nmax)
    Δ = perturbationevolution(R)

    @test length(t) == 2*nmax

    norms = log.(norm.(Δ))
    actual = norms[end]; i = length(norms) - 1
    while isinf(actual)
        actual = norms[i]
        i -= 1
    end

    λ_estimate = actual/t[i]
    @test abs(λ[1] - λ_estimate) < error_level
end

billiards_testset("Compare Lyapunovs and perturbation growth",
                  test_perturbationgrowth; caller=ergodic_tests)

ext: jl


text:
# Adjoints needed in polar formulation
#

"""
    AdjointStackObjective

An object for storing the adjoint stack in the adjoint objective computation

"""
struct AdjointStackObjective{VT}
    ∇fₓ::VT
    ∇fᵤ::VT
    ∂pg::VT
    ∂vm::VT
    ∂va::VT
    jvₓ::VT
    jvᵤ::VT
end

function put_active_power_injection!(fr, v_m, v_a, adj_v_m, adj_v_a, adj_P, ybus_re, ybus_im)
    @inbounds for c in ybus_re.colptr[fr]:ybus_re.colptr[fr+1]-1
        to = ybus_re.rowval[c]
        aij = v_a[fr] - v_a[to]
        cθ = ybus_re.nzval[c]*cos(aij)
        sθ = ybus_im.nzval[c]*sin(aij)
        adj_v_m[fr] += v_m[to] * (cθ + sθ) * adj_P
        adj_v_m[to] += v_m[fr] * (cθ + sθ) * adj_P

        adj_aij = -(v_m[fr]*v_m[to]*(ybus_re.nzval[c]*sin(aij)))
        adj_aij += v_m[fr]*v_m[to]*(ybus_im.nzval[c]*cos(aij))
        adj_aij *= adj_P
        adj_v_a[to] += -adj_aij
        adj_v_a[fr] += adj_aij
    end
end

function put_reactive_power_injection!(fr, v_m, v_a, adj_v_m, adj_v_a, adj_P, ybus_re, ybus_im)
    @inbounds for c in ybus_re.colptr[fr]:ybus_re.colptr[fr+1]-1
        to = ybus_re.rowval[c]
        aij = v_a[fr] - v_a[to]
        cθ = ybus_im.nzval[c]*cos(aij)
        sθ = ybus_re.nzval[c]*sin(aij)
        adj_v_m[fr] += v_m[to] * (sθ - cθ) * adj_P
        adj_v_m[to] += v_m[fr] * (sθ - cθ) * adj_P

        adj_aij = v_m[fr]*v_m[to]*(ybus_re.nzval[c]*cos(aij))
        adj_aij += v_m[fr]*v_m[to]*(ybus_im.nzval[c]*(sin(aij)))
        adj_aij *= adj_P
        adj_v_a[to] += -adj_aij
        adj_v_a[fr] += adj_aij
    end
end

function put!(
    polar::PolarForm{T, IT, VT, AT},
    ::State,
    ∂x::VT,
    ∂vmag::VT,
    ∂vang::VT,
) where {T, IT, VT, AT}
    npv = PS.get(polar.network, PS.NumberOfPVBuses())
    npq = PS.get(polar.network, PS.NumberOfPQBuses())
    # build vector x
    for (i, p) in enumerate(polar.network.pv)
        ∂x[i] = ∂vang[p]
    end
    for (i, p) in enumerate(polar.network.pq)
        ∂x[npv+i] = ∂vang[p]
        ∂x[npv+npq+i] = ∂vmag[p]
    end
end

function put!(
    polar::PolarForm{T, IT, VT, AT},
    ::Control,
    ∂u::VT,
    ∂vmag::VT,
    ∂pbus::VT,
) where {T, IT, VT, AT}
    npv = PS.get(polar.network, PS.NumberOfPVBuses())
    nref = PS.get(polar.network, PS.NumberOfSlackBuses())
    # build vector u
    for (i, p) in enumerate(polar.network.ref)
        ∂u[i] = ∂vmag[p]
    end
    for (i, p) in enumerate(polar.network.pv)
        ∂u[nref + i] = 0.0
        ∂u[nref + npv + i] = ∂vmag[p]
    end
end

@kernel function put_adjoint_kernel!(
    adj_u, adj_x, adj_vmag, adj_vang, adj_pg,
    index_pv, index_pq, index_ref, pv_to_gen,
)
    i = @index(Global, Linear)
    npv = length(index_pv)
    npq = length(index_pq)
    nref = length(index_ref)

    # PQ buses
    if i <= npq
        bus = index_pq[i]
        adj_x[npv+i] =  adj_vang[bus]
        adj_x[npv+npq+i] = adj_vmag[bus]
    # PV buses
    elseif i <= npq + npv
        i_ = i - npq
        bus = index_pv[i_]
        i_gen = pv_to_gen[i_]
        adj_u[nref + npv + i_] = adj_vmag[bus]
        adj_u[nref + i_] += adj_pg[i_gen]
        adj_x[i_] = adj_vang[bus]
    # SLACK buses
    elseif i <= npq + npv + nref
        i_ = i - npq - npv
        bus = index_ref[i_]
        adj_u[i_] = adj_vmag[bus]
    end
end

function put(
    polar::PolarForm{T, VT, AT},
    ::PS.Generator,
    ::PS.ActivePower,
    obj_autodiff::AdjointStackObjective,
    buffer::PolarNetworkState
) where {T, VT, AT}
    ngen = PS.get(polar.network, PS.NumberOfGenerators())
    nbus = PS.get(polar.network, PS.NumberOfBuses())
    npv = PS.get(polar.network, PS.NumberOfPVBuses())
    npq = PS.get(polar.network, PS.NumberOfPQBuses())
    nref = PS.get(polar.network, PS.NumberOfSlackBuses())

    index_pv = polar.indexing.index_pv
    index_ref = polar.indexing.index_ref
    index_pq = polar.indexing.index_pq
    pv_to_gen = polar.indexing.index_pv_to_gen
    ref_to_gen = polar.indexing.index_ref_to_gen

    # Get voltages. This is only needed to get the size of adjvmag and adjvang
    vmag = buffer.vmag
    vang = buffer.vang

    adj_pg = obj_autodiff.∂pg
    adj_x = obj_autodiff.∇fₓ
    adj_u = obj_autodiff.∇fᵤ
    adj_vmag = obj_autodiff.∂vm
    adj_vang = obj_autodiff.∂va
    fill!(adj_vmag, 0.0)
    fill!(adj_vang, 0.0)
    fill!(adj_x, 0.0)
    fill!(adj_u, 0.0)

    for i in 1:nref
        bus = index_ref[i]
        i_gen = ref_to_gen[i]
        # pg[i] = inj + polar.active_load[bus]
        adj_inj = adj_pg[i_gen]
        put_active_power_injection!(bus, vmag, vang, adj_vmag, adj_vang, adj_inj, polar.ybus_re, polar.ybus_im)
    end
    if isa(adj_x, Array)
        kernel! = put_adjoint_kernel!(CPU(), 1)
    else
        kernel! = put_adjoint_kernel!(CUDADevice(), 256)
    end
    ev = kernel!(adj_u, adj_x, adj_vmag, adj_vang, adj_pg,
                 index_pv, index_pq, index_ref, pv_to_gen,
                 ndrange=nbus)
    wait(ev)

    return
end

function ∂cost(polar::PolarForm, ∂obj::AdjointStackObjective, buffer::PolarNetworkState)
    pg = buffer.pg
    coefs = polar.costs_coefficients
    c3 = @view coefs[:, 3]
    c4 = @view coefs[:, 4]
    # Return adjoint of quadratic cost
    ∂obj.∂pg .= c3 .+ 2.0 .* c4 .* pg
    put(polar, PS.Generator(), PS.ActivePower(), ∂obj, buffer)
end


text:
module EchelleTools
using Reexport

include("constants.jl")

include("abstract_types.jl")

include("instruments.jl")

include("util.jl")
#=
export default_telluric_template, default_stellar_template
export make_model_wavelength_range
export ManifestFormatEntry
=#

include("psf.jl")

include("blaze/blaze.jl")

include("obs_data.jl")
export EchelleObservation, EchelleObservationSet
export make_average_observation, z_bc

include("model.jl")
export SpectralOrderModel, SpectralModel
export get_order

export ObservationModel, ObsSeriesModel
export eval_model, eval_model_baseline

#=
include("analyze.jl")
export mean_spectrum, median_spectrum
export lambda_ssb_from_obs, log_lambda_ssb_from_obs
export make_blaze_and_template_1d, make_blaze_and_template_2d
export make_stellar_basis, default_telluric_model
export make_template_1d
export make_spectral_order_model, make_spectral_model
export make_observations_model
# export model_spectrum_1d, model_spectrum_2d

include("pca.jl")
export interp_one_obs_order_for_pca_1d, fit_pca_obs_order_interpolated
# =#

#=
include("toy_data.jl")
export spectral_line, make_toy_data
=#

# Instrument specific files
include("HARPSN/HARPSN.jl")
#using HARPSN
#=
read_harpsn_manifest_file = HARPSN.read_harpsn_manifest_file
read_harpsn_manifest_files = HARPSN.read_harpsn_manifest_files
get_harpsn_data_filenames = HARPSN.get_harpsn_data_filenames
get_harpsn_data_filenames = HARPSN.get_harpsn_data_filenames
read_harpsn_data = HARPSN.read_harpsn_data
read_harpsn_blaze_list = HARPSN.read_blaze_list
export read_harpsn_manifest_file, read_harpsn_manifest_files, get_harpsn_data_filenames, read_harpsn_data
export get_harpsn_data_filenames_and_bsrv
export read_harpsn_blaze_list
=#

#=
include("HPF/HPF.jl")
read_hpf_manifest_file = HPF.read_hpf_manifest_file
read_hpf_manifest_files = HPF.read_hpf_manifest_files
get_hpf_data_filenames = HPF.get_hpf_data_filenames
read_hpf_data = HPF.read_hpf_data
export read_hpf_manifest_file, read_hpf_manifest_files, get_hpf_data_filenames, read_hpf_data
=#

end # module EchelleTools


text:
# mercer mapping

export MercerMap, apply, linear, rbf

immutable MercerMap
    X::Matrix           # original reference data
    kernel::Symbol      # kernel function
    K::Matrix           # (centered) kernel matrix of reference data
    d::Vector           # largest eigenvalues of kernel matrix
    V::Matrix           # eigenvectors of kernel matrix
end

# constructing mercer map
function MercerMap(X::Matrix, kernel::Symbol)
    kernelfcn = eval(kernel)
    K = center(kernelfcn(X, X))

    d, V = eig(Symmetric(K))

    # consider dimensions with eigenvalues > 1e-9
    i = d .> 1e-9
    V = V[:,i]
    d = sqrt(d[i])

    MercerMap(X, kernel, K, d, V)
end

# map reference data
apply(map::MercerMap) = map.V*diagm(map.d)
function apply(map::MercerMap, X::Matrix)
    kernelfcn = eval(map.kernel)
    KT = center(map.K, kernelfcn(map.X, X))
    KT'*map.V*diagm(1./map.d)
end

# center data to origin in feature space (M. Meila: Data Centering in Feature Space, Eq. 17)
function center(K::Matrix, KT::Matrix)
    m = size(K,1)
    M1 = 1/m*ones(size(K))
    N1 = ones(size(KT)) / m
    KT = KT - K*N1 - M1*KT + M1*K*N1
end
center(K::Matrix) = center(K, K)

# kernel functions
linear(X::Matrix, Y::Matrix) = full(X*Y')
rbf(X::Matrix, Y::Matrix; sigma::Float64=1.0) = full(exp(-dist2(X, Y)/(2*sigma^2)))

function dist2(X::Matrix, Y::Matrix)
    xx = sum(X.*X, 2)
    yy = sum(Y.*Y, 2)
    xy = X*Y'
    abs(repmat(xx, 1, size(yy, 1)) + repmat(yy', size(xx, 1), 1) - 2*xy)
end


text:
function grade(mid, fin, re)
    score = mid+fin
    (mid == -1 || fin == -1 || score < 30) && return "F"
    score >= 80 && return "A"
    65 <= score < 80 && return "B"
    50 <= score < 65 && return "C"
    if 30 <= score < 50
        if re >= 50
            return "C"
        else
            return "D"
        end
    end
end

while true
    mid, fin, re = parse.(Int, split(readline()))
    mid == fin == re == -1 && break
    println(grade(mid, fin, re))
end


text:
using StatsBase, StatsFuns, JLD,  Hwloc

counts = Hwloc.histmap(Hwloc.topology_load())
ncores = counts[:Core]
if nworkers()<ncores
    addprocs(ncores   - nprocs())
end

using DynNets
@everywhere using DynNets
@everywhere using StatsFuns
N_est = 500
Tvals = [50 300  ]#500#
N_ind = length(Tvals)[1]

N = 10 ;
GBAest =1; GWest = 10; useStartVal = false

groupsIndsEst = [[DynNets.distributeAinVecN(Array{Int}(1:GWest),N)];
                   [DynNets.distributeAinVecN(Array{Int}(1:GBAest),GWest)]]

est_times = Array{Float64}(zeros(N_est))
est_conv_flag = falses(N_est)
est_flag = falses(N_est)

sizesPar = [N,GBAest,GBAest]
for i =1:3
    ex1 = "estPar_$(i)  = SharedArray{Float64}(zeros(N_ind,N_est, $(sizesPar[i]) ))"
    eval(parse(ex1))
end

ObsAllT = Array{Array{Int16,3},1}(N_ind)
ParTvDgpAllT = Array{Array{Float64,3},1}(N_ind)

#define save path and list variables to save
save_fold = "./data/estimatesTest/"
    file_name = "EstTest_$(N)_$(Tvals)_$(N_est)_$(useStartVal)_Sin_Miss_$(GWest)_$(GBAest).jld"
    save_path = save_fold*file_name#

est_times = SharedArray{Float64}(zeros(N_ind,N_est))
est_conv_flag = SharedArray{Bool,2}(N_ind,N_est)
est_flag = SharedArray{Bool,2}(N_ind,N_est)
#
@progress "T"  for ind = 1:N_ind
#   ind=1

    T = Tvals[ind]
    ObsTmp = SharedArray{Int16,3}(Tvals[ind],N,N_est)
    ParTvDgpTmp = SharedArray{Float64,3}(Tvals[ind],N,N_est)

    @sync @parallel     for n = 1:N_est
#     n=1

        Y_T,fit_T  = DynNets.sampleDynNetsBin1(T,"FIT-MULTI-SIN";N=N)
        sampleDegsT = squeeze(sum(Y_T,2),2)

        ParTvDgpTmp[:,:,n] = fit_T
        ObsTmp[:,:,n] = sampleDegsT

        Model2Est =  DynNets.SdErgm1_UmBAgr(sampleDegsT,[zeros(GWest),zeros(GBAest),zeros(GBAest)],
                                                    groupsIndsEst)
        #estimate
        useStartVal ?   start_val = [sim_W, sim_B, sim_A] : start_val = zeros(3,3)
        est_times[ind,n] = @elapsed  hat_pars , conv_flag =
                            DynNets.estimate(Model2Est,start_values = start_val )

        estPar_1[ind,n,:] = hat_pars[1]
        estPar_2[ind,n,:] = hat_pars[2]
        estPar_3[ind,n,:] = hat_pars[3]
        est_conv_flag[ind,n] = conv_flag
        est_flag[ind,n] = true

   end

   ObsAllT[ind] = ObsTmp
   ParTvDgpAllT[ind] = ParTvDgpTmp

    # save

    @save(save_path,groupsIndsEst,ObsAllT,ParTvDgpAllT ,estPar_1,estPar_2,estPar_3,est_conv_flag,est_times)
    println([N,GWest,T,sum(est_times[ind,:]),sum(est_conv_flag[ind,:])/sum(est_flag[ind,:])])
end

##

 ##


text:
######### CmdStan program example  ###########

using CmdStan

ProjDir = dirname(@__FILE__)
cd(ProjDir)

bernoullimodel = "
data {
  int<lower=1> N;
  int<lower=0,upper=1> y[N];
}
parameters {
  real<lower=0,upper=1> theta;
}
model {
  theta ~ beta(1,1);
  y ~ bernoulli(theta);
}
";

tmpdir = ProjDir*"/tmp"

stanmodel = Stanmodel(Sample(save_warmup=true, num_warmup=1000,
  num_samples=2000, thin=1), name="bernoulli", model=bernoullimodel,
  printsummary=false, tmpdir=tmpdir);

println("\nTest for Atom completed, stanmodel should not have been displayed\n")


text:
function subleading(a, n)
 d = map(D -> 1:D, size(a)[1:end-1])
 sub(a, d..., n)

end

function subleading(a::Matrix, n)
 sub(a, 1:size(a,1), n)
end

function leading(a, n)
 d = map(D -> 1:D, size(a)[1:end-1])
 a[d..., n]

end

function leading(a::Matrix, n)
 a[1:size(a,1), n]
end


text:
"""
    moduleof

Returns the module associated with a method.
Note: this is not the module where it was defined (m.module)
but rather the module for which it's function was defined.
"""
moduleof(m::Method) = moduleof(functiontypeof(m))
moduleof(sig::DataType) = sig.name.module
moduleof(sig::UnionAll) = moduleof(sig.body)

functiontypeof(m::Method) = parameter_typeof(m.sig)[1]
parameter_typeof(sig::UnionAll) = parameter_typeof(sig.body)
parameter_typeof(sig::DataType) = sig.parameters


text:
using MultivariateFunctions
using Dates

tol = 10*eps()
const global_base_date = Date(2000,1,1)
StartDate = Date(2018, 7, 21)

x = Array{Date}(undef, 100)
for i in 1:100
    x[i] = StartDate +Dates.Day(2* (i-1))
end

function ff(x::Date)
    days_between = years_from_global_base(x)
    return log(days_between) + sqrt(days_between)
end
y = ff.(x)

splinex = create_quadratic_spline(x,y; dim_name = :x)
collect(keys(splinex.thresholds_))[1] == :x

spline = create_quadratic_spline(x,y)
# Test if interpolating
all(abs.(evaluate.(Ref(spline), x) .- y) .< tol)

# Testing third derivatives
third_derivatives = evaluate.(Ref(derivative(derivative(derivative(spline)))), x)
maximum(third_derivatives) < tol

# Testing Integrals
function analytic_integral(lhs,rhs)
    lhs_in_days = years_from_global_base(lhs)
    rhs_in_days = years_from_global_base(rhs)
    return rhs_in_days*log(rhs_in_days) - rhs_in_days + (2/3)*rhs_in_days^(3/2) - (lhs_in_days*log(lhs_in_days) - lhs_in_days)- (2/3)*lhs_in_days^(3/2)
end

lhs = StartDate
rhs = StartDate + Dates.Month(16)
numerical_integral  = integral(spline, lhs,rhs)

analytical = analytic_integral(lhs,rhs)
abs(  analytical - numerical_integral  ) < 0.001

#= Speed test of MultivariateFunctions vs Univariate vs SchumakerSpline
using SchumakerSpline
using UnivariateFunctions
x = collect(range(0.001, stop=10, length=1000))
y = log.(x)
@time schum = SchumakerSpline.Schumaker(x, y)
@time univ = UnivariateFunctions.create_quadratic_spline(x,y)
@time multiv = MultivariateFunctions.create_quadratic_spline(x,y)
@time SchumakerSpline.evaluate(schum, 7.8)
@time UnivariateFunctions.evaluate(univ, 7.8)
@time MultivariateFunctions.evaluate(multiv, 7.8)
@time MultivariateFunctions.evaluate(multiv, 7.8; variable = :default)
@time MultivariateFunctions.evaluate(multiv, Dict{Symbol,Float64}(:default => 7.8))
=#


text:
struct Assembler{T}
    I::Vector{Int}
    J::Vector{Int}
    V::Vector{T}
end

function Assembler(N)
    I = Int[]
    J = Int[]
    V = Float64[]
    sizehint!(I, N)
    sizehint!(J, N)
    sizehint!(V, N)

    Assembler(I, J, V)
end

"""
    start_assemble([N=0]) -> Assembler

Call before starting an assembly.

Returns an `Assembler` type that is used to hold the intermediate
data before an assembly is finished.
"""
function start_assemble(N::Int=0)
    return Assembler(N)
end

"""
    assemble!(a, Ke, edof)

Assembles the element matrix `Ke` into `a`.
"""
function assemble!(a::Assembler{T}, edof::AbstractVector{Int}, Ke::AbstractMatrix{T}) where {T}
    n_dofs = length(edof)
    append!(a.V, Ke)
    @inbounds for j in 1:n_dofs
        append!(a.I, edof)
        for i in 1:n_dofs
            push!(a.J, edof[j])
        end
    end
end

"""
    end_assemble(a::Assembler) -> K

Finalizes an assembly. Returns a sparse matrix with the
assembled values.
"""
function end_assemble(a::Assembler)
    return sparse(a.I, a.J, a.V)
end

"""
    assemble!(g, ge, edof)

Assembles the element residual `ge` into the global residual vector `g`.
"""
@propagate_inbounds function assemble!(g::AbstractVector{T}, edof::AbstractVector{Int}, ge::AbstractVector{T}) where {T}
    @boundscheck checkbounds(g, edof)
    @inbounds for i in 1:length(edof)
        g[edof[i]] += ge[i]
    end
end

abstract type AbstractSparseAssembler end

struct AssemblerSparsityPattern{Tv,Ti} <: AbstractSparseAssembler
    K::SparseMatrixCSC{Tv,Ti}
    f::Vector{Tv}
    permutation::Vector{Int}
    sorteddofs::Vector{Int}
end
struct AssemblerSymmetricSparsityPattern{Tv,Ti} <: AbstractSparseAssembler
    K::Symmetric{Tv,SparseMatrixCSC{Tv,Ti}}
    f::Vector{Tv}
    permutation::Vector{Int}
    sorteddofs::Vector{Int}
end

@inline getsparsemat(a::AssemblerSparsityPattern) = a.K
@inline getsparsemat(a::AssemblerSymmetricSparsityPattern) = a.K.data

start_assemble(f::Vector, K::Union{SparseMatrixCSC, Symmetric}) = start_assemble(K, f)
function start_assemble(K::SparseMatrixCSC, f::Vector=Float64[])
    fill!(K.nzval, 0.0)
    fill!(f, 0.0)
    AssemblerSparsityPattern(K, f, Int[], Int[])
end
function start_assemble(K::Symmetric, f::Vector=Float64[])
    fill!(K.data.nzval, 0.0)
    fill!(f, 0.0)
    AssemblerSymmetricSparsityPattern(K, f, Int[], Int[])
end

@propagate_inbounds function assemble!(A::AbstractSparseAssembler, dofs::AbstractVector{Int}, Ke::AbstractMatrix)
    assemble!(A, dofs, Ke, eltype(Ke)[])
end
@propagate_inbounds function assemble!(A::AbstractSparseAssembler, dofs::AbstractVector{Int}, fe::AbstractVector, Ke::AbstractMatrix)
    assemble!(A, dofs, Ke, fe)
end
@propagate_inbounds function assemble!(A::AssemblerSparsityPattern, dofs::AbstractVector{Int}, Ke::AbstractMatrix, fe::AbstractVector)
    _assemble!(A, dofs, Ke, fe, false)
end
@propagate_inbounds function assemble!(A::AssemblerSymmetricSparsityPattern, dofs::AbstractVector{Int}, Ke::AbstractMatrix, fe::AbstractVector)
    _assemble!(A, dofs, Ke, fe, true)
end

@propagate_inbounds function _assemble!(A::AbstractSparseAssembler, dofs::AbstractVector{Int}, Ke::AbstractMatrix, fe::AbstractVector, sym::Bool)
    if length(fe) != 0
        assemble!(A.f, dofs, fe)
    end

    K = getsparsemat(A)
    permutation = A.permutation
    sorteddofs = A.sorteddofs
    @boundscheck checkbounds(K, dofs, dofs)
    resize!(permutation, length(dofs))
    resize!(sorteddofs, length(dofs))
    copy!(sorteddofs, dofs)
    sortperm2!(sorteddofs, permutation)

    current_col = 1
    @inbounds for Kcol in sorteddofs
        maxlookups = sym ? current_col : length(dofs)
        current_idx = 1
        for r in nzrange(K, Kcol)
            Kerow = permutation[current_idx]
            if K.rowval[r] == dofs[Kerow]
                Kecol = permutation[current_col]
                K.nzval[r] += Ke[Kerow, Kecol]
                current_idx += 1
            end
            current_idx > maxlookups && break
        end
        if current_idx <= maxlookups
            error("some row indices were not found")
        end
        current_col += 1
    end
end

# Sort utilities

function sortperm2!(B, ii)
   @inbounds for i = 1:length(B)
      ii[i] = i
   end
   quicksort!(B, ii)
   return
end

function quicksort!(A, order, i=1,j=length(A))
    @inbounds if j > i
        if  j - i <= 12
           # Insertion sort for small groups is faster than Quicksort
           InsertionSort!(A, order, i, j)
           return A
        end

        pivot = A[div(i+j,2)]
        left, right = i, j
        while left <= right
            while A[left] < pivot
                left += 1
            end
            while A[right] > pivot
                right -= 1
            end
            if left <= right
                A[left], A[right] = A[right], A[left]
                order[left], order[right] = order[right], order[left]

                left += 1
                right -= 1
            end
        end  # left <= right

        quicksort!(A,order, i,   right)
        quicksort!(A,order, left,j)
    end  # j > i

    return A
end

function InsertionSort!(A, order, ii=1, jj=length(A))
    @inbounds for i = ii+1 : jj
        j = i - 1
        temp  = A[i]
        itemp = order[i]

        while true
            if j == ii-1
                break
            end
            if A[j] <= temp
                break
            end
            A[j+1] = A[j]
            order[j+1] = order[j]
            j -= 1
        end

        A[j+1] = temp
        order[j+1] = itemp
    end  # i
    return
end


text:
function test2()
  filename = "test2.asc"
  # exectuablepath = null string will not run LTspice.exe.  Test parsing only.
  sim = LTspiceSimulation(filename,executablepath="")
  @test LTspice.does_circuitfilearray_file_match(sim)
  show(IOBuffer(),sim)
  @test sim["c"]≈2.0
  @test sim["b"]≈8.0
  @test sim["a"]≈10.0
  @test sim["d"]≈100.0
  @test sim["x"]≈1.00394
  @test sim["z"]≈0.019685
  @test sim["y"]≈0.984252

  @test(length(measurementvalues(sim))==3)
  @test(length(parametervalues(sim))==4)

  @test(length(keys(sim))==7)
  @test(length(values(sim))==7)
  @test(length(sim)==7)
  show(IOBuffer(),sim)
end
test2()


text:
"""
    compute_energy(a4, phi2, phi0, phiM2)

Returns the normalized energy density ``ϵ / ϕ₀⁴``.
"""
function compute_energy(a4, phi2, phi0, oophiM2)
    if phi0 < 1.e-9
        return -3//4 * a4
    else
        return -3//4 * a4 / phi0^4 - phi2 / phi0^3 .- 1//4 * oophiM2 .+ 7//36
    end
end

"""
    compute_flux(f2)

Returns the normalized flux ``J / ϕ₀⁴``.
"""
function compute_flux(f2, phi0)
    if phi0 < 1.e-9
        return f2
    else
        return f2 / phi0^4
    end
end

"""
    compute_pressure_x(a4, b14, b24, phi2, phi0, oophiM2)

Returns the normalized pressure in the x direction ``pₓ / ϕ₀⁴``.
"""
function compute_pressure_x(a4, b14, b24, phi2, phi0, oophiM2)
    if phi0 < 1.e-9
        return -a4/4 - b14 - b24
    else
        return -a4/(4*phi0^4) - b14/phi0^4 - b24/phi0^4 + phi2/(3*phi0^3) .+
            (-5//108 + 1//4 * oophiM2)
    end
end

"""
    compute_pressure_y(a4, b14, b24, phi2, phi0, oophiM2)

Returns the normalized pressure in the y direction ``p_y / ϕ₀⁴``.
"""
function compute_pressure_y(a4, b14, b24, phi2, phi0, oophiM2)
    if phi0 < 1.e-9
        return -a4/4 + b14 - b24
    else
        return -a4/(4*phi0^4) + b14/phi0^4 - b24/phi0^4 + phi2/(3*phi0^3) .+
            (-5//108 + 1//4*oophiM2)
    end
end

"""
    compute_pressure_z(a4, b24, phi2, phi0, oophiM2)

Returns the normalized pressure in the z direction ``p_z / ϕ₀⁴``.
"""
function compute_pressure_z(a4, b24, phi2, phi0, oophiM2)
    if phi0 < 1.e-9
        return -a4/4 + 2*b24
    else
        return -a4/(4*phi0^4) + 2*b24/phi0^4 + phi2/(3*phi0^3) .+
            (-5//108 + 1//4*oophiM2)
    end
end

"""
    compute_pressure_xy(a4, b24, phi2, phi0, oophiM2)

Returns the normalized pressure in the xy direction ``p_xy / ϕ₀⁴``.
"""
function compute_pressure_xy(g4, phi0)
    if phi0 < 1.e-9
        return -g4
    else
        return -g4 / phi0^4
    end
end

"""
    compute_Ophi(phi2, phi0, oophiM2)

Returns the normalized scalar operator VEV ``O_ϕ / ϕ₀³``.
"""
function compute_Ophi(phi2, phi0, oophiM2)
    if phi0 < 1.e-9
        return -2 * phi2
    else
        return -2 * phi2 / phi0^3 .+ (1//3 - oophiM2)
    end
end

function compute_temperature(Du_A_uAH, uAH, phi0)
    if phi0 < 1.e-9
        return (-0.5 * uAH .^ 2 .* Du_A_uAH)/(2*pi)
    else
        return (-0.5 * uAH .^ 2 .* Du_A_uAH)/(2*pi*phi0)
    end
end

function compute_entropy(S_uAH, uAH, phi0)
    if phi0 < 1.e-9
        return pi*S_uAH.^3
    else
        return pi*(S_uAH/phi0).^3
    end
end

function get_energy(ts::OpenPMDTimeSeries; it::Int, verbose::Bool=false)
    a4,   chart = get_field(ts, it=it, field="a4", verbose=verbose)
    phi2, chart = get_field(ts, it=it, field="phi2", verbose=verbose)

    phi0 = try
        ts.params["phi0"]
    catch e
        if isa(e, KeyError)
            0.0   # if "phi0" is not found in the params Dict, set phi0 = 0
        else
            throw(e)
        end
    end

    oophiM2 = try
        ts.params["oophiM2"]
    catch e
        if isa(e, KeyError)
            0.0   # if "oophiM2" is not found in the params Dict, set oophiM2 = 0
        else
            throw(e)
        end
    end

    en = compute_energy(a4, phi2, phi0, oophiM2)

    en, chart
end

function get_Jx(ts::OpenPMDTimeSeries; it::Int, verbose::Bool=false)
    fx2, chart = get_field(ts, it=it, field="fx2", verbose=verbose)

    phi0 = try
        ts.params["phi0"]
    catch e
        if isa(e, KeyError)
            0.0   # if "phi0" is not found in the params Dict, set phi0 = 0
        else
            throw(e)
        end
    end

    Jx = compute_flux(fx2, phi0)

    Jx, chart
end

function get_Jy(ts::OpenPMDTimeSeries; it::Int, verbose::Bool=false)
    fy2, chart = get_field(ts, it=it, field="fy2", verbose=verbose)

    phi0 = try
        ts.params["phi0"]
    catch e
        if isa(e, KeyError)
            0.0   # if "phi0" is not found in the params Dict, set phi0 = 0
        else
            throw(e)
        end
    end

    Jy = compute_flux(fy2, phi0)

    Jy, chart
end

function get_px(ts::OpenPMDTimeSeries; it::Int, verbose::Bool=false)
    a4,   chart = get_field(ts, it=it, field="a4", verbose=verbose)
    b14,  chart = get_field(ts, it=it, field="b14", verbose=verbose)
    b24,  chart = get_field(ts, it=it, field="b24", verbose=verbose)
    phi2, chart = get_field(ts, it=it, field="phi2", verbose=verbose)

    phi0 = try
        ts.params["phi0"]
    catch e
        if isa(e, KeyError)
            0.0   # if "phi0" is not found in the params Dict, set phi0 = 0
        else
            throw(e)
        end
    end

    oophiM2 = try
        ts.params["oophiM2"]
    catch e
        if isa(e, KeyError)
            0.0   # if "oophiM2" is not found in the params Dict, set oophiM2 = 0
        else
            throw(e)
        end
    end

    px = compute_pressure_x(a4, b14, b24, phi2, phi0, oophiM2)

    px, chart
end

function get_py(ts::OpenPMDTimeSeries; it::Int, verbose::Bool=false)
    a4,   chart = get_field(ts, it=it, field="a4", verbose=verbose)
    b14,  chart = get_field(ts, it=it, field="b14", verbose=verbose)
    b24,  chart = get_field(ts, it=it, field="b24", verbose=verbose)
    phi2, chart = get_field(ts, it=it, field="phi2", verbose=verbose)

    phi0 = try
        ts.params["phi0"]
    catch e
        if isa(e, KeyError)
            0.0   # if "phi0" is not found in the params Dict, set phi0 = 0
        else
            throw(e)
        end
    end

    oophiM2 = try
        ts.params["oophiM2"]
    catch e
        if isa(e, KeyError)
            0.0   # if "oophiM2" is not found in the params Dict, set oophiM2 = 0
        else
            throw(e)
        end
    end

    py = compute_pressure_y(a4, b14, b24, phi2, phi0, oophiM2)

    py, chart
end

function get_pz(ts::OpenPMDTimeSeries; it::Int, verbose::Bool=false)
    a4,   chart = get_field(ts, it=it, field="a4", verbose=verbose)
    b24,  chart = get_field(ts, it=it, field="b24", verbose=verbose)
    phi2, chart = get_field(ts, it=it, field="phi2", verbose=verbose)

    phi0 = try
        ts.params["phi0"]
    catch e
        if isa(e, KeyError)
            0.0   # if "phi0" is not found in the params Dict, set phi0 = 0
        else
            throw(e)
        end
    end

    oophiM2 = try
        ts.params["oophiM2"]
    catch e
        if isa(e, KeyError)
            0.0   # if "oophiM2" is not found in the params Dict, set oophiM2 = 0
        else
            throw(e)
        end
    end

    pz = compute_pressure_z(a4, b24, phi2, phi0, oophiM2)

    pz, chart
end

function get_pxy(ts::OpenPMDTimeSeries; it::Int, verbose::Bool=false)
    g4,   chart = get_field(ts, it=it, field="g4", verbose=verbose)

    phi0 = try
        ts.params["phi0"]
    catch e
        if isa(e, KeyError)
            0.0   # if "phi0" is not found in the params Dict, set phi0 = 0
        else
            throw(e)
        end
    end

    pxy = compute_pressure_xy(g4, phi0)

    pxy, chart
end

function get_Ophi(ts::OpenPMDTimeSeries; it::Int, verbose::Bool=false)
    phi2,  chart = get_field(ts, it=it, field="phi2", verbose=verbose)

    phi0 = try
        ts.params["phi0"]
    catch e
        if isa(e, KeyError)
            0.0   # if "phi0" is not found in the params Dict, set phi0 = 0
        else
            throw(e)
        end
    end

    oophiM2 = try
        ts.params["oophiM2"]
    catch e
        if isa(e, KeyError)
            0.0   # if "oophiM2" is not found in the params Dict, set oophiM2 = 0
        else
            throw(e)
        end
    end

    Ophi = compute_Ophi(phi2, phi0, oophiM2)

    Ophi, chart
end

#Remember that when doing get_field with a given it we always take the file whose it is the closest to the specified one.
function get_temperature(ts_const::OpenPMDTimeSeries, ts_diag::OpenPMDTimeSeries; it::Int, verbose::Bool=false)
    cmax = 0
    i    = 0

    sigma, _ = get_field(ts_diag, it=it, field="sigma", verbose=verbose)
    uAH      = 1 ./sigma

    phi0 = try
        ts_const.params["phi0"]
    catch e
        if isa(e, KeyError)
            0.0   # if "phi0" is not found in the params Dict, set phi0 = 0
        else
            throw(e)
        end
    end

    while cmax == 0
        i += 1
        try #We find what is the deepest grid.
            get_field(ts_const, it=it, field="A c=$i", verbose=verbose)
        catch
            cmax = i-1
            @warn "maximum value for c is $(i-1)"
        end
    end
    A, chart = get_field(ts_const, it=it, field="A c=$cmax", verbose=verbose)
    u, x, y  = chart[:]
    Nx       = length(x)
    Ny       = length(y)
    Nu       = length(u)
    der      = ChebDeriv(1,u[1],u[end],length(u))
    Du       = der.D
    uinterp  = ChebInterpolator(u[1],u[end],length(u))
    Du_A_uAH = similar(sigma)
    for j in 1:Ny
        for i in 1:Nx
            Du_A = Du*A[:,i,j]
            Du_A_uAH[1,i,j] = uinterp(view(Du_A,:))(uAH[1,i,j])
        end
    end

    T = compute_temperature(Du_A_uAH, uAH, phi0)

    T, chart
end

function get_entropy(ts_const::OpenPMDTimeSeries, ts_diag::OpenPMDTimeSeries; it::Int, verbose::Bool=false)
    cmax = 0
    i    = 0

    sigma, _ = get_field(ts_diag, it=it, field="sigma", verbose=verbose)
    uAH      = 1 ./sigma

    phi0 = try
        ts_const.params["phi0"]
    catch e
        if isa(e, KeyError)
            0.0   # if "phi0" is not found in the params Dict, set phi0 = 0
        else
            throw(e)
        end
    end

    while cmax == 0
        i += 1
        try #We find what is the deepest grid.
            get_field(ts_const, it=it, field="S c=$i", verbose=verbose)
        catch
            cmax = i-1
            @warn "maximum value for c is $(i-1)"
        end
    end
    S, chart = get_field(ts_const, it=it, field="S c=$cmax", verbose=verbose)
    u, x, y = chart[:]
    Nx      = length(x)
    Ny      = length(y)
    Nu      = length(u)
    uinterp = ChebInterpolator(u[1],u[end],length(u))
    S_uAH   = similar(sigma)
    for j in 1:Ny
        for i in 1:Nx
            S_uAH[1,i,j] = uinterp(view(S, :,i,j))(uAH[1,i,j])
        end
    end
    entr = compute_entropy(S_uAH, uAH, phi0)

    entr, chart
end

#=
Computes the fluid velocity and the local frame (diagonal) stress tensor. p3 is just pz. I think that
ut will always be 1, if not we can always hand ux/ut and uy/ut, which diretly are the speeds measured
in the LAB frame, i.e. dx/dt and dy/dt.
=#
function compute_local_VEVs(e, Jx, Jy, px, py, pxy)
    T = [e -Jx -Jy; -Jx px pxy; -Jy pxy py]
    values, vectors = eigen(T)
    norms = -vectors[1,:].^2 + vectors[2,:].^2 + vectors[3,:].^2
    index = findfirst(==(1), norms.<0)
    ut = sign(vectors[1,index])*vectors[1,index]/(-norms[index])^0.5
    ux = sign(vectors[1,index])*vectors[2,index]/(-norms[index])^0.5
    uy = sign(vectors[1,index])*vectors[3,index]/(-norms[index])^0.5
    e_local = values[index]
    p1_local = values[1:end .!=index][1]
    p2_local = values[1:end .!=index][2]

    return ut, ux, uy, e_local, p1_local, p2_local
end


text:
type SDF_MYD_BW <: SDFDynamicSystem
    omg::Float64
    zta::Float64
    kpa::Float64
    phi::Float64
    alpha::Float64
    u0::Float64
    v0::Float64
    z0::Float64
    n::Int64
    bta::Float64

    SDF_MYD_BW(omg::Float64, zta::Float64, kpa::Float64, phi::Float64, alpha::Float64, n::Int64, bta::Float64) = new(omg, zta, kpa, phi, alpha, 0.0, 0.0, 0.0, n, bta)
    SDF_MYD_BW(omg::Float64, zta::Float64, kpa::Float64, phi::Float64, alpha::Float64) = new(omg, zta, kpa, phi, alpha, 0.0, 0.0, 0.0, 5, 0.5)
    SDF_MYD_BW(omg::Float64, zta::Float64, kpa::Float64, phi::Float64) = new(omg, zta, kpa, phi, 0.0, 0.0, 0.0, 0.0, 5, 0.5)
end

function systemEquation(s::SDF_MYD_BW, p::Function)
    function dydt(y::AbstractArray{Float64,1}, t::Float64)
        u = y[1]
        v = y[2]
        z = y[3]
        omg2 = s.omg*s.omg
        dvdt = p(t)-2.0*s.zta*s.omg*v-(s.alpha*s.kpa*omg2*u+(1.0-s.alpha)*s.phi*g*z)-omg2*u
        dzdt = v - (s.bta*abs(v)*abs(z)^(s.n-1)*z + (1.0-s.bta)*v*abs(z)^(s.n))
        dzdt *= s.kpa*omg2/(s.phi*g)
        [v; dvdt; dzdt]
    end
    dydt
end

function response(s::SDF_MYD_BW, p::Function, t::AbstractArray{Float64,1})
    y0 = [s.u0; s.v0; s.z0]

    f = systemEquation(s, p)

    y = ode4(f, y0, t)

    u = y[:,1]
    v = y[:,2]
    z = y[:,3]
    omg2 = s.omg*s.omg
    fdz = 2.0*s.zta*s.omg*v
    fdp = s.alpha*s.kpa*omg2*u+(1.0-s.alpha)*s.phi*g*z
    a = map(p,t)-fdz-fdp-s.omg*s.omg*u

    [u v a fdz fdp]
end

max_stars_repo_name: Panchatantra/DSR

max_issues_repo_name: Panchatantra/DSR

max_forks_repo_name: Panchatantra/DSR


text:
const _DEFAULT_PICTURE_KWARGS =
    OrderedDict(
        :clip => nothing
    )

struct Picture <: GraphicElement2D
    elements::Array{<:GraphicElement,1}
    options::Array{Any,1}
end

Picture(;kwargs...) = Picture(GraphicElement[],kwargs...)

Picture(elements::Array{<:GraphicElement,1};kwargs...) =
        Picture(elements,collect(kwargs))

Picture(element::GraphicElement;kwargs...) = Picture([element];kwargs...)

function Picture(P::Picture;kwargs...)
    Picture(P.elements,updatedvals(
                    _DEFAULT_PICTURE_KWARGS,kwargs))
end


text:
# n_x = 5
# n_h = 4
# n_y = 1

# W1 = randn(n_h, n_x)
# println(W1)

import Random
Random.seed!(210)

# function initialize_parameters(layer_dims)
#     parameters = Dict{String, Array{Float64}}()
#     for i = 2:length(layer_dims)
#         parameters[string("W", i-1)] = randn(layer_dims[i], layer_dims[i-1]) * 0.01
#         parameters[string("b", i-1)] = zeros(layer_dims[i], 1)
#     end
#     return parameters
# end

include("NeuralNetwork.jl")

parameters = initialize_parameters([5 10 1])

function printdict(dict)
    for (key, value) in dict
        println(key, " ", size(value))
    end
end
function printdict_full(dict)
    for (key, value) in dict
        println(key, " ", (value))
    end
end
printdict(parameters)

# using LinearAlgebra

# function sigmoid(X)
#     return 1 ./ (1 .+ exp.(-X))
# end
# function sigmoid(x)
#     return 1 / (1 + exp(-x))
# end

# function relu(x)
#     if x>0
#         return x
#     else
#         return 0
#     end
# end

# function forward_prop(X::Array{Float64}, parameters::Dict{String, Array{Float64}}, layer_dims::Array{Int}, activations=Nothing)
#     num_layers = length(layer_dims)

#     if activations === Nothing
#         activations = Dict()
#         for i = 1:num_layers-1
#             # println("relu")
#             activations[i] = relu
#         end
#         # println(sigmoid)
#         activations[num_layers] = sigmoid
#     end
#     # println(activations)
#     caches = Dict()

#     Z1 = parameters[string("W", 1)] * X .+ parameters[string("b", 1)]
#     A1 = activations[1].(Z1)
#     caches[string("Z", 1)] = Z1
#     caches[string("A", 1)] = A1

#     Ai = Nothing
#     for i = 2:num_layers-1
#         Zi = parameters[string("W", i)] * caches[string("A", i-1)] .+ parameters[string("b", i)]
#         Ai = activations[i].(Zi)
#         caches[string("Z", i)] = Zi
#         caches[string("A", i)] = Ai
#     end

#     return Ai, caches
# end

X = randn(5, 100)
a2, caches = forward_prop(X, parameters, [5 10 1], (relu, sigmoid))
# a2, caches = forward_prop(X, parameters, [5 10 1])
println(a2)
# println(caches)

# function cost_binary(Y, Ŷ)
#     @assert length(Y) == length(Ŷ)
#     m = length(Y)

#     cost = - sum(Y .* log.(Ŷ) .+ (1 .- Y) .* log.(1 .- Ŷ)) / m
# end

Y = [0.777085  0.248755  0.660713  0.37982  0.978839  0.182363  0.07204  0.855448  0.37393  0.40736]
@assert cost_binary(Y, Y) == 0.4982477807320199

# function sigmoid_back(x)
#     sig = sigmoid(x)
#     return sig*(1-sig)
# end

# function relu_back(x)
#     if x>0
#         return 1
#     else
#         return 0
#     end
# end

# function backward_prop(Y::Array{Float64}, Ŷ::Array{Float64}, parameters::Dict{String, Array{Float64}}, caches::Dict{String, Array{Float64}}, layer_dims::Array{Int}, activations=Nothing)

#     num_layers = length(layer_dims)
#     @assert length(Y) == length(Ŷ)
#     m = size(Y)[2]

#     if activations === Nothing
#         activations = Dict()
#         for i = 1:num_layers-1
#             activations[i] = relu_back
#         end
#         activations[num_layers] = sigmoid_back
#     else
#         activations_orig = activations
#         activations = []
#         for activation in activations_orig
#             push!(activations, @eval ($(Symbol("$activation", "_back"))))
#         end
#     end
#     # println(activations)

#     dA = sum(.- Y ./ Ŷ .+ (1 .- Y) ./ (1 .- Ŷ))

#     grads = Dict{String, Array{Float64}}()

#     for l in num_layers-1:-1:1
#         dZ = dA .* activations[l].(caches[string("Z", l)])
#         grads[string("dw", l)] = 1/m .* (dZ * transpose(caches[string("A", l-1)]))
#         grads[string("db", l)] = 1/m .* sum(dZ, dims=2)
#         dA = transpose(parameters[string("W", l)]) * dZ
#     end

#     return grads
# end

# grads = backward_prop(randn(1, 100), a2, parameters, caches, [5 10 1])
# printdict(grads)
# @assert size(grads["dw1"]) == size(parameters["W1"])
grads = backward_prop(randn(1, 100), a2, parameters, caches, [5 10 1], [relu, sigmoid])
printdict(grads)
@assert size(grads["dw1"]) == size(parameters["W1"])

# function update_parameters(parameters::Dict{String, Array{Float64}}, grads::Dict{String, Array{Float64}}, layer_dims::Array{Int}, learning_rate::Float64)
#     num_layers = length(layer_dims)
#     for l = 1:num_layers-1
#         parameters[string("W", l)] += learning_rate .* grads[string("dw", l)]
#         parameters[string("b", l)] += learning_rate .* grads[string("db", l)]
#     end
#     return parameters
# end

printdict_full(parameters)
update_parameters(parameters, grads, [5 10 1], 0.01)
printdict_full(parameters)

# function neural_network_dense(X, Y, layer_dims::Array{Int}, num_iterations::Int, learning_rate::Float64, activations=Nothing)
#     num_layers = length(layer_dims) # calculate number of layers
#
#     # if activations are not given, assume that all hidden layers have relu and output layer has sigmoid
#     if activations === Nothing
#         activations = Array{Function}(undef, num_layers-1)
#         for i = 1:num_layers-2
#             activations[i] = relu
#         end
#         activations[num_layers-1] = sigmoid
#     end
#     activations = Tuple(activations)
#     println(activations)
#     activations_back = []
#     for activation in activations
#         push!(activations_back, @eval ($(Symbol("$activation", "_back"))))
#     end
#     activations_back = Tuple(activations_back)
#     println(activations_back)
#
#     parameters = initialize_parameters(layer_dims)
#
#     for iteration = 1:num_iterations
#         Ŷ, caches = forward_prop(X, parameters, layer_dims, activations)
#         grads = backward_prop(Y, Ŷ, parameters, caches, layer_dims, activations_back)
#         parameters = update_parameters(parameters, grads, layer_dims, learning_rate)
#
#         if iteration % 1 == 0
#             cost = cost_binary(Y, Ŷ)
#             println("Cost at iteration $iteration is $cost")
#         end
#     end
#
#     return parameters
# end

X = randn(10, 100)
Y = randn(1, 100)
layer_dims = [10 10 1]
neural_network_dense(X, Y, layer_dims, 100, 0.05)


text:
question_css = """
<style>
div.question {
  padding-left: 20px;
  padding-right: 20px;
  padding-top: 10px;
  padding-bottom: 10px;
  border: 3px dotted lightgrey;
  border-radius: 15px;
  background: #F8F8F8;
}
</style>"""

fyi_css = html"""
<style> pluto-output div.admonition.info .admonition-title {
        background: rgb(161, 161, 161);
      }

      pluto-output div.admonition.info {
        background: rgba(161, 161, 161, 0.2);
        border: 5px solid rgb(161, 161, 161);
      }include("styles.jl")

<style> pluto-output div.admonition.bomb .admonition-title {
						background: rgb(226, 157, 148);
						animation:blinkingBox 1.5s infinite;
					}

					pluto-output div.admonition.bomb {
						background: rgba(226, 157, 148, 0.2);
						border: 5px solid rgb(226, 157, 148);
						animation:blinkingBox 1.5s infinite;
					}

			@keyframes blinkingBox{
				0%{     visibility: hidden;    }
				30%{    visibility: hidden; }
				31%{    visibility: visible; }
				99%{    visibility: visible;  }
				100%{   visibility: hidden;     }
			}
</style>"""

max_stars_repo_name: Beramos/PlutoTeaching

max_issues_repo_name: Beramos/PlutoTeaching

max_forks_repo_name: Beramos/PlutoTeaching


text:
using DrWatson
@quickactivate :FractalDimension # uses DynamicalSystems, PyPlot
include(srcdir("style.jl"))

using DelimitedFiles

name = "double_pendulum2.csv"
file = datadir("experimental", name)
M = readdlm(file, ',')

x = M[:, 5][1:15000]
y = M[:, 6][1:15000]

fig, axs = subplots(1,2; figsize = (10, 5))
axs[1].plot(x)
axs[2].plot(x,y)
axs[1].plot(y)

writedlm(datadir("experimental", "doublependulum.txt"), [x y])

estimate_embedding_plot(x)
estimate_embedding_plot(y)
# Mutual information does NOT have a clear minimum
# the value of 50 is given approximately
# Cao's method provides embedding dimension of 5

_, τ_vals_x = pecuzal_embedding(x .+ 1e-9randn(length(x)); τs = 0:100, w = 5)
_, τ_vals_y = pecuzal_embedding(y .+ 1e-9randn(length(x)); τs = 0:100, w = 5)
# Astonishingly Pecuzal method gives very similar τ for both input timeseries
τx = [0, 51, 25, 39, 12]
τx = [0, 60, 30, 46, 14]

# Let's see how it does for the dual input
A = Dataset(x .+ 1e-9randn(length(x)), y .+ 1e-9randn(length(x)))
_, τ_vals_d, ts_vals = pecuzal_embedding(A; τs = 0:100, w = 5)

# It gives a 6dimension embedding
τd = [0, 61, 31, 47, 31, 0]
jd = [2, 2, 2, 2, 1, 1]

# # PECUZAL method gives 4d embedding with τ = [0, 6, 3, 14]
#
# # %%
# writedlm(datadir("experimental", "roessler_all.txt"), M)
# writedlm(datadir("experimental", "roessler_embed.txt"), x2)
#
# # %% Okay, this will take tremendous amount of computation time,
# # using all 28 timeseries...
# # Y, τ_vals, ts_vals, = pecuzal_embedding(X1)

max_stars_repo_name: Datseris/FractalDimension

max_issues_repo_name: Datseris/FractalDimension

max_forks_repo_name: Datseris/FractalDimension


text:
struct Tendencies{NF<:AbstractFloat}
    vor_tend     ::Array{Complex{NF},3}      # Vorticity of horizontal wind field
    div_tend     ::Array{Complex{NF},3}      # Divergence of horizontal wind field
    Tabs_tend    ::Array{Complex{NF},3}      # Absolute temperature [K]
    logp0_tend   ::Array{Complex{NF},2}      # Log of surface pressure [log(Pa)]
    humid_tend   ::Array{Complex{NF},3}      # Specific humidity [g/kg]
end

struct GridVariables{NF<:AbstractFloat}
    u       ::Array{NF,3}       # Zonal velocity [m/s]
    v       ::Array{NF,3}       # Meridional velocity [m/s]
    Tabs    ::Array{NF,3}       # Absolute temperature [K]
    logp0   ::Array{NF,2}       # Log of surface pressure [log(Pa)]
    geopot  ::Array{NF,3}       #
    humid   ::Array{NF,3}      # Specific humidity [g/kg]
end

"""Struct holding the diagnostic variables."""
struct DiagnosticVariables{NF<:AbstractFloat}
    tendencies  ::Tendencies{NF}
    gridvars    ::GridVariables{NF}
end


text:
#**************************************************************************************
# test_Output_Interface.jl
# ======================== part of the GeoEfficiency.jl package.
#
#
#
#**************************************************************************************

@testset "GeoEfficiency.checkResultsDirs" begin
	@test G.checkResultsDirs() === nothing
end #testset_checkResultsDirs

@debug("max_batch")
@testset "max_batch" begin
	@test max_batch(-1) 	== G._max_batch
	@test max_batch(100) 	== G._max_batch
	@test max_batch() 		== G.max_display
end

@debug("calc - CylDetector, WellDetector")
let pnt::Point = Point(1),
	cylDet::CylDetector = Detector(5, 10),
	wellDet::WellDetector = Detector(5, 4, 3.2, 2)

	@testset "calc - CylDetector, WellDetector" for
	SrcRadius = Real[1, 1//2, e, pi, 1.0],
	SrcLength = Real[1, 1//2, e, pi, 1.0]

		@test calc(cylDet,  (pnt, SrcRadius , SrcLength))    === nothing
		@test calc(wellDet, (pnt, SrcRadius , SrcLength))    === nothing

		@test_throws 	G.NotImplementedError  	H.@console(calc(cylDet), 0, 15, 0)
		@test_throws 	G.InValidGeometry  		H.@console(calc(cylDet), -1, 0, 0)
	end #testset_calc
end #let

@debug("calcN")
@testset "calcN - $console_input" for
console_input = ["4 0 1 2 Q",
				"4 0 1 3 " * "d " * "4 0 1 4 Q",
				"4 0 1 5 " * "n " * "11 6 0 " * "4 0 1 6 Q"]

	@test H.exec_console_unattended(calcN, console_input, Fn_ARGs =[Detector(5, 10)])  === nothing
	@test H.exec_console_unattended(calcN,  console_input, Fn_ARGs =[Detector(eps())])  === nothing
	@test H.exec_console_unattended(calcN, "13 7 0 " * console_input)      === nothing
	@test H.exec_console_unattended(calcN, "13 7 0 Q " * console_input)      === nothing
	#@test H.exec_console_unattended(calcN, "-13 -7 0 Q " * console_input)      === nothing
end #testset_calcN

@debug("writecsv_head")
@testset "writecsv_head" begin
	# `writecsv_head` 		tests in the `reading_from _CSV` testset in `test_Input_Batch`.
	# `writecsv_head_any` 	tests in the `reading_from _CSV` testset in `test_Input_Batch`.
end #testset_writecsv_head

@debug("GeoEfficiency._batch")
@testset "_batch, point is $isSrcPoint" for	isSrcPoint = [true, false]

	rtrn = G._batch(Val(isSrcPoint),  CylDetector(eps(0.1)), [0.0], [0.0], [0.0], [0.0])
	@test typeof(rtrn) <: Tuple{Detector, Matrix{Float64}, String}
	@test rtrn[2][end] ≈ 0.5
	rm(rtrn[end]; force=true)

	rtrn = G._batch(Val(isSrcPoint), CylDetector(eps(0.1)), [0.0], [1.0], [0.0],[0.0])
	@test typeof(rtrn) <: Tuple{Detector, Matrix{Float64}, String}
	@test rtrn[2][end] |> isnan
	rm(rtrn[end]; force=true)

end #testset_GeoEfficiency._batch

@debug("batch")
@testset "batch" begin

local paths::Vector{String} = String[]
local every_path::Vector{String} = String[]
max_batch(-1)		# display the batch calculation results on the console

	@testset "batch(<:Detector)" begin
		acylDetector1::CylDetector = CylDetector(eps(0.1))

		path::String = batch(acylDetector1, [0.0])
		@test occursin( G.id(acylDetector1), path)
		@test readdlm(path,',')[2,end] ≈ 0.5
		push!(every_path, path)

		path = batch(acylDetector1, [0.0], [0.0], [0.0],[0.0],false)
		@test occursin(G.id(acylDetector1), path)
		@test readdlm(path,',')[2,end] ≈ 0.5
		push!(every_path, path)
	end #testset_batch(<:Detector)

max_batch(Inf)		# prevent the display of batch calculation results on the console
	let acylDetector2::CylDetector = CylDetector(eps(0.2))

		paths = batch([acylDetector2], [0.0]) # in fact `paths` is a one element vector
		@test occursin.(G.id(acylDetector2), paths) |> any
		append!(every_path, paths)

		paths = batch([acylDetector2], [0.0], [0.0], [0.0],[0.0],false)
		@test occursin.(G.id(acylDetector2), paths) |> any
		append!(every_path, paths)
	end #let

	let aBDetector = BoreDetector(eps(0.5), eps(0.4), eps(0.2))

		paths = batch([aBDetector], [0.0])
		@test occursin.(G.id(aBDetector), paths) |> any
		append!(every_path, paths)
	end #let

	let aWDetector = WellDetector(eps(0.5), eps(0.4), eps(0.2), eps(0.1))

		paths = batch([aWDetector], [0.0])
		@test occursin.(G.id(aWDetector), paths) |> any
		append!(every_path, paths)

		paths = batch([aWDetector], [0.0], [0.0], [0.0],[0.0],false)
		@test occursin.(G.id(aWDetector), paths) |> any
		append!(every_path, paths)
	end #let

	let acylDetector3::CylDetector 	= CylDetector(eps(0.3)),
		aBDetector::BoreDetector	= BoreDetector(6, eps(), eps(0.5)),
		aWDetector::WellDetector 	= WellDetector(eps(), eps(0.4), eps(0.1), eps(0.1))

		paths = batch([acylDetector3, aWDetector], [0.0])
		@test occursin.(G.id(aWDetector), paths) |> any
		append!(every_path, paths)

		paths = batch([acylDetector3, aWDetector], [0.0], [0.0], [0.0],[0.0],false)
		@test occursin.(G.id(acylDetector3), paths) |> any
		@test occursin.(G.id(aWDetector)  , paths) |> any
		append!(every_path, paths)

		paths = batch([acylDetector3, aBDetector], [0.0])
		@test occursin.(G.id(acylDetector3), paths) |> any
		append!(every_path, paths)

		paths = batch([acylDetector3, aBDetector], [0.0], [0.0], [0.0],[0.0],false)
		@test occursin.(G.id(acylDetector3), paths) |> any
		@test occursin.(G.id(aBDetector)  , paths) |> any
		append!(every_path, paths)
	end #let

	let aBDetector::BoreDetector = BoreDetector(eps(0.5), eps(0.4), eps(0.1)),
		aWDetector::WellDetector = WellDetector(eps(0.5), eps(0.4), eps(0.1), eps(0.1))

		paths = batch([aBDetector, aWDetector], [0.0])
		@test occursin.(G.id(aBDetector), paths) |> any
		@test occursin.(G.id(aWDetector), paths) |> any
		append!(every_path, paths)

		paths = batch([aBDetector, aWDetector], [0.0], [0.0], [0.0],[0.0],false)
		@test occursin.(G.id(aBDetector), paths) |> any
		@test occursin.(G.id(aWDetector), paths) |> any
		append!(every_path, paths)
	end #let

	let acylDetector = CylDetector(6, eps(0.65)),
		aBDetector::BoreDetector  = BoreDetector(7, eps(), eps(0.5)),
		aWDetector::WellDetector  = WellDetector(8, eps(), eps(0.5), eps(0.1))

		temppaths::Vector{String} = batch([acylDetector, aBDetector, aWDetector], [0.0])
		@test occursin.(G.id(acylDetector), temppaths) |> any
		@test occursin.(G.id(aBDetector), temppaths) |> any
		@test occursin.(G.id(aWDetector), temppaths) |> any
	chmod.(temppaths, 0o100444)		#make path read only

		paths = batch([acylDetector, aBDetector, aWDetector], [0.0])
		@test occursin.("_" * G.id(acylDetector), paths) |> any
		@test occursin.("_" * G.id(aBDetector)  , paths) |> any
		@test occursin.("_" * G.id(aWDetector)  , paths) |> any
	append!(every_path, paths)
	try
		append!(every_path, chmod.(temppaths, 0o777))   #
	catch err
		append!(every_path, temppaths)
		chmod.(temppaths, 0o777)
	end

		temppaths = batch([acylDetector, aBDetector, aWDetector], [0.0], [0.0], [0.0],[0.0], false)
		@test occursin.(G.id(acylDetector), temppaths) |> any
		@test occursin.(G.id(aBDetector)  , temppaths) |> any
		@test occursin.(G.id(aWDetector)  , temppaths) |> any
	chmod.(temppaths, 0o100444)

		paths = batch([acylDetector, aBDetector, aWDetector], [0.0], [0.0], [0.0],[0.0], false)
		@test occursin.("_" * G.id(acylDetector), paths) |> any
		@test occursin.("_" * G.id(aBDetector)  , paths) |> any
		@test occursin.("_" * G.id(aWDetector)  , paths) |> any
	append!(every_path, paths)
	try
		append!(every_path, chmod.(temppaths, 0o777))   #
	catch err
		append!(every_path, temppaths)
		chmod.(temppaths, 0o777)
	end
	end #let

	let acylDetector::CylDetector = CylDetector(5, eps()),
		aBDetector::BoreDetector  = BoreDetector(5, eps(), eps(0.5)),
		aWDetector::WellDetector  = WellDetector(5, eps(), eps(0.5), eps(0.1))

		@test append!(every_path, batch([eps() 0], [0.0]))|> eltype === String
		@test append!(every_path, batch([eps() 0], [0.0], [0.0], [0.0], [0.0], false)) |> eltype === String
		@test append!(every_path, batch([1.0 eps()], [0.0]))|> eltype === String
		@test append!(every_path, batch([1.0 eps()], [0.0], [0.0], [0.0], [0.0], false)) |> eltype === String
		@test append!(every_path, batch([1//2 eps()], [0.0]))|> eltype === String
		@test append!(every_path, batch([1//2 eps()], [0.0], [0.0], [0.0], [0.0], false)) |> eltype === String
		@test append!(every_path, batch([1//2 eps(0.5)], [0.0]))|> eltype === String
		@test append!(every_path, batch([1//2 eps(0.5)], [0.0], [0.0], [0.0], [0.0], false)) |> eltype === String
		@test append!(every_path, batch([e pi], [0.0]))|> eltype === String
		@test append!(every_path, batch([e pi], [0.0], [0.0], [0.0], [0.0], false)) |> eltype === String

		@test append!(every_path, batch([5.0 4 3.1], [0.0]))|> eltype === String
		@test append!(every_path, batch([5.0 4 3.1], [0.0], [0.0],[0.0],[0.0],false)) |> eltype === String
		@test append!(every_path, batch([5.0 4 3//1], [0.0]))|> eltype === String
		@test append!(every_path, batch([5.0 4 3//1], [0.0], [0.0],[0.0],[0.0],false)) |> eltype === String
		@test append!(every_path, batch([5.0 4 pi], [0.0]))|> eltype === String
		@test append!(every_path, batch([5.0 4 pi], [0.0], [0.0],[0.0], [0.0], false)) |> eltype === String

		@test append!(every_path, batch([5.0 4 3 2], [0.0]))|> eltype === String
		@test append!(every_path, batch([5.0 4 3 2], [0.0], [0.0], [0.0],[0.0],false))|> eltype === String

		#=@test append!(every_path, batch([acylDetector, aWDetector], [0.0]))|> eltype === String
		@test append!(every_path, batch([acylDetector, aWDetector], [0.0], [0.0], [0.0],[0.0],false))|> eltype === String

		@test append!(every_path, batch([acylDetector, aBDetector], [0.0]))|> eltype === String
		@test append!(every_path, batch([acylDetector, aBDetector], [0.0], [0.0], [0.0],[0.0],false))|> eltype === String

		@test append!(every_path, batch([aBDetector, aWDetector], [0.0]))|> eltype === String
		@test append!(every_path, batch([aBDetector, aWDetector], [0.0], [0.0], [0.0],[0.0],false))|> eltype === String
		=#
		@test append!(every_path, batch([acylDetector, aBDetector, aWDetector], [0.0]))|> eltype === String
		@test append!(every_path, batch([acylDetector, aBDetector, aWDetector], [0.0], [0.0], [0.0],[0.0],false))|> eltype === String
	end #let

	G.detector_info_from_csvFile()
	if  [0.0] != G.read_from_csvFile(G.srcHeights, G.dataDir)
		setSrcToPoint(true);
		@test_skip append!(every_path, batch())|> eltype === String

		if [0.0] != G.read_from_csvFile(G.srcRadii, G.dataDir)
			setSrcToPoint(false);
			@test_skip append!(every_path, batch())|> eltype === String
		end #if
	end	#if

	#rm.(batch([aWDetector], [0.0]))
	#rm.(batch([aWDetector], [0.0], [0.0], [0.0],[0.0],false))
	#=for cr = 0.2:0.1:0.7
		@test append!(every_path, batch([Detector(11, eps(cr))], [0.0], [0.0], [0.0],[0.0],false))|> eltype === String
		@test append!(every_path, batch([Detector(22, eps(cr))], collect(0.0:0.1:10), [0.0], [0.0],[0.0],false))|> eltype === String
		@test append!(every_path, batch([Detector(33, eps(cr))], [0.0]))|> eltype === String
	end #for
	=#
try
	rm.(every_path, force=true)
catch err
	@error err
end #try
end #testset_batch

max_stars_repo_name: DrKrar/GeoEfficiency

max_issues_repo_name: DrKrar/GeoEfficiency

max_forks_repo_name: DrKrar/GeoEfficiency


text:
"""
    FChain(layers...)

Your friendly `Chain` but uses `foldl` instead of `tuple-peeling` for faster
Zygote gradient generation. This might not work when higher order gradients are
needed -- so use carefully.

See: https://github.com/FluxML/Zygote.jl/issues/1126
"""
struct FChain{T}
    layers::T
    function FChain(xs...)
        xs = flatten_model(xs)
        return new{typeof(xs)}(xs)
    end
    FChain(xs::AbstractVector) = FChain(xs...)
end

function flatten_model(layers::Union{AbstractVector,Tuple})
    new_layers = []
    for l in layers
        f = flatten_model(l)
        if f isa Tuple || f isa AbstractVector
            append!(new_layers, f)
        elseif f isa FChain
            append!(new_layers, f.layers)
        elseif f isa Chain
            append!(new_layers, f.layers...)
        else
            push!(new_layers, f)
        end
    end
    return layers isa AbstractVector ? new_layers : Tuple(new_layers)
end

flatten_model(x) = x

flatten_model(x::Chain) = Chain(flatten_model(x.layers))

Flux.@forward FChain.layers Base.getindex, Base.length, Base.first, Base.last, Base.iterate, Base.lastindex, Base.keys

Flux.@functor FChain

# (c::FChain)(x) = foldl((y,f) -> f(y), c.layers; init=x) # NO, this forgets the gradient for x

(c::FChain)(x) = foldl((y, f) -> f(y), (x, c.layers...))

(c::FChain)(x...) = foldl((y, f) -> f(y...), (x, c.layers...))

# (c::FChain{<:AbstractVector})(x) = foldl((y, f) -> f(y), vcat([x], c.layers))

# (c::FChain{<:AbstractVector})(x...) = foldl((y, f) -> f(y...), vcat([x], c.layers))

Base.getindex(c::FChain, i::AbstractArray) = FChain(c.layers[i]...)

function Base.show(io::IO, c::FChain)
    print(io, "FChain(")
    Flux._show_layers(io, c.layers)
    print(io, ")")
end

Flux._show_layers(io, layers::AbstractVector) = join(io, layers, ", ")


text:
function outputFaceByNumber(k::Int64)
  (k == 1)  && return SVector{3, Float64}(-1,  0,  0)
  (k == 2)  && return SVector{3, Float64}(+1,  0,  0)
  (k == 3)  && return SVector{3, Float64}( 0, -1,  0)
  (k == 4)  && return SVector{3, Float64}( 0, +1,  0)
  (k == 5)  && return SVector{3, Float64}( 0,  0, -1)
  (k == 6)  && return SVector{3, Float64}( 0,  0, +1)
  error("a cube has 6 faces")
end

function outputCornerByNumber(k::Int64)
  (k == 1) && return SVector{3, Float64}(-1, -1, -1)
  (k == 2) && return SVector{3, Float64}(+1, -1, -1)
  (k == 3) && return SVector{3, Float64}(-1, +1, -1)
  (k == 4) && return SVector{3, Float64}(+1, +1, -1)
  (k == 5) && return SVector{3, Float64}(-1, -1, +1)
  (k == 6) && return SVector{3, Float64}(+1, -1, +1)
  (k == 7) && return SVector{3, Float64}(-1, +1, +1)
  (k == 8) && return SVector{3, Float64}(+1, +1, +1)
  error("a cube has 8 corners")
end

function outputEdgeByNumber(k::Int64)
  (k == 1)  && return SVector{3, Float64}( 0, -1, -1)
  (k == 2)  && return SVector{3, Float64}( 0, +1, -1)
  (k == 3)  && return SVector{3, Float64}( 0, -1, +1)
  (k == 4)  && return SVector{3, Float64}( 0, +1, +1)
  (k == 5)  && return SVector{3, Float64}(-1,  0, -1)
  (k == 6)  && return SVector{3, Float64}(+1,  0, -1)
  (k == 7)  && return SVector{3, Float64}(-1,  0, +1)
  (k == 8)  && return SVector{3, Float64}(+1,  0, +1)
  (k == 9)  && return SVector{3, Float64}(-1, -1,  0)
  (k == 10) && return SVector{3, Float64}(+1, -1,  0)
  (k == 11) && return SVector{3, Float64}(-1, +1,  0)
  (k == 12) && return SVector{3, Float64}(+1, +1,  0)
  error("a cube has 12 edges")
end

function face_corner_test(tt::TT_Cache, dir_1::SVector{3,Float64}, dir_2::SVector{3,Float64})
  aabb_1 = OBB(zero_three, box_1_extent, I3)
  aabb_2 = OBB(zero_three, box_2_extent, I3)
  vec_1 = dir_1 .* aabb_1.e
  vec_2 = dir_2 .* aabb_2.e
  R_total = RotMatrix(rotation_between(dir_2, dir_1))
  total_sep_dist = vec_1 + R_total * vec_2
  update_TT_Cache!(tt, total_sep_dist * (1 - tolerance), R_total)
  is_sect_pos = BB_BB_intersect(tt, aabb_1, aabb_2)
  update_TT_Cache!(tt, total_sep_dist * (1 + tolerance), R_total)
  is_sect_neg = BB_BB_intersect(tt, aabb_1, aabb_2)
  return is_sect_pos, is_sect_neg
end

function edge_edge_test(tt::TT_Cache, dir_1::SVector{3,Float64}, theta_axis::Float64, R_box_2::Rotation{3, Float64})
  aabb_1 = OBB(zero_three, one_three, I3)
  aabb_2 = OBB(zero_three, one_three, I3)
  vec_1 = dir_1 .* aabb_1.e
  R_total = RotMatrix(UnitQuaternion(AngleAxis(theta_axis, vec_1[1], vec_1[2], vec_1[3])) * R_box_2)
  total_sep_dist = vec_1 * 2
  update_TT_Cache!(tt, total_sep_dist * (1 - tolerance), R_total)
  is_sect_pos = BB_BB_intersect(tt, aabb_1, aabb_2)
  update_TT_Cache!(tt, total_sep_dist * (1 + tolerance), R_total)
  is_sect_neg = BB_BB_intersect(tt, aabb_1, aabb_2)
  return is_sect_pos, is_sect_neg
end

tt = TT_Cache()
box_1_extent = SVector(1.0, 2.0, 3.0)
box_2_extent = SVector(2.1, 2.2, 2.3)
one_three    = SVector(1.0, 1.0, 1.0)
zero_three   = SVector(0.0, 0.0, 0.0)
n_face       = 6
n_corner     = 8
n_edge       = 12
tolerance    = 1.0e-6
I3           = one(SMatrix{3,3,Float64,9})

@testset "OBB OBB Face Face Test" begin
  for k_face = 1:n_face
    for k_corner = 1:n_corner
      dir_face = outputFaceByNumber(k_face)
      dir_corner = outputCornerByNumber(k_corner)
      should_be_positive, should_be_negative = face_corner_test(tt, dir_face, dir_corner)
      @test should_be_positive && !should_be_negative
      should_be_positive, should_be_negative = face_corner_test(tt, dir_corner, dir_face)
      @test should_be_positive && !should_be_negative
    end
  end
end

@testset "OBB OBB Edge Edge Test" begin
  for k_edge_1 = 1:n_face
    for theta_axis = rand(15)*2*pi
      for theta_extra = 0:(pi/2):(2*pi)
        dir_edge = outputEdgeByNumber(k_edge_1)
        should_be_positive, should_be_negative = edge_edge_test(tt, dir_edge, theta_axis, RotX(theta_extra))
        @test should_be_positive && !should_be_negative
        should_be_positive, should_be_negative = edge_edge_test(tt, dir_edge, theta_axis, RotY(theta_extra))
        @test should_be_positive && !should_be_negative
        should_be_positive, should_be_negative = edge_edge_test(tt, dir_edge, theta_axis, RotZ(theta_extra))
        @test should_be_positive && !should_be_negative
      end
    end
  end
end


text:
export
Vector2, Vector2i, Vector2f, Vector2s, Vector2d,
Point2, Point2i, Point2f, Point2s, Point2d,
Vector3, Vector3i, Vector3f, Vector3s, Vector3d,
Point3, Point3i, Point3f, Point3s, Point3d,
Normal3, Normal3i, Normal3f, Normal3s, Normal3d

const CartesianPair{T,U} = CartesianTuple{2,T,U}
const CartesianTriple{T,U} = CartesianTuple{3,T,U}

const VectorLike{N,T,U<:VecLike} = CartesianTuple{N,T,U}
const Vect{N,T} = CartesianTuple{N,T,Vec}
const Point{N,T} = CartesianTuple{N,T,Pnt}
const Normal{N,T} = CartesianTuple{N,T,Nrm}

const Vector2{T} = Vect{2,T}
const Point2{T} = Point{2,T}
const Vector3{T} = Vect{3,T}
const Point3{T} = Point{3,T}
const Normal3{T} = Normal{3,T}

for type ∈ (:Vector2, :Point2, :Vector3, :Point3, :Normal3)
    @eval begin
        const $(Symbol("$(type)i")) = $(type){Int}
        const $(Symbol("$(type)f")){T<:AbstractFloat} = $(type){T}
        const $(Symbol("$(type)s")) = $(type){Float32}
        const $(Symbol("$(type)d")) = $(type){Float64}
    end
end


text:
if VERSION < v"1"
  println("julia $VERSION is not supprted")
  return
end

global evalon = true
#global evalon = false

global trycount = 0
global succount = 0

max_stars_repo_name: sazare/cheapViews

max_issues_repo_name: sazare/cheapViews

max_forks_repo_name: sazare/cheapViews


text:
# ---
# title: 903. Valid Permutations for DI Sequence
# id: problem903
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Divide and Conquer, Dynamic Programming
# link: <https://leetcode.com/problems/valid-permutations-for-di-sequence/description/>
# hidden: true
# ---
#
# We are given `S`, a length `n` string of characters from the set `{'D', 'I'}`.
# (These letters stand for "decreasing" and "increasing".)
#
# A  _valid permutation_  is a permutation `P[0], P[1], ..., P[n]` of integers
# `{0, 1, ..., n}`, such that for all `i`:
#
#   * If `S[i] == 'D'`, then `P[i] > P[i+1]`, and;
#   * If `S[i] == 'I'`, then `P[i] < P[i+1]`.
#
# How many valid permutations are there?  Since the answer may be large,
# **return your answer modulo`10^9 + 7`**.
#
#
#
# **Example 1:**
#
#
#
#     Input: "DID"
#     Output: 5
#     Explanation:
#     The 5 valid permutations of (0, 1, 2, 3) are:
#     (1, 0, 3, 2)
#     (2, 0, 3, 1)
#     (2, 1, 3, 0)
#     (3, 0, 2, 1)
#     (3, 1, 2, 0)
#
#
#
#
# **Note:**
#
#   1. `1 <= S.length <= 200`
#   2. `S` consists only of characters from the set `{'D', 'I'}`.
#
#
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
"""
    CooperativeIDM <: DriverModel{LaneFollowingAccel}
The cooperative IDM (c-IDM) is a rule based driver model for merging scenarios.
It controls the longitudinal actions of vehicles on the main lane.
A cooperation level `c` controls how the vehicles reacts to the merging vehicle.
When `c=0` the vehicle ignores the merging vehicle. When `c=1` the vehicle considers the merging
vehicle as its front car when TTC_ego > TTC_mergin_vehicle. When it is not considering the merging vehicle,
the car follows the IntelligentDriverModel.

# Fields
    - `env::MergingEnvironment = MergingEnvironment(main_lane_angle = 0.0, merge_lane_angle = pi/6)` the merging environment
    - `idm::IntelligentDriverModel = IntelligentDriverModel(v_des = env.main_lane_vmax, d_cmf = 2.0, d_max=2.0, T = 1.5, s_min = 2.0, a_max = 2.0)` the default IDM
    - `c::Float64 = 0.0` the cooperation level
    - `fov::Float64 = 20.0` [m] A field of view, the merging vehicle is not considered if it is further than `fov`
"""
@with_kw mutable struct CooperativeIDM <: DriverModel{LaneFollowingAccel}
    env::MergingEnvironment = MergingEnvironment(main_lane_angle = 0.0, merge_lane_angle = pi/6)
    idm::IntelligentDriverModel = IntelligentDriverModel(v_des = env.main_lane_vmax,
                                                         d_cmf = 2.0,
                                                         d_max=2.0,
                                                         T = 1.5,
                                                         s_min = 2.0,
                                                         a_max = 2.0)
    c::Float64 = 0.0 # cooperation level
    fov::Float64 = 20.0 # when to consider merge car
    # internals
    a::Float64 = 0.0
    a_merge::Float64 = 0.0
    a_idm::Float64 = 0.0
    other_acc::Float64 = 0.0
    s_des::Float64 = idm.s_min
    dist_at_merge::Float64 = 0.0
    ego_ttm::Float64 = 0.0
    veh_ttm::Float64 = 0.0
    front_car::Bool = false
    consider_merge::Bool = false
end

Base.rand(model::CooperativeIDM) = LaneFollowingAccel(model.a)
function AutomotiveDrivingModels.reset_hidden_state!(model::CooperativeIDM)
    reset_hidden_state!(model.idm)
    model.a = 0.0
    model.a_merge = 0.0
    model.a_idm = 0.0
    model.other_acc = 0.0
    model.s_des = model.idm.s_min
    model.dist_at_merge = 0.0
end

function AutomotiveDrivingModels.set_desired_speed!(model::CooperativeIDM, vdes::Float64)
    set_desired_speed!(model.idm, vdes)
end

function AutomotiveDrivingModels.observe!(model::CooperativeIDM, scene::Scene, roadway::Roadway, egoid::Int64)
    ego_ind = findfirst(egoid, scene)
    # @printf("OBSERVE COOPERATIVE IDM \n")
    ego = scene[ego_ind]
    fore = get_neighbor_fore_along_lane(scene, ego_ind, roadway, VehicleTargetPointFront(), VehicleTargetPointRear(), VehicleTargetPointFront())
    if fore.ind == nothing # uses the first vehicle on the lain as neighbor

        vehmin, vehind = findfirst_lane(scene, main_lane(model.env))
        headway = get_end(main_lane(model.env)) - ego.state.posF.s + vehmin.state.posF.s
        track_longitudinal!(model.idm, ego.state.v, vehmin.state.v, headway)
        # @printf("No neighbor, ID: %d | neighbor %d | headway %2.1f \n", egoid, scene[vehind].id, headway)
    else
        # @printf("ID: %d | neighbor %d \n", egoid, scene[fore.ind].id)
        observe!(model.idm, scene, roadway, egoid)
    end
    a_idm = model.idm.a
    model.a_idm = a_idm
    veh = find_merge_vehicle(model.env, scene)
    if veh == nothing || veh.state.posF.s < model.fov
        # println("No merge vehicle")
        model.a = model.a_idm
    else
        model.other_acc = 0.0
        model.a = 0.0
        model.a
        ego_ttm = time_to_merge(model.env, ego, model.a)
        veh_ttm = time_to_merge(model.env, veh, model.other_acc)
        model.ego_ttm = ego_ttm
        model.veh_ttm = veh_ttm
        if ( ego_ttm < 0.0 || ego_ttm < veh_ttm || veh_ttm == Inf)
            # println("Ego TTM < Merge TTM, ignoring")
            ego_ttm < veh_ttm
            model.a = model.a_idm
            model.consider_merge = false
            model.front_car = false
        else
            model.consider_merge = true
            if veh_ttm < model.c*ego_ttm
                model.front_car = true
                headway = distance_projection(model.env, veh) - distance_projection(model.env, ego)
                headway -= veh.def.length
                v_oth = veh.state.v
                v_ego = ego.state.v
                # @show "tracking front car"
                track_longitudinal!(model.idm, v_ego, v_oth, headway)
                model.a_merge = model.idm.a
                model.a = min(model.a_merge, model.a_idm)
            else
                model.a = model.a_idm
                model.front_car = false
            end
        end
    end
    return model
end

"""
    findfirst_lane(scene::Scene, lane::Lane)
find the first vehicle on the lane (in terms of longitudinal position)
"""
function findfirst_lane(scene::Scene, lane::Lane)
    s_min = Inf
    vehmin = nothing
    vehind = nothing
    for (i, veh) in enumerate(scene)
        if veh.state.posF.roadind.tag == lane.tag && veh.state.posF.s < s_min
            s_min = veh.state.posF.s
            vehmin = veh
            vehind = i
        end
    end
    return vehmin, vehind
end

function get_end(lane::Lane)
    return lane.curve[end].s
end


text:
"""
    Genser.serialize([io::IO,] value, M::Module; kwargs...)

Serialize a value with an explicit serialization module.

If `io` is not provided, a `Vector{UInt8}` of the contents is returned.

The keyword arguments are passed to the serialization module.
"""
function serialize(io::IO, value, M::Module; kwargs...)
    if !isdefined(M, :serialize_genser)
        throw(ArgumentError("$M cannot serialize data"))
    end
    genserdata = togenser(value)

    M.serialize_genser(io, genserdata; kwargs...)
end

function serialize(io::IO, value, ::Type{M}; kwargs...) where M <: MIME
    Mod = genser_type_for_mime(M)
    serialize(io, value, Mod::Module; kwargs...)
end

"""
    Genser.serialize([io::IO,] value, mime::MIME; kwargs...)

Serialize a value using the registered serialization module for a mime type.
"""
@inline function serialize(io::IO, value, m::MIME; kwargs...)
    serialize(io, value, typeof(m); kwargs...)
end

"""
    Genser.serialize([io::IO,] value, mime::AbstractString; kwargs...)

Serialize a value using the registered serialization module for a mime type provided as a string.
"""
@inline function serialize(io::IO, value, m::AbstractString; kwargs...)
    serialize(io, value, MIME(m); kwargs...)
end

@inline function serialize(value, m; kwargs...)
    io = IOBuffer()
    serialize(io, value, m; kwargs...)
    take!(io)
end


text:
module PassiveTracerFlows

include("traceradvdiff.jl")

end # module


text:
#=

  Broken weights problem in Julia ConstraintSolver.jl

  From
  http://www.mathlesstraveled.com/?p=701
  """
  Here's a fantastic problem I recently heard. Apparently it was first
  posed by Claude Gaspard Bachet de Méziriac in a book of arithmetic problems
  published in 1612, and can also be found in Heinrich Dorrie’s 100
  Great Problems of Elementary Mathematics.

      A merchant had a forty pound measuring weight that broke
      into four pieces as the result of a fall. When the pieces were
      subsequently weighed, it was found that the weight of each piece
      was a whole number of pounds and that the four pieces could be
      used to weigh every integral weight between 1 and 40 pounds. What
      were the weights of the pieces?

  Note that since this was a 17th-century merchant, he of course used a
  balance scale to weigh things. So, for example, he could use a 1-pound
  weight and a 4-pound weight to weigh a 3-pound object, by placing the
  3-pound object and 1-pound weight on one side of the scale, and
  the 4-pound weight on the other side.
  """

  Model created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#

using ConstraintSolver, JuMP
using Cbc, GLPK, Ipopt
using Printf
const CS = ConstraintSolver
include("constraints_utils.jl")

function broken_weights(n=4,m=40,print_solutions=true,all_solutions=true)

    cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
    glpk_optimizer = optimizer_with_attributes(GLPK.Optimizer)
    ipopt_optimizer = optimizer_with_attributes(Ipopt.Optimizer)

    model = Model(optimizer_with_attributes(CS.Optimizer,   # "all_solutions"=> all_solutions,
                                                            "all_optimal_solutions"=>all_solutions,
                                                            "logging"=>[],

                                                            "traverse_strategy"=>:BFS,
                                                            # "traverse_strategy"=>:DFS,
                                                            # "traverse_strategy"=>:DBFS,

                                                            # "branch_split"=>:Smallest,
                                                            # "branch_split"=>:Biggest,
                                                            "branch_split"=>:InHalf,

                                                            # https://wikunia.github.io/ConstraintSolver.jl/stable/options/#branch_strategy-(:Auto)
                                                            "branch_strategy" => :IMPS, # default
                                                            # "branch_strategy" => :ABS, # Activity Based Search
                                                            # "activity.decay" => 0.999, # default 0.999
                                                            # "activity.max_probes" => 10, # default, 10
                                                            # "activity.max_confidence_deviation" => 20, # default 20

                                                            # "simplify"=>false,
                                                            # "simplify"=>true, # default

                                                            "time_limit"=>6,

                                                            # "backtrack" => false, # default true
                                                            # "backtrack_sorting" => false, # default true

                                                            # "lp_optimizer" => cbc_optimizer,
                                                            # "lp_optimizer" => glpk_optimizer,
                                                            # "lp_optimizer" => ipopt_optimizer,
                                        ))

    @variable(model, 1 <= weights[1:n] <= m, Int)
    @variable(model, -1 <= x[1:m,1:n] <= 1, Int)

    @constraint(model, weights in CS.AllDifferent())
    @constraint(model, m == sum(weights))

    increasing(model,weights)
    @constraint(model, m == sum(weights))

    table = resize_matrix([ [i,j,i * j] for i in -m:m, j in -m:m])
    for j in 1:m
        # scalar_product(model,x[j,:],weights,j) # This is a non-linear constraint!
        # This reformulation works.
        t = @variable(model, [1:n], CS.Integers(-m:m))
        for i in 1:n
            @constraint(model, [x[j,i], weights[i], t[i]] in CS.TableSet(table))
        end
        @constraint(model, j == sum(t))
    end

    @objective(model,Min,weights[n])

    # Solve the problem
    optimize!(model)

    status = JuMP.termination_status(model)
    # println("status:$status")
    if status == MOI.OPTIMAL
        num_sols = MOI.get(model, MOI.ResultCount())
        println("num_sols:$num_sols\n")
        if print_solutions
            for sol in 1:num_sols
                println("solution #$sol")
                weights_val = convert.(Integer,JuMP.value.(weights; result=sol))
                x_val = convert.(Integer,JuMP.value.(x; result=sol))
                println("weights:$weights_val")
                # display(x_val)
                for j in 1:m
                    @printf "%2d: " j
                    for i in 1:n
                        @printf "%2d" x_val[j,i]
                    end
                    println()
                end
            end
        end
    else
        println("status:$status")
    end

    return status
end

@time broken_weights(4,40)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Desription
# ==============================================================================
#
#   Tests related to the operations with Euler angles.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# File: ./src/angle.jl
# ====================

# Functions: *
# ------------

@testset "Operations with Euler angles: *"  begin
    # We do not need comprehensive test here. `*` first converts the Euler
    # angles to quaternions, and then it performs the multiplication. At the
    # end, it converts the quaternion back to Euler angles.

    for T in (Float32, Float64)
        ea1 = EulerAngles(deg2rad(T(45)), 0, 0, :ZYX)
        ea2 = EulerAngles(0, deg2rad(T(15)), 0, :YZY)
        ea3 = ea2 * ea1
        @test eltype(ea3) === T
        @test ea3.a1 ≈ T(0)
        @test ea3.a2 ≈ deg2rad(T(60))
        @test ea3.a3 ≈ T(0)
        @test ea3.rot_seq == :YZY
    end
end


text:
#module gigaword_64k

using Revise, Underscores
using FileTrees, Glob, DataStructures
using Gumbo, WordTokenizers
using AbstractTrees, Test
using JDF, DataFrames
using IterTools, StatsBase, Lazy
import Lazy: list, tail, first

#export read_and_wc, process_part_of_tree

# cd /users/yh31/scratch/projects/gigaword_64k
# /users/yh31/scratch/datasets/entity_linking/raw_data/gigaword/giga5/data

# [yh31@node1323 data]$ du -shc ./*
# 4.9G    ./afp_eng
# 7.8G    ./apw_eng
# 263M    ./cna_eng
# 1.7G    ./ltw_eng
# 8.8G    ./nyt_eng
# 112M    ./wpb_eng
# 2.5G    ./xin_eng
# 26G     total

# there's about 12 files per year in afp_eng

# ==========
# CONSTANTS
# ==========

# =====
# UTILS
# =====

showall(x) = show(stdout, "text/plain", x)

is_wanted(s::String) = length(s) >= 7 || all(isletter.(collect(s)))
# A reason to be more stringent here: would allow us to have smaller dictionaries!
# imo it's OK to be stringent about how much to filter at this stage, because ultimately what we're doing is filtering out stuff from other copora based on what's in the exclusion list that were making. Being stringent here just means that we'd have a smaller exclusion list than we'd otherwise have. And we can always apply more filters later.

# naive implementation of most common things
most_common(c::Accumulator) = most_common(c, length(c))
most_common(c::Accumulator, k) = sort(collect(c), by=kv->kv[2], rev=true)[1:k]

"Load JDF and get it into DataFrame form"
jdf_to_df(jdf) = JDF.load(jdf) |> DataFrame

"""
String, String, String -> array

Returns list of dirs that are in `base` or any of its subdirectories
"""
function get_dir_paths(base::String, extension::String)
    dirlist = []

    for (root, dirs, files) in walkdir(base)
        for dir in dirs
            if endswith(dir, "jdf")
                push!(dirlist, joinpath(root, dir))
            end
        end
    end

    return dirlist
end

# need foldl because of how merge works

@rec lzfoldl(f::Function, v, xs::List) = begin
    isempty(xs) ? v : lzfoldl(f, f(v, first(xs)), tail(xs))
end

lzfoldl(f::Function, xs::List) = begin
  isempty(xs) ? f() : lzfoldl(f, first(xs), tail(xs))
end

# ==================================
# MAIN functions (including helpers)
# ==================================

"Returns a df that consists of pairs in `acc`"
function acc_to_df(acc)
    # don't restrict this to accumulators,
    # because the result of nlargest on accumulators is not an accumulator
    df = DataFrame([String, Int], [:word, :freq])

    for (k, v) in acc
        push!(df, (k, v))
    end
    return df
end

"Converts a df (with same format as those produced by `acc_to_df`) to acc"
function df_to_acc(df :: DataFrame)
    acc = counter(String)
    map(eachrow(df)) do row acc[row.word] = row.freq end

    return acc
end

"""
Tokenizes doc and returns word counts for all words ∈ doc
"""
function count_words_from_file(doc :: HTMLDocument)
    doc_counter = counter(String)

    "Helper function"
    function inc_counter_with_doc()
        for elt in PreOrderDFS(doc.root)
            if isa(elt, HTMLText)
                #=
                1. tokenize the text
                    (note: each use of tokenize returns a 1-dim array of strings)
                2. increment counter with lowercase repn of each word
                =#
                @_ tokenize(elt.text) |>
                filter(is_wanted, __) |>
                foldl((lst, str) -> inc!(doc_counter, lowercase(str)), __; init = [])
            end
        end
    end

    inc_counter_with_doc()
    return doc_counter
end

"Composes file reading and count_words_from_file"
read_and_wc(filepath::String) = read(filepath, String) |>
                 parsehtml |>
                 count_words_from_file
# TO DO: might not need underscores.jl here!

"Returns year from filename in gigaword corpus"
year_from_fnm(fnm :: String) = split(fnm, "_")[3][1:4]

"Picks min(`k`, size of array) elements from array"
pick_elts(arr, k :: Integer) = length(arr) <= k ? arr : sample(arr, k, replace=false)

"Given FileTree for dir, return FileTree with __only__ the files we want to sample from"
function process_part_of_tree(path_of_tree :: String, path_output :: String, n_items :: Integer)
    # `n_items` is number of items we want to get per year

    data_tree = FileTree(path_of_tree)

    # 1. Get the list of names of files in tree
    filenames = name.(files(data_tree))

    # 2. For each year that's represented, randomly choose min(array_size, n_items) of the files in that year

    pick_elts_from_arr(arr) = pick_elts(arr, n_items)

    chosen_files = @_ filenames |>
            collect(IterTools.groupby(year_from_fnm, __)) |>
            map(pick_elts_from_arr, __) |>
            vcat(__...) |> Set(__)

    # 3. Get the resulting tree with `filter`
    filtered_tree = filter(x->x.name ∈ chosen_files, data_tree, dirs=false)

    # 4. load it; count words; save it
    loaded_tree = FileTrees.load(filtered_tree; lazy = true) do file
                       @_ file |>
                       string(FileTrees.path(__)) |>
                       read_and_wc(__) |>
                       acc_to_df(__)
                  end

    out_tree = FileTrees.rename(loaded_tree, path_output)

    FileTrees.save(out_tree) do acc_df
        savejdf(string(FileTrees.path(acc_df)) * ".jdf", FileTrees.get(acc_df))
    #       @warn "failed to save $(string(FileTrees.path(acc_df)))"
    #   end
    end
end

"Reduces jdfs in `base` (or some subdir thereof) to a final acc"
function reduce_jdfs(basepath::String)
    jdf_paths = @lazy get_dir_paths(basepath, ".jdf")
    final_acc = @_ lazymap(df_to_acc(jdf_to_df(_)), jdf_paths) |>
                lzfoldl(merge!(_1, _2) , __)
                # need to use foldl if using merge! and not merge
    return final_acc
end

```
Reduces all the JDFs and saves two final JDFs,
one with the top 64k words and one with the top 70k words
```
function reduce_jdfs_and_save(input_path::String, output_path::String; ns_tuple = (70_000, 64_000))
    final_acc = reduce_jdfs(input_path)
    final_acc_len = length(final_acc)

    n1 = min(ns_tuple[1], final_acc_len); n2 = min(ns_tuple[2], final_acc_len)

    final_n1_acc = @_ final_acc |> nlargest(__, n1)
    final_df_larger = final_n1_acc |> acc_to_df
    savejdf(joinpath(output_path, "final_df_$n1.jdf"), final_df_larger)

    final_df_smaller = @_ final_acc |>
                        nlargest(__, n2) |>
                        acc_to_df
    savejdf(joinpath(output_path, "final_df_$n2.jdf"), final_df_smaller)
end

#=

function reduce_files()

end

counter_dfs = FileTree(path_of_tree)

# trying to understand the tutorial...
# http://shashi.biz/FileTrees.jl/tree-manipulation/
testtree = maketree("dir"=>
                [string(i)=>[(name=string(j), value=(i,j)) for j in 1:5] for i=1:5])

vcated = mapsubtrees(testtree, glob"*") do subtree
    reducevalues(vcat, subtree)
end

=#

#end

#= Notes for future

# Re FileTrees.jl

## Not sure quid difference between `filtered_tree` and `data_tree[filtered_tree]`
typeof(filtered_tree) # FileTree
typeof(data_tree[filtered_tree]) # FileTree

filtered_tree == data_tree[filtered_tree]
# false, hmm, not sure why
files(filtered_tree) == files(data_tree[filtered_tree])
# this is true

# for testing
filtered_filenames = name.(files(filtered_tree))
Set(filtered_filenames) == chosen_files

* EzXML won't work with SGML.
* Transducers.jl didn't seem to be any faster than base foldl and filter, though that might change in the future

using Transducers

function test_no_transducer(doc)
    for elt in PreOrderDFS(doc.root)
        if isa(elt, HTMLText)

            1. tokenize the text
                (note: each use of tokenize returns a 1-dim array of strings)
            2. filter out punctuation, numbers, etc
            3. increment counter accordingly

            @_ tokenize(elt.text) |>
            filter(is_wanted, __) |>
            foldl((lst, str) -> inc!(doc_counter, str), __; init = [])
        end
    end
end

function test_transducer(doc)
    for elt in PreOrderDFS(doc.root)
        if isa(elt, HTMLText)
            tokenize(elt.text) |>
            Filter(is_wanted) |>
            foldxl((lst, str) -> inc!(doc_counter, str); init = [])
            # to do: figure out why we don't need underscores for this!
        end
    end
end

@time test_no_transducer(doc)
# 200308: 9.063 s
# 200304: 82.035s
# reduce version: 82.059720 seconds (211.31 M allocations: 13.842 GiB, 4.55% gc time)
# re-doing it again with foldl
#  77.735108 seconds (212.57 M allocations: 13.857 GiB, 4.71% gc time)

@time test_transducer(doc)
# 200308: 9.11s!
# 200304: 75.447160 seconds (211.13 M allocations: 13.748 GiB, 3.60% gc time)
=#