text:
# Copyright (c) Microsoft Corporation.
# Licensed under the MIT license.

module BaTests

dir = @__DIR__
ba_test_implementations = Tuple{String, Float64}[
    ("$dir/../../../../src/julia/modules/Zygote/ZygoteBA.jl", 1e-8)
]

map!(tup -> (abspath(tup[1]), tup[2]), ba_test_implementations, ba_test_implementations)

using Test
include("../../../../src/julia/shared/load.jl")
include("../../../../src/julia/runner/load.jl")
include("load.jl")
import ADPerfTest
import TestLoader
using BAData
using TestUtils

@testset "BA Module Tests" begin

@testset "BA Module Test ($(basename(module_path)))" for (module_path, tolerance) in ba_test_implementations
    module_dir, module_filename = splitdir(module_path)
    module_name, module_ext = splitext(module_filename)
    # Assert that the path to the module is correct
    ext_r = @test module_ext == ".jl"
    ispath_r = @test (ispath(module_path) && !isdir(module_path))
    if isa(ext_r, Test.Pass) && isa(ispath_r, Test.Pass)
        need_modify_load_path = !(module_dir ∈ LOAD_PATH)
        if need_modify_load_path
            push!(LOAD_PATH, module_dir)
        end
        try
            test = TestLoader.get_ba_test(module_name)
            # Module loads
            @test isa(test, ADPerfTest.Test{BAInput, BAOutput})
            input = load_ba_input("$dir/../../../../data/ba/test.txt")
            test.prepare!(test.context, input)
            test.calculate_objective!(test.context, 1)
            output = empty_ba_output()
            test.output!(output, test.context)
            correct_reproj_err_elem = [ -2.69048849235189402e-01, 2.59944792677901881e-01 ]
            correct_w_err_elem = 8.26092651515999976e-01
            # Objective is calculated correctly
            @test (2, 10) == size(output.reproj_err)
            @test (10,) == size(output.w_err)
            for i ∈ 1:10
               @test correct_reproj_err_elem[1] ≈ output.reproj_err[1, i] atol=tolerance
               @test correct_reproj_err_elem[2] ≈ output.reproj_err[2, i] atol=tolerance
            end
            for i ∈ 1:10
               @test correct_w_err_elem ≈ output.w_err[i] atol=tolerance
            end
            test.calculate_objective!(test.context, 3)
            test.output!(output, test.context)
            # Objective is calculated correctly for time == 3
            @test (2, 10) == size(output.reproj_err)
            @test (10,) == size(output.w_err)
            for i ∈ 1:10
               @test correct_reproj_err_elem[1] ≈ output.reproj_err[1, i] atol=tolerance
               @test correct_reproj_err_elem[2] ≈ output.reproj_err[2, i] atol=tolerance
            end
            for i ∈ 1:10
               @test correct_w_err_elem ≈ output.w_err[i] atol=tolerance
            end

            correct_jacobian_nrows = 30
            correct_jacobian_ncols = 62
            correct_jacobian_rows = [ 0, 15, 30, 45, 60, 75, 90, 105, 120, 135, 150, 165, 180, 195, 210, 225, 240, 255, 270, 285, 300, 301, 302, 303, 304, 305, 306, 307, 308, 309, 310 ]
            correct_jacobian_cols = [ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 22, 23, 24, 52, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 22, 23, 24, 52, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 25, 26, 27, 53, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 25, 26, 27, 53, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 28, 29, 30, 54, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 28, 29, 30, 54, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 31, 32, 33, 55, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 31, 32, 33, 55, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 34, 35, 36, 56, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 34, 35, 36, 56, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 37, 38, 39, 57, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 37, 38, 39, 57, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 40, 41, 42, 58, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 40, 41, 42, 58, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 43, 44, 45, 59, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 43, 44, 45, 59, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 46, 47, 48, 60, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 46, 47, 48, 60, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 49, 50, 51, 61, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 49, 50, 51, 61, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61 ]
            correct_jacobian_vals = [ 2.28877202208246757e+02, 6.34574811495545418e+02, -7.82222866259340549e+02, 2.42892615607159668e+00, -1.17828079628011313e+01, 2.54169312487743460e+00, -1.03657084958518086e+00, 4.17022000000000004e-01, 0.00000000000000000e+00, -3.50739521096005205e+02, -9.12107773668008576e+02, -2.42892615607159668e+00, 1.17828079628011313e+01, -2.54169312487743460e+00, -6.45167039712987389e-01, -1.20542435994996879e+02, -3.85673240766460424e+02, 9.75476291403326456e+01, -1.78372108529576567e+00, 4.15466799433126077e+00, 2.04025718029898906e+00, 3.49176397433145880e-01, 0.00000000000000000e+00, 4.17022000000000004e-01, 1.18149147704414503e+02, 3.07250108960343255e+02, 1.78372108529576567e+00, -4.15466799433126077e+00, -2.04025718029898906e+00, 6.23335921553064054e-01, 2.28877202208246757e+02, 6.34574811495545418e+02, -7.82222866259340549e+02, 2.42892615607159668e+00, -1.17828079628011313e+01, 2.54169312487743460e+00, -1.03657084958518086e+00, 4.17022000000000004e-01, 0.00000000000000000e+00, -3.50739521096005205e+02, -9.12107773668008576e+02, -2.42892615607159668e+00, 1.17828079628011313e+01, -2.54169312487743460e+00, -6.45167039712987389e-01, -1.20542435994996879e+02, -3.85673240766460424e+02, 9.75476291403326456e+01, -1.78372108529576567e+00, 4.15466799433126077e+00, 2.04025718029898906e+00, 3.49176397433145880e-01, 0.00000000000000000e+00, 4.17022000000000004e-01, 1.18149147704414503e+02, 3.07250108960343255e+02, 1.78372108529576567e+00, -4.15466799433126077e+00, -2.04025718029898906e+00, 6.23335921553064054e-01, 2.28877202208246757e+02, 6.34574811495545418e+02, -7.82222866259340549e+02, 2.42892615607159668e+00, -1.17828079628011313e+01, 2.54169312487743460e+00, -1.03657084958518086e+00, 4.17022000000000004e-01, 0.00000000000000000e+00, -3.50739521096005205e+02, -9.12107773668008576e+02, -2.42892615607159668e+00, 1.17828079628011313e+01, -2.54169312487743460e+00, -6.45167039712987389e-01, -1.20542435994996879e+02, -3.85673240766460424e+02, 9.75476291403326456e+01, -1.78372108529576567e+00, 4.15466799433126077e+00, 2.04025718029898906e+00, 3.49176397433145880e-01, 0.00000000000000000e+00, 4.17022000000000004e-01, 1.18149147704414503e+02, 3.07250108960343255e+02, 1.78372108529576567e+00, -4.15466799433126077e+00, -2.04025718029898906e+00, 6.23335921553064054e-01, 2.28877202208246757e+02, 6.34574811495545418e+02, -7.82222866259340549e+02, 2.42892615607159668e+00, -1.17828079628011313e+01, 2.54169312487743460e+00, -1.03657084958518086e+00, 4.17022000000000004e-01, 0.00000000000000000e+00, -3.50739521096005205e+02, -9.12107773668008576e+02, -2.42892615607159668e+00, 1.17828079628011313e+01, -2.54169312487743460e+00, -6.45167039712987389e-01, -1.20542435994996879e+02, -3.85673240766460424e+02, 9.75476291403326456e+01, -1.78372108529576567e+00, 4.15466799433126077e+00, 2.04025718029898906e+00, 3.49176397433145880e-01, 0.00000000000000000e+00, 4.17022000000000004e-01, 1.18149147704414503e+02, 3.07250108960343255e+02, 1.78372108529576567e+00, -4.15466799433126077e+00, -2.04025718029898906e+00, 6.23335921553064054e-01, 2.28877202208246757e+02, 6.34574811495545418e+02, -7.82222866259340549e+02, 2.42892615607159668e+00, -1.17828079628011313e+01, 2.54169312487743460e+00, -1.03657084958518086e+00, 4.17022000000000004e-01, 0.00000000000000000e+00, -3.50739521096005205e+02, -9.12107773668008576e+02, -2.42892615607159668e+00, 1.17828079628011313e+01, -2.54169312487743460e+00, -6.45167039712987389e-01, -1.20542435994996879e+02, -3.85673240766460424e+02, 9.75476291403326456e+01, -1.78372108529576567e+00, 4.15466799433126077e+00, 2.04025718029898906e+00, 3.49176397433145880e-01, 0.00000000000000000e+00, 4.17022000000000004e-01, 1.18149147704414503e+02, 3.07250108960343255e+02, 1.78372108529576567e+00, -4.15466799433126077e+00, -2.04025718029898906e+00, 6.23335921553064054e-01, 2.28877202208246757e+02, 6.34574811495545418e+02, -7.82222866259340549e+02, 2.42892615607159668e+00, -1.17828079628011313e+01, 2.54169312487743460e+00, -1.03657084958518086e+00, 4.17022000000000004e-01, 0.00000000000000000e+00, -3.50739521096005205e+02, -9.12107773668008576e+02, -2.42892615607159668e+00, 1.17828079628011313e+01, -2.54169312487743460e+00, -6.45167039712987389e-01, -1.20542435994996879e+02, -3.85673240766460424e+02, 9.75476291403326456e+01, -1.78372108529576567e+00, 4.15466799433126077e+00, 2.04025718029898906e+00, 3.49176397433145880e-01, 0.00000000000000000e+00, 4.17022000000000004e-01, 1.18149147704414503e+02, 3.07250108960343255e+02, 1.78372108529576567e+00, -4.15466799433126077e+00, -2.04025718029898906e+00, 6.23335921553064054e-01, 2.28877202208246757e+02, 6.34574811495545418e+02, -7.82222866259340549e+02, 2.42892615607159668e+00, -1.17828079628011313e+01, 2.54169312487743460e+00, -1.03657084958518086e+00, 4.17022000000000004e-01, 0.00000000000000000e+00, -3.50739521096005205e+02, -9.12107773668008576e+02, -2.42892615607159668e+00, 1.17828079628011313e+01, -2.54169312487743460e+00, -6.45167039712987389e-01, -1.20542435994996879e+02, -3.85673240766460424e+02, 9.75476291403326456e+01, -1.78372108529576567e+00, 4.15466799433126077e+00, 2.04025718029898906e+00, 3.49176397433145880e-01, 0.00000000000000000e+00, 4.17022000000000004e-01, 1.18149147704414503e+02, 3.07250108960343255e+02, 1.78372108529576567e+00, -4.15466799433126077e+00, -2.04025718029898906e+00, 6.23335921553064054e-01, 2.28877202208246757e+02, 6.34574811495545418e+02, -7.82222866259340549e+02, 2.42892615607159668e+00, -1.17828079628011313e+01, 2.54169312487743460e+00, -1.03657084958518086e+00, 4.17022000000000004e-01, 0.00000000000000000e+00, -3.50739521096005205e+02, -9.12107773668008576e+02, -2.42892615607159668e+00, 1.17828079628011313e+01, -2.54169312487743460e+00, -6.45167039712987389e-01, -1.20542435994996879e+02, -3.85673240766460424e+02, 9.75476291403326456e+01, -1.78372108529576567e+00, 4.15466799433126077e+00, 2.04025718029898906e+00, 3.49176397433145880e-01, 0.00000000000000000e+00, 4.17022000000000004e-01, 1.18149147704414503e+02, 3.07250108960343255e+02, 1.78372108529576567e+00, -4.15466799433126077e+00, -2.04025718029898906e+00, 6.23335921553064054e-01, 2.28877202208246757e+02, 6.34574811495545418e+02, -7.82222866259340549e+02, 2.42892615607159668e+00, -1.17828079628011313e+01, 2.54169312487743460e+00, -1.03657084958518086e+00, 4.17022000000000004e-01, 0.00000000000000000e+00, -3.50739521096005205e+02, -9.12107773668008576e+02, -2.42892615607159668e+00, 1.17828079628011313e+01, -2.54169312487743460e+00, -6.45167039712987389e-01, -1.20542435994996879e+02, -3.85673240766460424e+02, 9.75476291403326456e+01, -1.78372108529576567e+00, 4.15466799433126077e+00, 2.04025718029898906e+00, 3.49176397433145880e-01, 0.00000000000000000e+00, 4.17022000000000004e-01, 1.18149147704414503e+02, 3.07250108960343255e+02, 1.78372108529576567e+00, -4.15466799433126077e+00, -2.04025718029898906e+00, 6.23335921553064054e-01, 2.28877202208246757e+02, 6.34574811495545418e+02, -7.82222866259340549e+02, 2.42892615607159668e+00, -1.17828079628011313e+01, 2.54169312487743460e+00, -1.03657084958518086e+00, 4.17022000000000004e-01, 0.00000000000000000e+00, -3.50739521096005205e+02, -9.12107773668008576e+02, -2.42892615607159668e+00, 1.17828079628011313e+01, -2.54169312487743460e+00, -6.45167039712987389e-01, -1.20542435994996879e+02, -3.85673240766460424e+02, 9.75476291403326456e+01, -1.78372108529576567e+00, 4.15466799433126077e+00, 2.04025718029898906e+00, 3.49176397433145880e-01, 0.00000000000000000e+00, 4.17022000000000004e-01, 1.18149147704414503e+02, 3.07250108960343255e+02, 1.78372108529576567e+00, -4.15466799433126077e+00, -2.04025718029898906e+00, 6.23335921553064054e-01, -8.34044000000000008e-01, -8.34044000000000008e-01, -8.34044000000000008e-01, -8.34044000000000008e-01, -8.34044000000000008e-01, -8.34044000000000008e-01, -8.34044000000000008e-01, -8.34044000000000008e-01, -8.34044000000000008e-01, -8.34044000000000008e-01 ]

            test.calculate_jacobian!(test.context, 1)
            test.output!(output, test.context)
            # Jacobian is calculated correctly
            @test correct_jacobian_nrows == output.jacobian.nrows
            @test correct_jacobian_ncols == output.jacobian.ncols
            @test correct_jacobian_rows == output.jacobian.rows
            @test correct_jacobian_cols == output.jacobian.cols
            @test correct_jacobian_vals ≈ output.jacobian.vals atol=tolerance

            test.calculate_jacobian!(test.context, 3)
            test.output!(output, test.context)
            # Jacobian is calculated correctly for time == 3
            @test correct_jacobian_nrows == output.jacobian.nrows
            @test correct_jacobian_ncols == output.jacobian.ncols
            @test correct_jacobian_rows == output.jacobian.rows
            @test correct_jacobian_cols == output.jacobian.cols
            @test correct_jacobian_vals ≈ output.jacobian.vals atol=tolerance

            # Objective runs multiple times
            @test can_objective_run_multiple_times!(test.context, test.calculate_objective!)
            # Jacobian runs multiple times
            @test can_objective_run_multiple_times!(test.context, test.calculate_jacobian!)
        finally
            pop!(LOAD_PATH)
        end
    end
end

end

end

max_stars_repo_name: dsyme/ADBench

max_issues_repo_name: dsyme/ADBench

max_forks_repo_name: dsyme/ADBench


text:
# A SVM solver based on the dual coordinate descent method.

# The naif implementation.
function svm_dualcoordinatedescent(X::AbstractMatrix,
                                   Y::AbstractVector,
                                   C::Real,
                                   ϵ::Real;
                                   maxiter::Integer=10000)
  α = C/2*ones(Y)
  w = X'*(Y.*α)
  Q = Diagonal(Y)*X*X'*Diagonal(Y)
  primal(w) = norm(w)^2/2 + C*sum(max.(0, ones(Y)-Y.*(X*w)))
  dual(α) = -α⋅(Q*α)/2 + sum(α)
  wmin = copy(w)
  wminvalue = primal(w)
  dualitygaps = [wminvalue-dual(α)]
  n = 0
  while dualitygaps[end] > ϵ && n < maxiter
    for i = 1:length(Y)
      ∇ᵢ = Q[i,:]⋅α - 1
      αᵢ = α[i]
      α[i] = min(max(α[i]-∇ᵢ/Q[i, i], 0), C)
      w += (α[i]-αᵢ)*Y[i]*X[i,:]
    end
    wvalue = primal(w)
    if wvalue < wminvalue
      wmin, wminvalue = copy(w), wvalue
    end
    push!(dualitygaps, wminvalue-dual(α))
    n += 1
  end
  wmin, α, dualitygaps
end

# Add random permutation of indices inside each outer iteration and a
# more efficient way to compute the gradient of each univariate function.
function svm_dualcoordinatedescentopt(X::AbstractMatrix,
                                      Y::AbstractVector,
                                      C::Integer,
                                      iter::Integer)
  α = C/2*ones(Y)
  w = X'*(Y.*α)
  Qdiag = sum(X.^2, 2)
  primal(w) = norm(w)^2/2 + C*sum(max.(0, ones(Y)-Y.*(X*w)))
  wmin = copy(w)
  wminvalue = primal(w)
  for k = 1:iter
    for i = randperm(length(Y))
      ∇ᵢ = Y[i]*w⋅X[i,:] - 1
      p∇ᵢ = α[i]==0?min(∇ᵢ, 0):α[i]==C?max(∇ᵢ, 0):∇ᵢ
      if p∇ᵢ ≠ 0
        αᵢ = α[i]
        α[i] = min(max(α[i]-∇ᵢ/Qdiag[i], 0), C)
        w += (α[i]-αᵢ)*Y[i]*X[i,:]
      end
    end
    wvalue = primal(w)
    if wvalue < wminvalue
      wmin, wminvalue = copy(w), wvalue
    end
  end
  wmin, α
end

ext: jl


text:
################################################################################
# Observations

@doc raw"""
    mutable struct Observation{T <: Real}

An *observation* is a sequence of time-ordered data (TOD) that has
been acquired by some detector. It implements the following fields:

- `time`: an array of `N` elements, representing the time of the
  sample. Being defined as an `AbstractArray`, a range (e.g.,
  `0.0:0.1:3600.0`) can be used to save memory
- `pixidx`: an array of `N` elements, each being a pixel index in a
  map. In no way the Healpix pixelization scheme is enforced; in fact,
  the map must not even be spherical.
- `psi_angle`: an array of `N` elements, each containing the
  orientation of the polarization angle with respect to some fixed
  reference system (e.g., the Ecliptic plane). The angles must be in
  **radians**.
- `tod`: an array of `N` elements, containing the actual data measured
  by the instrument and already calibrated to some physical units
  (e.g., K_CMB)
- `sigma_squared`: an array of `N` elements, containing the squared
  RMS of each sample. To save memory, you should usually use a
  `RunLengthArray` here.
- `flagged`: a Boolean array of `N` elements, telling whether the
  sample should be discarded (`true`) or not (`false`) during the
  map-making process.
- `name`: a string representing the detector. This is used only for
  debugging purposes.

To create an observation, the constructor takes *keyword arguments*
instead of parameters; in this way, the code should be more
readable. Also, sensible defaults will be provided for the missing
fields:

```julia
# We do not specify "sigma_squared" nor "flagged", so they will be
# initialized to a vector of ones and a vector of `false`
obs = Observation{Float64}(
    time = 0.0:0.1:3600.0,
    pixidx = 1:3601,
    psi_angle = zeros(3601),
    tod = randn(3601),
)
```

"""
mutable struct Observation{T <: Real}
    time::AbstractVector{T}
    pixidx::Vector{Int}
    psi_angle::Vector{T}

    tod::Vector{T}
    sigma_squared::AbstractVector{T}
    flagged::AbstractVector{Bool}

    name::String

    function Observation{T}(
        ;
        time = T[],
        pixidx = Int[],
        psi_angle = T[],
        tod = T[],
        sigma_squared = T[],
        flagged = Bool[],
        name = "unnamed",
    ) where {T <: Real}
        @assert !isempty(pixidx)

        nsamples = length(pixidx)

        int_time = isempty(time) ? (0:(nsamples - 1)) : time
        int_tod = isempty(tod) ? zeros(T, nsamples) : tod
        int_psi_angle = isempty(psi_angle) ? zeros(T, nsamples) : psi_angle
        int_sigma_squared = isempty(sigma_squared) ? RunLengthArray{Int,T}(nsamples, 1.0) : sigma_squared
        int_flagged = isempty(flagged) ? zeros(Bool, nsamples) : flagged

        @assert length(int_time) == nsamples
        @assert length(int_psi_angle) == nsamples
        @assert length(int_tod) == nsamples
        @assert length(int_sigma_squared) == nsamples
        @assert length(int_flagged) == nsamples

        new(
            int_time,
            pixidx,
            int_psi_angle,
            int_tod,
            int_sigma_squared,
            int_flagged,
            name,
        )
    end
end

function Base.show(io::IO, obs::Observation)
    print(io, "Observation(name => \"$(obs.name)\", samples = $(length(obs.time)))")
end

function Base.show(io::IO, ::MIME"text/plain", obs::Observation)

    isempty(obs.time) && println(io, "Observation for detector %(obs.name): empty")

    println(
        io,
        @sprintf(
            """Observation for detector %s:
Number of samples....................... %d
Start time.............................. %f
Duration................................ %f
Number of flagged samples............... %d (%.1f%%)
Average level of the TOD................ %e
Average noise^2 level................... %e
""",
            obs.name,
            length(obs.time),
            obs.time[1],
            obs.time[end] - obs.time[1],
            length(obs.flagged[obs.flagged]),
            100 * length(obs.flagged[obs.flagged]) / length(obs.flagged),
            mean(obs.tod),
            mean(obs.sigma_squared),
        )
    )
end


text:
using RobustDMD, Random, Printf

# Generate DMD Synthetic Data
#------------------------------------------------------------------------------
# dimensions
m = 500;    # temporal dimension
n = 1000;    # spatial dimension
k = 3;      # number of modes
T = Float64;
sigma = 1e-3;
# generate data
X, t, at, Bt = simDMD(m, n, k, T; seed=123456);
Xdat = X + sigma*Random.randn!(copy(X))
# specify the loss
function lossFunc(r)
	return sum(abs2, r)
end
#
function lossGrad(r)
	conj!(r)
end
#
# Create object
#------------------------------------------------------------------------------
params = DMDParams(k, Xdat, t, lossFunc, lossGrad);
#
# Initial Guess
#------------------------------------------------------------------------------
# create initial guess using exact dmd
a0, B0 = dmdexactestimate(m, n, k, Xdat, t);
#
# Create Object
#------------------------------------------------------------------------------
params = DMDParams(k, Xdat, t, lossFunc, lossGrad, inner_directl2=true);
copyto!(params.a, a0);
copyto!(params.B, B0);
#Random.randn!(params.ar);
#
# Apply Solver
#------------------------------------------------------------------------------
tau = 10;
eta = 1e1;
itm = 2000;
tol = 1e-6;
ptf = 10;

function prox_stab(ar)
    # project for numerical stability
    k = length(ar) >> 1;
    @inbounds for i = 1:k
        ar[i<<1-1] = min(2.0, ar[i<<1-1]);
    end
end

opts = DMD_SVRG_Options(tau, eta, itm=itm, tol=tol, ptf=ptf,
                        prox=prox_stab);

@time obj_his, err_his = solveDMD_withSVRG(params, opts);

@printf "%7.4e relative residual after SVRG\n" sqrt(dmdobjective(params)/lossFunc(X))

params2 = DMDParams(k,Xdat,t,lossFunc,lossGrad);
copyto!(params2.a,at)
@printf "%7.4e relative residual with true evals\n" sqrt(dmdobjective(params2)/lossFunc(X))
copyto!(params2.a,a0)
@printf "%7.4e relative residual with initial guess\n" sqrt(dmdobjective(params2)/lossFunc(X))

@printf "%7.4e abs err evals after SVRG\n" besterrperm(params.a,at)
@printf "%7.4e abs err evals initial guess\n" besterrperm(a0,at)


text:
include("Graph.jl")

using Random
using SparseArrays

function powerLaw(n; alp = 2.5, xmin = 1)
    Random.seed!(round(Int, time() * 10000))
    x = rand(n)
    for i = 1 : n
        x[i] = xmin * ((1 - x[i])^(-1.0/(alp - 1.0)))
    end
    xm = x[argmax(x)]
    x ./= xm
    return x
end

function Uniform(n)
    Random.seed!(round(Int, time() * 10000))
    x = rand(n)
    return x
end

function Exponential(n; lmd = 1, xmin = 1)
    Random.seed!(round(Int, time() * 10000))
    x = rand(n)
    for i = 1 : n
        x[i] = xmin - (1.0/lmd)*log(1-x[i])
    end
    xm = x[argmax(x)]
    x ./= xm
    return x
end

function getW(G)
	L = zeros(G.n, G.n)
	for (ID, u, v, w) in G.E
		L[u, u] += w
		L[v, v] += w
		L[u, v] -= w
		L[v, u] -= w
	end
	for i = 1 : G.n
		L[i, i] += 1.0
	end
	return inv(L)
end

function getL(G)
	L = zeros(G.n, G.n)
	for (ID, u, v, w) in G.E
		L[u, u] += w
		L[v, v] += w
		L[u, v] -= w
		L[v, u] -= w
	end
	return L
end

function getSparseIpL(G)
    d = ones(G.n)
    for (ID, u, v, w) in G.E
        d[u] += w
        d[v] += w
    end
    Is = zeros(Int32, G.m*2+G.n)
    Js = zeros(Int32, G.m*2+G.n)
    Vs = zeros(G.m*2+G.n)
    for (ID, u, v, w) in G.E
        Is[ID] = u
        Js[ID] = v
        Vs[ID] = -w
        Is[ID + G.m] = v
        Js[ID + G.m] = u
        Vs[ID + G.m] = -w
    end
    for i = 1 : G.n
        Is[G.m + G.m + i] = i
        Js[G.m + G.m + i] = i
        Vs[G.m + G.m + i] = d[i]
    end
    return sparse(Is, Js, Vs, G.n, G.n)
end

function getSparseL(G)
    d = zeros(G.n)
    for (ID, u, v, w) in G.E
        d[u] += w
        d[v] += w
    end
    Is = zeros(Int32, G.m*2+G.n)
    Js = zeros(Int32, G.m*2+G.n)
    Vs = zeros(G.m*2+G.n)
    for (ID, u, v, w) in G.E
        Is[ID] = u
        Js[ID] = v
        Vs[ID] = -w
        Is[ID + G.m] = v
        Js[ID + G.m] = u
        Vs[ID + G.m] = -w
    end
    for i = 1 : G.n
        Is[G.m + G.m + i] = i
        Js[G.m + G.m + i] = i
        Vs[G.m + G.m + i] = d[i]
    end
    return sparse(Is, Js, Vs, G.n, G.n)
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "ArcadeLearningEnvironment"
version = v"0.6.1"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/mgbellemare/Arcade-Learning-Environment/archive/v0.6.1.tar.gz", "8059a4087680da03878c1648a8ceb0413a341032ecaa44bef4ef1f9f829b6dde"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/Arcade-Learning-Environment-*/
atomic_patch -p1 ../patches/fix-dlext-macos.patch
atomic_patch -p1 ../patches/cmake-install-for-windows.patch
mkdir build && cd build
cmake -DCMAKE_INSTALL_PREFIX=$prefix \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DCMAKE_BUILD_TYPE=Release \
    -DUSE_SDL=OFF \
    -DBUILD_EXAMPLES=OFF \
    -DBUILD_CPP_LIB=OFF \
    -DBUILD_CLI=OFF \
    -DCMAKE_CXX_FLAGS="-I${includedir}" \
    -DCMAKE_SHARED_LINKER_FLAGS_INIT="-L${libdir}" \
    ..
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()
platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct("libale_c", :libale_c)
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
    Dependency(PackageSpec(name="Zlib_jll", uuid="83775a58-1f1d-513f-b197-d71354ab007a"))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
#########################################################################
##############################  SNAP  ###################################
#########################################################################
function force(c::Configuration, p::SNAP)
    A = get_snap(c, p)
    force = A[2:end-6, 1:(end-1)] * p.β[1:(end-1)]
    n = length(force)
    num_tuples = Int(n/3)
    force = [[force[3*(j-1)+k] for k = 1:3] for j = 1:num_tuples]
    return force
end

function force(r::Vector{Configuration}, p::SNAP)
    n = length(r)
    f = [[zeros(3) for j = 1:length(r[i].Atoms)] for i = 1:n]
    for i = 1:n
        f[i] = force(r[i], p)
    end
    return f
end


text:
using Dates
using Base.Threads

"""
    find_read_positions(files::Vector{<:Integer}, filebase::String,
                        blocks::AbstractVector{String}, parttype::Integer,
                        keylist::Vector{<:Integer}, key_info::Vector{InfoLine},
                        verbose::Bool)

Helper function to get positions and length of particle blocks in files.
"""
function find_read_positions(files::Vector{<:Integer}, filebase::String,
                             blocks::AbstractVector{String}, parttype::Integer,
                             keylist::Vector{<:Integer}, key_info::Vector{InfoLine},
                             verbose::Bool)

    # store number of file
    N_files = length(files)

    # allocate arrys to store reading information
    file_offset_key      = Vector{Vector{<:Integer}}(undef, N_files)
    file_part_per_key    = Vector{Vector{<:Integer}}(undef, N_files)
    file_block_positions = Vector{Dict{String,<:Integer}}(undef, N_files)

    @inbounds for i = 1:N_files

        filename = select_file(filebase, files[i])

        # read header of the file
        h = head_to_obj(filename)

        if h.npart[parttype+1] == 0
            error("No particles of type $parttype in file!")
        end

        filename_keyfile = filename * ".key"

        # Vector of field names (KEY, NKEY, OKEY)
        fields = getfield.(key_info, :block_name)

        # Reads the following in, respectively:
        # - key file data
        # - number of particles associated with PH key
        # - offsets in the blocks to get to the relevant particles
        keys_in_file, part_per_key, offset_key = read_key_block.(filename_keyfile,
                                                                 (key_info,),
                                                                 (fields,),
                                                                 ["KEY", "NKEY", "OKEY"],
                                                                 parttype)

        if verbose
            @info "Calculating index list..."
            t1 = Dates.now()
        end

        index_list = get_index_list(keylist, keys_in_file)

        if verbose
            t2 = Dates.now()
            @info "Index list done. Took: $(t2 - t1)"
            @info "Reading $(length(index_list)) key segments..."
        end

        # sort the offset arrays for simplification step
        sorted_offset = sortperm(offset_key[index_list])

        # save sorted arrays
        offset_key   = @views offset_key[index_list[sorted_offset]]
        part_per_key = @views part_per_key[index_list[sorted_offset]]

        # check if blocks can be joined
        use_block, part_per_key = join_blocks(offset_key, part_per_key)

        if verbose
            @info "Reduced independent blocks from $(length(offset_key)) to $(count(use_block))"
        end

        # store the arrays for later reading
        file_offset_key[i]   = offset_key[use_block]
        file_part_per_key[i] = part_per_key[use_block]

        file_block_positions[i] = get_block_positions(filename)

        # double-check if blocks are present
        blocks_in_file = String.(keys(file_block_positions[i]))
        for blockname in blocks
            if !block_present(filename, blockname, blocks_in_file)
                error("Block $blockname not present in file $filename !")
            end
        end

    end # for i = 1:N_files

    return file_offset_key, file_part_per_key, file_block_positions
end

"""
    read_key_block( filename_keyfile::String, key_info::Vector{InfoLine},
                    fields::Vector{String}, key::String, parttype::Integer)

Returns key block for different `key`. Valid values of `key`: KEY, NKEY, OKEY
"""
function read_key_block(filename_keyfile::String, key_info::Vector{InfoLine},
                        fields::Vector{String}, key::String, parttype::Integer)
    read_block(filename_keyfile, key,
               info = key_info[findfirst(==(key), fields)],
               parttype = parttype)
end

"""
    get_index_bounds(ids::Vector{<:Integer}, low_bounds::Vector{<:Integer}, high_bounds::Vector{<:Integer})

Returns sorted `Vector` of indices `i`, for which `low_bounds[i] ≤ ids[j] ≤ high_bounds[i]` for any `j`.
All parameters `ids`, `low_bounds`, and `high_bounds` have to already be sorted.
"""
function get_index_bounds(ids::Vector{<:Integer}, low_bounds::Vector{<:Integer}, high_bounds::Vector{<:Integer})

    nids    = length(ids)
    nbounds = length(low_bounds)
    ind_all = Vector{Int}(undef, min(nids, nbounds))

    icountall = 1
    icountids = 1
    @inbounds for icountbounds = 1:length(low_bounds)

        while icountids ≤ nids && ids[icountids] < low_bounds[icountbounds]
            icountids += 1
        end

        if icountids ≤ nids && ids[icountids] ≤ high_bounds[icountbounds]
            ind_all[icountall] = icountbounds
            icountall += 1
            icountids += 1
        end

    end

    return ind_all[1:(icountall-1)]
end

"""
    find_files_for_keys(filebase::String, nfiles::Integer, keylist::Vector{<:Integer})

Selects the files in which the particles associated with the given Peano-Hilbert keys are stored.
"""
function find_files_for_keys(filebase::String, nfiles::Integer, keylist::Vector{<:Integer})

    file_key_index = filebase * ".key.index"

    # if index file does not exist all key files need to be read
    if !isfile(file_key_index)
        return collect(0:nfiles-1)
    end

    # get the data from the index file (low and high lists are sorted)
    low_list, high_list, file_list = read_key_index(file_key_index)

    # get files that contain the wanted particles (parameters have to be sorted)
    index_bounds = get_index_bounds(sort(keylist), low_list, high_list)

    file_sort = file_list[index_bounds] |> unique! |> sort!

    return Int64.(file_sort)
end

"""
    get_index_list(list_to_find::Vector{<:Integer}, list_to_check::Vector{<:Integer})

Finds the indices at which `list_to_check` contains elements from `list_to_find`.
If both either of the lists are not sorted it uses a `Dict` lookup, otherwise it uses a `Vector` forward-search.
"""
function get_index_list(list_to_find::Vector{<:Integer}, list_to_check::Vector{<:Integer})

    # check if the lists are sorted
    if ( issorted(list_to_find) && issorted(list_to_check) )
        return get_index_list_arr(list_to_find, list_to_check)
    else
        return get_index_list_set(list_to_find, list_to_check)
    end
end

"""
    get_index_list_arr(list_to_find::Vector{<:Integer}, list_to_check::Vector{<:Integer})

Get positions in `list_to_check` where `list_to_check` matches `list_to_find`. Uses forward-searching in sorted array.
Both arrays have to be sorted.
"""
@inline function get_index_list_arr(list_to_find::Vector{<:Integer}, list_to_check::Vector{<:Integer})
    if isempty(list_to_find) || isempty(list_to_check)
        return Int64[]
    end

    narr1 = length(list_to_find)
    narr2 = length(list_to_check)

    ind_all   = zeros(Int64, narr1)

    icountall     = 1
    icountarr1    = 1
    icountarr2    = 1

    lend = false

    while !lend

        if list_to_check[icountarr2] == list_to_find[icountarr1]

            ind_all[icountall] = icountarr2
            icountall  += 1
            icountarr1 += 1
            icountarr2 += 1
        else  # list_to_check[icountnotarr2] == list_to_find[icountnotar1]
            if list_to_check[icountarr2] < list_to_find[icountarr1]
                icountarr2    += 1
            else # list_to_check[icountnotarr2] < list_to_find[icountnotar1]
                icountarr1    += 1
            end # list_to_check[icountnotarr2] < list_to_find[icountnotar1]
        end # list_to_check[icountnotarr2] == list_to_find[icountnotar1]
        if (icountarr2 > narr2 ) || (icountarr1 > narr1)
            lend = true
        end
    end # while !lend

    if icountall > 1
        return ind_all[1:icountall-1]
    else
        return Int64[]
    end
end

"""
    get_index_list_dict(list_to_find::Vector{<:Integer}, list_to_check::Vector{<:Integer})

Get positions in `list_to_check` where `list_to_check` matches `list_to_find`. Uses a `Dict` for lookup -> slower than the array search, but works on unsorted arrays.
"""
@inline function get_index_list_dict(list_to_find::Vector{<:Integer}, list_to_check::Vector{<:Integer})

    dict = Dict((n, i) for (i, n) in enumerate(list_to_check))
    result = Vector{Int}(undef, length(list_to_find))
    len = 0

    for k in list_to_find
        i = get(dict, k, nothing)
        if !isnothing(i)
            len += 1
            @inbounds result[len] = i
        end
    end
    return resize!(result, len)
end

"""
    get_index_list_set(list_to_find::Vector{<:Integer}, list_to_check::Vector{<:Integer})

Get positions in `list_to_check` where `list_to_check` matches `list_to_find`.
Uses a `Set` for lookup -> slower than the array search `get_index_list_arr`, but works on unsorted arrays
like `get_index_list_dict`, just faster.
"""
@inline function get_index_list_set(list_to_find::Vector{<:Integer}, list_to_check::Vector{<:Integer})

    set = Set(list_to_find)

    return findall(in.(list_to_check, (set,)))
end

"""
   join_blocks(offset_key, part_per_key)

Joins neighboring blocks to simplify read-in.
"""
@inline function join_blocks(offset_key::AbstractVector{<:Integer}, part_per_key::AbstractVector{<:Integer})

    use_block = trues(length(offset_key))

    icount = 1

    @inbounds for i = 2:length(offset_key)-1

        if offset_key[i] == offset_key[icount] + part_per_key[icount]
            part_per_key[icount] += part_per_key[i]
            use_block[i] = false
        else
            icount += 1
        end # if

    end # for

    return use_block, part_per_key
end


text: include("interleaved.jl")


text:
const _pkg_root = dirname(dirname(@__FILE__))
const _pkg_deps = joinpath(_pkg_root,"deps")
const _pkg_assets = joinpath(_pkg_root,"assets")

!isdir(_pkg_assets) && mkdir(_pkg_assets)

download("https://github.com/jgthms/bulma/releases/download/0.7.2/bulma-0.7.2.zip", joinpath(_pkg_assets, "bulma.zip"))


text:
# Spectrum calculator for AtmosGCM

struct AtmosGCMSpectraDiagnosticsParams <: DiagnosticsGroupParams
    nor::Float64
end

"""
    setup_atmos_spectra_diagnostics(
        ::AtmosGCMConfigType,
        interval::String,
        out_prefix::String;
        writer = NetCDFWriter(),
        interpol = nothing,
        nor = 1.0,
    )

Create the "AtmosGCMSpectra" `DiagnosticsGroup` which contains the
following diagnostic variables:

- spectrum_1d: 1D power spectrum of kinetic energy ( 1/2*u^2 + 1/2*v^2)
- spectrum_2d: 2D power spectrum of kinetic energy ( 1/2*u^2 + 1/2*v^2)

These are output with `m`, `lat`, `level` and `m_t`, `n`, `level`
dimensions, i.e. zonal wavenumber, latitude, level and (truncated) zonal
wavenumber, total wavenumber, level, as well as a (unlimited) `time`
dimension at the specified `interval`. The spectrum is computed after the
velocity fields have been interpolated (`interpol` _must_ be specified).
"""
function setup_atmos_spectra_diagnostics(
    ::AtmosGCMConfigType,
    interval::String,
    out_prefix::String;
    writer = NetCDFWriter(),
    interpol = nothing,
    nor = 1.0,
)
    @assert !isnothing(interpol)

    return DiagnosticsGroup(
        "AtmosGCMSpectra",
        Diagnostics.atmos_gcm_spectra_init,
        Diagnostics.atmos_gcm_spectra_fini,
        Diagnostics.atmos_gcm_spectra_collect,
        interval,
        out_prefix,
        writer,
        interpol,
        AtmosGCMSpectraDiagnosticsParams(nor),
    )
end

function get_spectra(mpicomm, mpirank, Q, bl, interpol, nor)
    FT = eltype(Q)
    if array_device(Q) isa CPU
        ArrayType = Array
    else
        ArrayType = CUDA.CuArray
    end

    istate = ArrayType{FT}(undef, interpol.Npl, number_states(bl, Prognostic()))
    interpolate_local!(interpol, Q.realdata, istate)

    # TODO: get indices here without hard-coding them
    _ρu, _ρv, _ρw = 2, 3, 4
    project_cubed_sphere!(interpol, istate, (_ρu, _ρv, _ρw))

    all_state_data = accumulate_interpolated_data(mpicomm, interpol, istate)

    if mpirank == 0

        dims = dimensions(interpol)
        lat = dims["lat"][1]
        lon = dims["long"][1]
        level = dims["level"][1] .- FT(planet_radius(Settings.param_set))

        mass_weight = ones(FT, length(level)) # TODO: interpolate on pressure levs and do mass weighted calculations

        var_grid_u = all_state_data[:, :, :, 2] ./ all_state_data[:, :, :, 1]
        spectrum_1d_u, m = power_spectrum_1d(
            AtmosGCMConfigType(),
            var_grid_u,
            level,
            lat,
            lon,
            mass_weight,
        )
        spectrum_2d_u, m_and_n, _, __ =
            power_spectrum_2d(AtmosGCMConfigType(), var_grid_u, mass_weight)

        var_grid_v = all_state_data[:, :, :, 3] ./ all_state_data[:, :, :, 1]
        mass_weight = ones(FT, length(level)) # TODO: interpolate on pressure levs and do mass weighted calculations
        spectrum_1d_v, m = power_spectrum_1d(
            AtmosGCMConfigType(),
            var_grid_v,
            level,
            lat,
            lon,
            mass_weight,
        )
        spectrum_2d_v, m_and_n, _, __ =
            power_spectrum_2d(AtmosGCMConfigType(), var_grid_v, mass_weight)

        spectrum_1d = 0.5 .* (spectrum_1d_u + spectrum_1d_v)
        spectrum_2d = 0.5 .* (spectrum_2d_u + spectrum_2d_v)

        return spectrum_1d, m, spectrum_2d, m_and_n
    end
    return nothing, nothing, nothing, nothing
end

function atmos_gcm_spectra_init(dgngrp, currtime)
    Q = Settings.Q
    bl = Settings.dg.balance_law
    mpicomm = Settings.mpicomm
    mpirank = MPI.Comm_rank(mpicomm)
    FT = eltype(Q)
    interpol = dgngrp.interpol
    nor = dgngrp.params.nor

    # get the 1d and 2d spectra and their associated wavenumbers
    spectrum_1d, m, spectrum_2d, m_and_n =
        get_spectra(mpicomm, mpirank, Q, bl, interpol, nor)

    if mpirank == 0
        # Setup dimensions
        interpol_dims = dimensions(interpol)
        lat = interpol_dims["lat"]
        level = interpol_dims["level"]
        level = (level[1] .- FT(planet_radius(Settings.param_set)), level[2])
        num_fourier = ((2 * length(lat[1]) - 1) / 3) # number of truncated zonal wavenumbers
        num_spherical = (num_fourier + 1) # number of total wavenumbers (n)
        dims = OrderedDict(
            "m" => (collect(FT, 1:1:length(lat[1])) .- 1, Dict()),
            "m_t" => (collect(FT, 0:1:num_fourier), Dict()),
            "n" => (collect(FT, 0:1:num_spherical), Dict()),
            "level" => level,
            "lat" => lat,
        )
        # Setup variables (NB: 2d spectrum is on the truncated wavenumber grid)
        vars = OrderedDict(
            "spectrum_1d" => (("m", "lat", "level"), FT, Dict()),
            "spectrum_2d" => (("m_t", "n", "level"), FT, Dict()),
        )

        dprefix = @sprintf(
            "%s_%s-%s",
            dgngrp.out_prefix,
            dgngrp.name,
            Settings.starttime,
        )
        dfilename = joinpath(Settings.output_dir, dprefix)
        init_data(dgngrp.writer, dfilename, dims, vars)
    end

    return nothing
end

function atmos_gcm_spectra_collect(dgngrp, currtime)
    Q = Settings.Q
    bl = Settings.dg.balance_law
    mpicomm = Settings.mpicomm
    mpirank = MPI.Comm_rank(mpicomm)
    FT = eltype(Q)
    interpol = dgngrp.interpol
    nor = dgngrp.params.nor

    spectrum_1d, _, spectrum_2d, _ =
        get_spectra(mpicomm, mpirank, Q, bl, interpol, nor)

    if mpirank == 0
        varvals = OrderedDict(
            "spectrum_1d" => spectrum_1d,
            "spectrum_2d" => spectrum_2d,
        )
        append_data(dgngrp.writer, varvals, currtime)
    end

    MPI.Barrier(mpicomm)
    return nothing
end

function atmos_gcm_spectra_fini(dgngrp, currtime) end

max_forks_repo_name: jm-c/CLIMA


text:
"""
General map utilities for rectangular pixels.
"""

using WCS  # wraps the WCSLIB C library
using FFTW

"""
Rectangular map type, contains an AbstractArray and a WCS object, but behaves like the
AbstractArray it contains for array operations.

It only implements the subset of Base.Array operations which are common on maps.
You should work with the data directly using `enmap_instance.data` if you need
additional Array functions.
"""
struct Enmap{T,N,AA<:AbstractArray} <: AbstractArray{T,N}
    data::AA  # some kind of abstract array
    wcs::WCSTransform  # WCS object from WCS.jl
end

# getters
Base.parent(A::Enmap) = A.data
get_wcs(A::Enmap) = A.wcs

# constructor
Enmap(A::AbstractArray{T,N}, wcs::WCSTransform) where {T,N} =
    Enmap{T,N,typeof(A)}(A, wcs)

# array details
parenttype(::Type{Enmap{T,N,AA}}) where {T,N,AA} = AA
parenttype(A::Enmap) = parenttype(typeof(A))
Base.eachindex(::IndexCartesian, A::Enmap) = CartesianIndices(axes(A))
Base.IndexStyle(::Type{EM}) where {EM<:Enmap} = IndexStyle(parenttype(EM))
Base.size(A::Enmap) = size(parent(A))
Base.size(A::Enmap, d) = size(parent(A), d)

# create a similar map by making a deep copy of the WCS
function Base.similar(A::Enmap, ::Type{T}, dims::Dims) where T
    B = Enmap(similar(parent(A), T, dims), deepcopy(get_wcs(A)))
end

# send all the get and set to the data inside the enmap
@inline function Base.getindex(A::Enmap{T,N,AA}, I::Vararg{Int,N}) where {T,N,AA}
    parent(A)[I...]
end
@inline function Base.getindex(A::Enmap{T,N,AA}, i::Int) where {T,N,AA}
    parent(A)[i]
end
@inline function Base.setindex!(A::Enmap{T,N,AA}, val, I::Vararg{Int,N}) where {T,N,AA}
    parent(A)[I...] = val
end
@inline function Base.setindex!(A::Enmap{T,N,AA}, val, i::Int) where {T,N,AA}
    parent(A)[i] = val
end

# area functions
function extent_intermediate(shape::Tuple, wcs::WCSTransform)
    return abs.(deg2rad.(shape .* wcs.cdelt))
end

function area(shape::Tuple, wcs::WCSTransform)
    return prod(extent_intermediate(shape, wcs))
end

function pixsize(shape::Tuple, wcs::WCSTransform)
    return area(shape, wcs) / prod(shape)
end

function pixsize(m::Enmap)
    return pixsize(size(m), m.wcs)
end

function get_fft_norm(x::Enmap)
    return sqrt(pixsize(x) / prod(size(x.data)))
end

# normalizations made to agree with enmap
function enfft(x::Enmap{T,N,AA}) where {T,N,AA}
    return Enmap(FFTW.fft(x.data) .* get_fft_norm(x),
        deepcopy(x.wcs))
end

function enifft(x::Enmap{T,N,AA}) where {T,N,AA}
    return Enmap(FFTW.ifft(x.data) .* (1.0 / get_fft_norm(x)), deepcopy(x.wcs))
end

"""
Physically normalized enmap FFT
"""
function enfft!(x::Enmap{T,N,AA}) where {T,N,AA}
    FFTW.fft!(x.data)
    x .*= get_fft_norm(x)
end

"""
Physically normalized enmap inverse FFT
"""
function enifft!(x::Enmap{T,N,AA}) where {T,N,AA}
    FFTW.ifft!(x.data)
    x .*= (1.0/get_fft_norm(x))
end

max_stars_repo_name: xzackli/XGPaint

max_issues_repo_name: xzackli/XGPaint

max_forks_repo_name: xzackli/XGPaint


text:
module Laszip

using Printf

export
	create_empty_header, create_empty_point, create_empty_vlr,
	laszip_get_point_pointer, laszip_get_header_pointer,
	laszip_get_version,
	laszip_create,
	laszip_get_error,
	laszip_get_warning,
	laszip_clean,
	laszip_destroy,
	laszip_get_header_pointer,
	laszip_get_point_pointer,
	laszip_get_point_count,
	laszip_set_header,
	laszip_check_for_integer_overflow,
	laszip_auto_offset,
	laszip_set_point,
	laszip_set_coordinates,
	laszip_get_coordinates,
	laszip_set_geokeys,
	laszip_set_geodouble_params,
	laszip_set_geoascii_params,
	laszip_add_vlr,
	laszip_open_writer,
	laszip_write_point,
	laszip_close_writer,
	laszip_open_reader,
	laszip_seek_point,
	laszip_read_point,
	laszip_close_reader,
# Non C-wrapper functions
	xyz2laz, laz2xyz, msgerror, las2dat, dat2las, create_header, create_point

include("laszip_h.jl")
include("laszip_dll.jl")
include("las2dat.jl")
include("dat2las.jl")

"""
	Prints an laszip error message and error out

	msgerror(lzobj::Ptr{Void}, extramsg::AbstractString="")

	Where:
		"lzobj" is a pointer to laszip_readert|writer created by unsafe_load(laszip_create(arg))
		"extramsg" is an optional extra message to be printed before the laszip error message.
"""
function msgerror(lzobj::Ptr{Cvoid}, extramsg::AbstractString="")
	pStr = pointer([pointer(lpad("",128,"        "))])		# Create a 1024 bytes string and get its pointer
	laszip_get_error(lzobj, pStr)
	Str = unsafe_string(unsafe_load(pStr))
	if (isempty(extramsg))
		error(Str)
	else
		error(extramsg * "\n\t" * Str)
	end
end
function msgerror(lzobj::Ptr{Ptr{Cvoid}}, extramsg::AbstractString="")
	lzobj = unsafe_load(lzobj)
	msgerror(lzobj, extramsg)
end

# ---------------------------------------------------------------------------------------------
function create_header(project_ID_GUID_data_4="", system_identifier="", generating_software="")
# Initialize an empty (or almost) header type
	if (isempty(project_ID_GUID_data_4))
		project_ID_GUID_data_4 = Laszip.Array_8_UInt8(fill(Base.zero(UInt8),8)...)
	end
	if (isempty(system_identifier))
		system_identifier = Laszip.Array_32_UInt8(fill(Base.zero(UInt8),32)...)
	end
	if (isempty(generating_software))
		generating_software = Laszip.Array_32_UInt8(fill(Base.zero(UInt8),32)...)
	end
	number_of_points_by_return = Laszip.Array_5_UInt32(fill(Base.zero(UInt32),5)...)
	extended_number_of_points_by_return = Laszip.Array_15_Culonglong(fill(Base.zero(Culonglong),15)...)
	user_data_in_header = pointer("")
	vlrs = pointer([laszip_vlr])
	user_data_after_header = pointer("")
	Laszip.laszip_header(0, 0, 0, 0, 0, project_ID_GUID_data_4, '1', '2', system_identifier, generating_software,
	                     0, 0, 0, 0, 0, '0', 0, 0, number_of_points_by_return, 1., 1., 1., 0., 0., 0.,
	                     0., 0., 0., 0., 0., 0., 0, 0, 0, 0, extended_number_of_points_by_return,
	                     0, user_data_in_header, vlrs, 0, user_data_after_header)
end

# ---------------------------------------------------------------------------------------------
function create_point()
# Initialize an empty point type
	rgb = Laszip.Array_4_UInt16(fill(Base.zero(UInt16),4)...)
	dummy = Laszip.Array_7_Cuchar(fill(Base.zero(Cuchar),7)...)
	wave_packet = Laszip.Array_29_Cuchar(fill(Base.zero(Cuchar),29)...)
	laszip_point(0, 0, 0, 0, '0', '0', '0', '0', '0', '0', '0', '0', 0, '0', 0, 0, '0', '0', '0', '0', '0', '0',
	             dummy, 0., rgb, wave_packet, 0, pointer(""))
end

end # module

max_stars_repo_name: joa-quim/Laszip

max_issues_repo_name: joa-quim/Laszip

max_forks_repo_name: joa-quim/Laszip


text:
## PART 4

println("\nPart 4");

bvals = [b-3:b+3];
println("\n\tEstimating σ for f1..f3 with bin around b=$b");

σf1jack = Float64[];
σf2jack = Float64[];
σf3jack = Float64[];

for bb in bvals
    jestb = jacknife(v, N2, bb);
    push!(σf1jack, jacknifeσ(f1, jestb));
    push!(σf2jack, jacknifeσ(f2, jestb));
    push!(σf3jack, jacknifeσ(f3, jestb));
end

σfjack = cat(2, σf1jack, σf2jack, σf3jack);

println("\n\tGenerating plots in σfvsb to show σf dependence on b.");
using Gadfly
σfvsb = Plot[];
for i in [1:3]
    push!(σfvsb, Gadfly.plot(
        x = bvals,
        y = σfjack[:,i],
        Geom.line,
        Geom.point,
        Guide.xlabel("Bin size b"),
        Guide.ylabel("Estimated Variance"),
        Guide.title("Estimated Variance of function f$i vs. Bin Size")
    ));
end

println("\n\tCompare σfjack to naive σ̂ from part 1:\n")
println("\t\tσ̂f1N2/σfjack[4,1] = $(σ̂f1N2/σfjack[4,1])");
println("\t\tσ̂f2N2/σfjack[4,2] = $(σ̂f2N2/σfjack[4,2])");
println("\t\tσ̂f3N2/σfjack[4,3] = $(σ̂f3N2/σfjack[4,3])");


text:
@with_kw struct PredictionOverlay <: SceneOverlay
    prediction::Array{Float64} = Array{Float64}()
    color::Colorant = RGBA(1., 0.0, 0.0, 0.2)
end

function AutoViz.render!(rendermodel::RenderModel, overlay::PredictionOverlay, scene::Scene, roadway::R) where R

    predictions = overlay.prediction
    ego = scene[findfirst(1, scene)]
    for i=1:size(predictions)[1]
      pos = predictions[i,:]
      ped = Vehicle(VehicleState(VecSE2(pos[1], pos[2], 1.57), 0.), AutomotivePOMDPs.PEDESTRIAN_DEF, 1)
      add_instruction!(rendermodel, render_vehicle, (ego.state.posG.x+ped.state.posG.x+ego.def.length/2-AutomotivePOMDPs.PEDESTRIAN_DEF.width/2, ego.state.posG.y+ped.state.posG.y, ped.state.posG.θ, ped.def.length, ped.def.width, overlay.color))
    end
    return rendermodel
end

function animate_record(scenes::Vector{Scene},dt::Float64, env::CrosswalkEnv, ego_vehicle::Vector{Vehicle}, sensor::GaussianSensor, sensor_o::Vector{Vector{Vehicle}},
                             risk::Vector{Float64}, ttc::Vector{Float64}, collision_rate::Vector{Float64}, brake_request::Vector{Bool}, prediction::Vector{Array{Float64}}, cam=FitToContentCamera(0.0))

    duration = length(scenes)*dt::Float64
    fps = Int(1/dt)

    function render_rec(t, dt)
        frame_index = Int(floor(t/dt)) + 1
        text_to_visualize = [string("v-ego: ", round(ego_vehicle[frame_index].state.v*3.6,digits=1), "km/h", " TTC: ", round(ttc[frame_index],digits=1), "s",
                                " Collision: ",  collision_rate[frame_index], " Risk: ",  round(risk[frame_index],digits=1),
                                " Brake request: ",  brake_request[frame_index])]

        sensor_overlay = GaussianSensorOverlay(sensor=sensor, o=sensor_o[frame_index])
        occlusion_overlay = OcclusionOverlay(obstacles=env.obstacles)
        prediction_overlay = PredictionOverlay(prediction=prediction[frame_index])
        text_overlay = TextOverlay(text=text_to_visualize,pos=VecE2(75.,8.),incameraframe=true,color=colorant"white",font_size=20)
        return render(scenes[frame_index], env, [prediction_overlay, sensor_overlay, occlusion_overlay, text_overlay, IDOverlay()], cam=cam)

    end
    return duration, fps, render_rec
end

@with_kw mutable struct ObservationCallback

    sensor_observations::Vector{Vector{Vehicle}}
    ego_vehicle::Vector{Vehicle}
    ego_a::Vector{Float64}
    collision::Vector{Bool}

    collision_rate::Vector{Float64}
    ttc::Vector{Float64}
    risk::Vector{Float64}
    brake_request::Vector{Bool}
    prediction_obstacle::Vector{Array{Float64}}

end

function AutomotiveDrivingModels.run_callback(
        callback::ObservationCallback,
        scenes::Vector{Scene},
        actions::Union{Nothing, Vector{Frame{A}}},
        roadway::R,
        models::Dict{I,M},
        tick::Int) where {A<:EntityAction,I,R,M<:DriverModel}

    push!(callback.sensor_observations, models[1].sensor_observations)
    push!(callback.ego_vehicle, models[1].ego_vehicle)
    push!(callback.ego_a, models[1].a_current)
    collision = is_crash(scenes[tick])
    push!(callback.collision, collision)

    push!(callback.collision_rate, models[1].collision_rate)
    push!(callback.ttc, models[1].ttc)
    push!(callback.risk, models[1].risk)
    push!(callback.brake_request, models[1].brake_request)
    push!(callback.prediction_obstacle, models[1].prediction_obstacle)

    return collision
end


text:
include("utils.jl")

obs(env::PendulumEnv, ::Nothing) = obs(env, [0.0, 0.0])
obs(env::PendulumEnv, (θ, θ̄)) = Flux.data(θ)

struct PendulumDrawParams <: AbstractDrawParams
    screen_height::UInt32
    screen_width::UInt32
    world_width::Float32
    scale::Float32
    arm_length::Float32
    arm_width::Float32
    axle_radius::Float32
end

PendulumDrawParams() =
    PendulumDrawParams(
        500,       # screen_height
        500,       # screen_width
        44f-1,     # world_width
        500/44f-1, # scale (screen_width / world_width)
        1f0,       # arm_length
        2f-1,      # arm_width
        5f-2       # axle_radius
    )

function Ctx(env::PendulumEnv, mode::Symbol=:no_render)
    if mode == :human_pane
        draw_params = PendulumDrawParams()
        viewer = CairoRGBSurface(draw_params.screen_width, draw_params.screen_height)

        CairoCtx(draw_params, viewer)
    elseif mode == :human_window
        draw_params = PendulumDrawParams()
        viewer = CairoRGBSurface(draw_params.screen_width, draw_params.screen_height)

        canvas = @GtkCanvas()
        canvas.backcc = CairoContext(viewer)
        win = GtkWindow(canvas, "Pendulum",
                draw_params.screen_width, draw_params.screen_height; resizable=false)
        show(canvas)
        visible(win, false)
        signal_connect(win, "delete-event") do widget, event
            ccall((:gtk_widget_hide_on_delete, Gtk.libgtk), Bool, (Ptr{GObject},), win)
        end

        GtkCtx(draw_params, canvas, win)
    elseif mode == :rgb
        draw_params = PendulumDrawParams()
        viewer = CairoRGBSurface(draw_params.screen_width, draw_params.screen_height)

        RGBCtx(draw_params, viewer)
    elseif mode == :no_render
        return
    else
        error("Unrecognized mode in Ctx(): $(mode)")
    end
end

function drawcanvas!(env::PendulumEnv, viewer::CairoContext, params::PendulumDrawParams)
    # Background
    set_source_rgb(viewer, 1, 1, 1)
    rectangle(viewer, 0, 0, params.screen_width, params.screen_height)
    fill(viewer)

    # Move to center of screen
    translate_dist = Pair(params.screen_width/2, params.screen_height/2)
    translate(viewer, translate_dist.first, translate_dist.second)

    # Rotate
    rotate(viewer, env.state[2] * env.dt)

    # Arm Start Circle
    set_source_rgb(viewer, 8f-1, 3f-1, 3f-1)
    move_to(viewer, 0, 0)
    arc(viewer, 0, 0, params.scale * params.arm_width/2, π, 2*π)

    # Arm Side 1
    rel_line_to(viewer, 0, params.scale * params.arm_length)

    # Arm End Circle
    arc(viewer, 0, params.scale * params.arm_length, params.scale * params.arm_width/2, 0, π)

    # Arm Side 2
    rel_line_to(viewer, 0, -params.scale * params.arm_length)

    # Fill arm
    fill(viewer)

    # Axle
    set_source_rgb(viewer, 0, 0, 0)
    circle(viewer, 0, 0, params.scale * params.axle_radius)
    fill(viewer)

    # Undo translation and rotation
    rotate(viewer, -env.state[2] * env.dt)
    translate(viewer, -translate_dist.first, -translate_dist.second)
end


text:
"""
arraydata = formatter(filename::String, nx::Int, ny::Int; HeaderLines::Int=0)
### Arguments
- filename : file path and file name
- nx, ny   : number of x- and y- grids
### Keyword argument
- HeaderLines : number of header lines. Default is zero (no header).
### Return value
- arraydata: mesh data formatted from fixed length array
"""
function formatter(filename::String, nx::Int, ny::Int; HeaderLines::Int=0)
    txt = readlines(filename)
    arraydata = [parse(Float64, txt[i+HeaderLines][8(j-1)+1:8j]) for i=1:length(txt), j=1:10]
    arraydata = permutedims(arraydata, [2,1])
    arraydata = reshape(vec(arraydata), (nx, ny))
    arraydata = reverse(permutedims(arraydata, [2,1]), dims=1)
    return arraydata
end


text:
module Detection

using Parameters
using ImageFiltering
using DataFrames

export PeakMesh, extract_sources

"""
    Detection.SourceFinder

Abstract super type for source detection algorithms used with [`extract_sources`](@ref).
"""
abstract type SourceFinder end

"""
    extract_sources(::SourceFinder, data, [error]; sorted=true)

Uses `method` to find and extract point-like sources.

Returns a `DataFrame` with positions and information related to the `method`. For instance, using `PeakMesh` returns a table a column for the peak values.

`data` is assumed to be background-subtracted. If `error` is provided it will be propagated into the detection algorithm. If `sorted` is `true` the sources will be sorted by their amplitude.

# See Also
* [Source Detection Algorithms](@ref)
"""
extract_sources

"""
    PeakMesh(box_size=(3, 3), nsigma=3.0)

Detect sources by finding peaks above a threshold in grids across the image.

This creates a pixel-wise threshold for sources by calculating `error * nsigma` when used with [`extract_sources`](@ref)
The peaks are found by searching the image in boxes of size `box_size`. If the maximum value in that box is greater than the threshold set above, the point is extracted.
"""
@with_kw struct PeakMesh <: SourceFinder
    box_size::NTuple{2,<:Integer} = (3, 3)
    nsigma::Float64 = 3
    PeakMesh(box_size::NTuple{2,<:Integer}, nsigma) = new(box_size, nsigma)
    PeakMesh(box_size::Integer, nsigma) = new((box_size, box_size), nsigma)
end

function extract_sources(alg::PeakMesh, data::AbstractMatrix{T}, error = zeros(T, size(data)), sort = true) where T
    threshold = @. error * alg.nsigma
    data_max = mapwindow(maximum, data, alg.box_size, border = Fill(zero(T)))

    rows = NamedTuple{(:x, :y, :value),Tuple{Int,Int,T}}[]
    @inbounds for idx in CartesianIndices(data)
        if data[idx] > threshold[idx] && data[idx] == data_max[idx]
            push!(rows, (x = idx.I[2], y = idx.I[1], value = data[idx]))
        end
    end
    sort && sort!(rows, by = row->row.value, rev = true)
    return DataFrame(rows)
end

end # module Detection


text:
using Statistics

function meantime(times::Array, dlm::String=":")
    c = π / (12 * 60 * 60)
    a = map(x -> parse.(Int, x), split.(times, dlm))
    ϕ = collect(3600t[1] + 60t[2] + t[3] for t in a)
    d = angle(mean(exp.(c * im * ϕ))) / 2π # days
    if d < 0 d += 1 end
    # Convert to h:m:s
    h = trunc(Int, d * 24)
    m = trunc(Int, d * 24 * 60) - h * 60
    s = trunc(Int, d * 24 * 60 * 60) - h * 60 * 60 - m * 60
    return "$h:$m:$s"
end

times = String["23:00:17", "23:40:20", "00:12:45", "00:17:19"]
mtime = meantime(times)
println("Times:")
println.(times)
println("Mean: $mtime")

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# exp_xform.jl

using MIRT: exp_xform
using Test: @test, @test_throws, @inferred

# test for given data type
function exp_xform_test( ; T::DataType = ComplexF32)

    modes = (:element, :row, :column)

    isinteractive() && (@info "1D tests with $T")
    N = 500
    M = 6000
    D = 3
    X = randn(T, N)
    U = randn(T, D, N)
    V = randn(T, D, M)
    y1 = @inferred exp_xform(X, U, V ; mode = :matrix)

    for mode in modes
        y2 = @inferred exp_xform(X, U, V ; mode=mode)
        @test y1 ≈ y2
    end

    isinteractive() && (@info "2D tests with $T")
    N = 10000
    M = 100
    D = 20
    L = 10
    X = randn(T, N, L)
    U = randn(T, D, N)
    V = randn(T, D, M)
    y1 = @inferred exp_xform(X, U, V ; mode = :matrix)

    for mode in modes
        y2 = @inferred exp_xform(X, U, V ; mode=mode)
        @test y1 ≈ y2
    end

    true
end

@test_throws String exp_xform(ones(2,2), ones(2,2), ones(2,2) ; mode=:bad)

# for T in (ComplexF32, ComplexF64)
	@test exp_xform_test() # T=T
# end


text:
module Thesis

using Reexport
@reexport using DataFrames, CSV, StatsPlots

end # module


text:
"""
    t_matrix(CapsuleParticle{T,2,Acoustic{T,2},Sphere{T,2}}, Acoustic{T,2}, ω, order)

The T-matrix for a 2D circlular capsule particle in an acoustic medium.
"""
function t_matrix(cap::CapsuleParticle{T,2,Acoustic{T,2},Sphere{T,2}}, medium::Acoustic{T,2}, ω::T, M::Integer)::Diagonal{Complex{T}} where T <: AbstractFloat

    k = ω / medium.c
    k0 = ω / cap.inner.medium.c
    k1 = ω / cap.outer.medium.c
    a0 = outer_radius(cap.inner)
    a1 = outer_radius(cap.outer)
    q  = (medium.ρ*medium.c)/(cap.outer.medium.ρ*cap.outer.medium.c)
    q0 = (cap.inner.medium.ρ*cap.inner.medium.c)/(cap.outer.medium.ρ*cap.outer.medium.c)

    Yn(n::Integer)   = hankelh1(n,k1*a1)*besselj(n,k1*a0) - hankelh1(n,k1*a0)*besselj(n,k1*a1)
    Yddn(n::Integer) = diffhankelh1(n,k1*a1)*diffbesselj(n,k1*a0) - diffhankelh1(n,k1*a0)*diffbesselj(n,k1*a1)

    Ydn(n::Integer,x::Complex{T},y::Complex{T}) = hankelh1(n,x)*diffbesselj(n,y) - diffhankelh1(n,y)*besselj(n,x)

    function Tn(n::Integer)::Complex{T}
        numer = Ydn(n,k*a1,k*a1) * (q0*besselj(n,k0*a0)*Ydn(n,k1*a1,k1*a0) - Yn(n)*diffbesselj(n,k0*a0))
        denom = diffbesselj(n,k0*a0)*
            (q*hankelh1(n,k*a1)*Ydn(n,k1*a0,k1*a1) + diffhankelh1(n,k*a1)*Yn(n)) +
            q0*besselj(n,k0*a0)*
            (q*hankelh1(n,k*a1)*Yddn(n) - diffhankelh1(n,k*a1)*Ydn(n,k1*a1,k1*a0))

         return (numer - denom*besselj(n,k*a1)) / (denom*hankelh1(n,k*a1))
     end

    # Get Tns for positive m
    Tns = map(Tn,0:M)

    return Diagonal(vcat(reverse(Tns), Tns[2:end]))
end

function internal_field(x::AbstractArray{T}, p::CapsuleParticle{T,2,Acoustic{T,2},Sphere{T,2}}, source::RegularSource{T,Acoustic{T,2}}, ω::T, scattering_coefficients::AbstractVector) where T

    Nh = Int((length(scattering_coefficients) - one(T))/T(2.0)) #shorthand

    if iszero(p.outer.medium.c) || isinf(abs((p.outer.medium.c)))
        return zero(Complex{T})
    elseif norm(x - origin(p)) > outer_radius(p)
        @warn "point $x not insider particle $p. Returning zero."
        return zero(Complex{T})
    end

    k = ω / source.medium.c
    k0 = ω / p.inner.medium.c
    k1 = ω / p.outer.medium.c
    a0 = outer_radius(p.inner)
    a1 = outer_radius(p.outer)
    q0 = (p.inner.medium.ρ*p.inner.medium.c) / (p.outer.medium.ρ*p.outer.medium.c)
    q = (source.medium.ρ*source.medium.c) / (p.outer.medium.ρ*p.outer.medium.c)

    Yn(n::Integer) = hankelh1(n,k1*a1)*besselj(n,k1*a0) - hankelh1(n,k1*a0)*besselj(n,k1*a1)
    Yddn(n::Integer) = diffhankelh1(n,k1*a1)*diffbesselj(n,k1*a0) - diffhankelh1(n,k1*a0)*diffbesselj(n,k1*a1)
    Ydn(n::Integer,x::Complex{T},y::Complex{T}) = hankelh1(n,x)*diffbesselj(n,y) - diffhankelh1(n,y)*besselj(n,x)

    denom(n::Integer) = q0 * besselj(n,a0*k0) *
            (q*besselj(n,a1*k)*Yddn(n) - diffbesselj(n,a1*k)*Ydn(n,a1*k1,a0*k1)) + diffbesselj(n,a0*k0) * (q*besselj(n,k*a1)*Ydn(n,a0*k1,a1*k1) + diffbesselj(n,k*a1)*Yn(n))
    forces = scattering_coefficients .* [Ydn(n,k*a1,k*a1)/denom(n) for n = -Nh:Nh];

    if norm(x - origin(p)) <= outer_radius(p.inner)
        coefs = - q0 * forces .* Ydn.(-Nh:Nh, a0*k1, a0*k1)
        basis = regular_basis_function(p.inner, ω)

        return sum(basis(Nh,x-origin(p)) .* coefs)
    else # calculate field in the mid region
        J_coefs = forces .* [
            q0*besselj(n, a0*k0)*diffhankelh1(n, a0*k1) - hankelh1(n, a0*k1)*diffbesselj(n, a0*k0)
        for n = -Nh:Nh]

        H_coefs = forces .* [
            besselj(n, a0*k1)*diffbesselj(n,a0*k0) - q0*besselj(n, a0*k0)*diffbesselj(n,a0*k1)
        for n = -Nh:Nh]

        J_basis = regular_basis_function(p.outer, ω)
        H_basis = outgoing_basis_function(p.outer.medium, ω)
        return sum(
            J_basis(Nh,x-origin(p)) .* J_coefs .+ H_basis(Nh,x-origin(p)) .* H_coefs
        )
    end
end


text:
module EquationsSolver

using Symbolics
import Symbolics:@variables,Equation

include("base.jl")

abstract type EquationsProblem end
struct LinearProblem <: EquationsProblem
    eqs::Vector{Equation}
    vars::Vector{Num}
    maxiters::Int
end

struct NonlinearProblem <: EquationsProblem
    eqs::Vector{Num}
    vars::Dict
    maxiters::Int
    abstol::Float64
end

function LinearProblem(eqs::Any, vars::Dict, maxiters=10000)
    checked_eqs = check_eqs(eqs)
    res = check_vars(checked_eqs, vars)
    return LinearProblem(eqs, collect(keys(vars)), maxiters)
end
function NLProblem(eqs::Any, vars::Dict; maxiters=10000,abstol = 1.0E-6)
    eqs = check_eqs(eqs)
    res = check_vars(eqs, vars)
    return NonlinearProblem(eqs, vars, maxiters,abstol)
end

export LinearProblem
export NLProblem
export solve
export @variables,Equation
include("solver.jl")
end


text:
"""
    diff(a::AbstractVector, b::AbstractVector)

Compute the difference between vector `a` and vector `b`, and return a `VectorDifference`.
"""
Base.diff(a::AbstractVector, b::AbstractVector) =
    diff(a, b, collect(1:size(a, 1)), collect(1:size(b, 1)))

"""
    diff(a::AbstractVector, b::AbstractVector, ia::AbstractVector, ib::AbstractVector)

Like [`diff`](@ref), but provide integer vectors that number the rows of the vectors `a`
and `b`. The vectors `ia` and `ib` represent the row numbers of `a` and `b` respectively.
The position of each vector element refers to the row index of `a` or `b`.
"""
function Base.diff(
    a::AbstractVector, b::AbstractVector, ia::AbstractVector, ib::AbstractVector
)
    if length(a) == 0 || length(b) == 0
        T = promote_type(eltype(a), eltype(b))
        vab = sparse(view(T[], :))
        va = view(a, :)
        vb = view(b, :)
        ei = Int[]
        length(a) == 0 && return VectorDifference(ei, ei, ei, vab, vb, va)
        length(b) == 0 && return VectorDifference(ei, ei, ei, vab, va, vb)
    end

    mapa = Dict(zip(ia, 1:length(ia)))
    mapb = Dict(zip(ib, 1:length(ib)))

    # Compute modified indices and values
    modinds = intersect(ia, ib)
    mapped_modindsa = getindex.(Ref(mapa), modinds)
    mapped_modindsb = getindex.(Ref(mapb), modinds)
    modvals = sparse(view(a, mapped_modindsa) - view(b, mapped_modindsb))

    # Compute added indices and values
    addinds = setdiff(ib, ia)
    mapped_addinds = getindex.(Ref(mapb), addinds)
    addvals = view(b, mapped_addinds)

    # Compute removed indices and values
    reminds = setdiff(ia, ib)
    mapped_reminds = getindex.(Ref(mapa), reminds)
    remvals = view(a, mapped_reminds)

    return VectorDifference(modinds, addinds, reminds, modvals, addvals, remvals)
end

"""
    diff(A::AbstractMatrix, B::AbstractMatrix)

Compute the difference between matrix `A` and matrix `B`, and return a `MatrixDifference`.
"""
function Base.diff(A::AbstractMatrix, B::AbstractMatrix)
    ia = collect(1:size(A, 1))
    ja = collect(1:size(A, 2))
    ib = collect(1:size(B, 1))
    jb = collect(1:size(B, 2))
    return diff(A, B, ia, ja, ib, jb)
end

"""
    diff(A::AbstractMatrix, B::AbstractMatrix, ia::AbstractVector, ja::AbstractVector, ib::AbstractVector, jb::AbstractVector)

Like [`diff`](@ref), but provide integer vectors that number the rows and columns of the
matrices `A` and `B`. The vector `ia` represents the row numbers of `A`, and the vector `jb`
represents the column numbers of `B` etc. The position of each vector element refers to the
row index (or column index respectively) of `A` or `B`.
"""
function Base.diff(
    A::AbstractMatrix,
    B::AbstractMatrix,
    ia::AbstractVector,
    ja::AbstractVector,
    ib::AbstractVector,
    jb::AbstractVector,
)
    if size(A) == (0, 0) || size(B) == (0, 0)
        T = promote_type(eltype(A), eltype(B))
        vab = sparse(view(T[], :))
        va = view(A, :)
        vb = view(B, :)
        ei = Int[]
        size(A) == (0, 0) &&
            return MatrixDifference((ei, ei), (ei, ei), (ei, ei), vab, vb, va)
        size(B) == (0, 0) &&
            return MatrixDifference((ei, ei), (ei, ei), (ei, ei), vab, va, vb)
    end

    mapia = Dict(zip(ia, 1:length(ia)))
    mapja = Dict(zip(ja, 1:length(ja)))
    mapib = Dict(zip(ib, 1:length(ib)))
    mapjb = Dict(zip(jb, 1:length(jb)))

    # Compute modified indices and values
    modindsi = intersect(ia, ib)
    modindsj = intersect(ja, jb)
    mapped_modindsia = getindex.(Ref(mapia), modindsi)
    mapped_modindsja = getindex.(Ref(mapja), modindsj)
    mapped_modindsib = getindex.(Ref(mapib), modindsi)
    mapped_modindsjb = getindex.(Ref(mapjb), modindsj)
    modvals = sparse(
        vec(
            view(A, mapped_modindsia, mapped_modindsja) -
            view(B, mapped_modindsib, mapped_modindsjb),
        ),
    )

    # Compute added indices and values
    addinds = (setdiff(ib, ia), setdiff(jb, ja))
    mapped_indicesb = CartesianIndices(B)
    mapped_modindicesb =
        CartesianIndex.(Iterators.product(mapped_modindsib, mapped_modindsjb))
    mapped_addindices = setdiff(mapped_indicesb, mapped_modindicesb)
    addvals = view(B, mapped_addindices)

    # Compute removed indices and values
    reminds = (setdiff(ia, ib), setdiff(ja, jb))
    mapped_indicesa = CartesianIndices(A)
    mapped_modindicesa =
        CartesianIndex.(Iterators.product(mapped_modindsia, mapped_modindsja))
    mapped_remindices = setdiff(mapped_indicesa, mapped_modindicesa)
    remvals = view(A, mapped_remindices)

    return MatrixDifference(
        (modindsi, modindsj), addinds, reminds, modvals, addvals, remvals
    )
end

"""
    diff(a::AbstractDict, b::AbstractDict)

Compute the difference between dictionary `a` and dictionary `b`, and return a
`DictDifference`.
"""
function Base.diff(a::AbstractDict, b::AbstractDict)
    keysa = keys(a)
    keysb = keys(b)
    modkeys = intersect(keysa, keysb)
    addkeys = setdiff(keysb, keysa)
    remkeys = setdiff(keysa, keysb)

    # Compute modified values
    modvals = Dict()
    for k in modkeys
        typeof(a[k]) != typeof(b[k]) &&
            throw(ArgumentError("type of values of common keys in `a` and `b` must match"))
        modvals[k] = typeof(a[k]) <: Number ? a[k] - b[k] : diff(a[k], b[k])
    end

    # Compute added values
    addvals = Dict()
    for k in addkeys
        addvals[k] = b[k]
    end

    # Compute removed values
    remvals = Dict()
    for k in remkeys
        remvals[k] = a[k]
    end

    return DictDifference(modvals, addvals, remvals)
end

"""
    diff(a::NamedTuple, b::NamedTuple)

Compute the difference between named tuple `a` and named tuple `b`, and return a
`NamedTupleDifference`.
"""
function Base.diff(a::NamedTuple, b::NamedTuple)
    modnames = intersect(keys(a), keys(b))
    addnames = setdiff(keys(b), keys(a))
    remnames = setdiff(keys(a), keys(b))

    # Compute modified values
    modvalues = []
    for n in modnames
        typeof(a[n]) != typeof(b[n]) &&
            throw(ArgumentError("type of values of common names in `a` and `b` must match"))
        v = typeof(a[n]) <: Number ? a[n] - b[n] : diff(a[n], b[n])
        push!(modvalues, v)
    end
    modvals = (; zip(modnames, modvalues)...)

    # Compute added values
    addvalues = []
    for n in addnames
        push!(addvalues, b[n])
    end
    addvals = (; zip(addnames, addvalues)...)

    # Compute removed values
    remvalues = []
    for n in remnames
        push!(remvalues, a[n])
    end
    remvals = (; zip(remnames, remvalues)...)

    return NamedTupleDifference(modvals, addvals, remvals)
end

"""
    diff(a::AbstractSet, b::AbstractSet)

Compute the difference between set `a` and set `b`, and return a `SetDifference`.

# Examples
```jldoctest
julia> diff(Set([1, 2, 3, 3]), Set([4, 2, 1]))
(Set([1, 2]), Set([4]), Set([3]))
```
"""
Base.diff(a::AbstractSet, b::AbstractSet) =
    SetDifference(intersect(a, b), setdiff(b, a), setdiff(a, b))


text:
function allbindings(pat, bs)
  if isa(pat, QuoteNode)
    return allbindings(pat.value, bs)
  end
  return isbinding(pat) || (isslurp(pat) && pat ≠ :__) ? push!(bs, bname(pat)) :
  isa(pat, TypeBind) ? push!(bs, pat.name) :
  isa(pat, OrBind) ? (allbindings(pat.pat1, bs); allbindings(pat.pat2, bs)) :
  istb(pat) ? push!(bs, tbname(pat)) :
  isexpr(pat, :$) ? bs :
  isa(pat, Expr) ? map(pat -> allbindings(pat, bs), [pat.head, pat.args...]) :
  bs
end

allbindings(pat) = (bs = Any[]; allbindings(pat, bs); bs)

function bindinglet(bs, body)
  ex = :(let $(esc(:env)) = env, $((:($(esc(b)) = get(env, $(Expr(:quote, b)), nothing)) for b in bs)...)
           $body
         end)
  return ex
end

function makeclause(pat, yes, els = nothing)
  bs = allbindings(pat)
  pat = subtb(subor(pat))
  quote
    env = trymatch($(Expr(:quote, pat)), ex)
    if env != nothing
      $(bindinglet(bs, esc(yes)))
    else
      $els
    end
  end
end

function clauses(ex)
  line = nothing
  clauses = []
  for l in ex.args
    isline(l) && (line = l; continue)
    env = trymatch(:(pat_ => yes_), l)
    env == nothing && error("Invalid match clause $l")
    pat, yes = env[:pat], env[:yes]
    push!(clauses, (pat, :($line;$yes)))
  end
  return clauses
end

macro match(ex, lines)
  @assert isexpr(lines, :block)
  result = quote
    ex = $(esc(ex))
  end

  @static if VERSION < v"0.7.0-"
    body = foldr((clause, body) -> makeclause(clause..., body),
                 nothing, clauses(lines))
  else
    body = foldr((clause, body) -> makeclause(clause..., body),
                 clauses(lines); init=nothing)
  end

  push!(result.args, body)
  return result
end

macro capture(ex, pat)
  bs = allbindings(pat)
  pat = subtb(subor(pat))
  quote
    $([:($(esc(b)) = nothing) for b in bs]...)
    env = trymatch($(esc(Expr(:quote, pat))), $(esc(ex)))
    if env == nothing
      false
    else
      $([:($(esc(b)) = get(env, $(esc(Expr(:quote, b))), nothing)) for b in bs]...)
      true
    end
  end
end


text:
# ---
# title: 968. Binary Tree Cameras
# id: problem968
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Dynamic Programming, Tree, Depth-first Search
# link: <https://leetcode.com/problems/binary-tree-cameras/description/>
# hidden: true
# ---
#
# Given a binary tree, we install cameras on the nodes of the tree.
#
# Each camera at a node can monitor **its parent, itself, and its immediate
# children**.
#
# Calculate the minimum number of cameras needed to monitor all nodes of the
# tree.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2018/12/29/bst_cameras_01.png)
#
#
#
#     Input: [0,0,null,0,0]
#     Output: 1
#     Explanation: One camera is enough to monitor all nodes if placed as shown.
#
#
# **Example 2:**
#
# ![](https://assets.leetcode.com/uploads/2018/12/29/bst_cameras_02.png)
#
#
#
#     Input: [0,0,null,0,null,0,null,null,0]
#     Output: 2
#     Explanation: At least two cameras are needed to monitor all nodes of the tree. The above image shows one of the valid configurations of camera placement.
#
#
#
# **Note:**
#
#   1. The number of nodes in the given tree will be in the range `[1, 1000]`.
#   2. **Every** node has value 0.
#
#
## @lc code=start
using LeetCode

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


text:
using Asap, JSON, Statistics

# Taken from NCCR compas.dev example files
test1 = JSON.parsefile("examples/compas_fea_truss_frame.json")

# nodes (in mm)
ns = test1["node_list"]

# elements
es = test1["element_list"]

# material properties
mp = test1["material_properties"]
E = Float64(mp["youngs_modulus"]) # kn/cm2
A = mp["cross_sec_area"] #cm^2
G = Float64(mp["shear_modulus"]) #kN/cm^2
Iz = mp["Iz"] #cm^4
Iy = mp["Iy"] # cm^4
J = mp["Jx"] #cm^4

# Create nodes
nodes = Vector{Node}()

for i = 1:length(ns)
    n = ns[i]
    position = [n["point"]["X"], n["point"]["Y"], n["point"]["Z"]] ./ 10

    if n["is_grounded"] > 0
        fixtype = :fixed
    else
        fixtype = :free
    end
    push!(nodes, Node(position, :frame, fixtype))
end

# Create elements
elements = Vector{Element}()
for e in es
    indices = e["end_node_ids"] .+ 1
    push!(elements, Element(nodes, indices, E, A, G, Iz, Iy, J))
end

# Create self-weight loads
lengths = [e.length for e in elements]
P = mean(lengths) * A * mp["density"]/100^3 / 2
Pvec = [0, 0, -P, 0, 0, 0]

loads = Vector{Load}()
for node in nodes
    push!(loads, Load(nodes, node.position, Pvec))
end

structure = Structure(nodes, elements, loads)

analyze(structure)

geo = Geometry(structure)
structurePlot(geo; scaleToForce = true)

max_stars_repo_name: keithjlee/Asap

max_issues_repo_name: keithjlee/Asap

max_forks_repo_name: keithjlee/Asap


text:
using Spec
using ..Tagging, ..Traits, ..IDS
export appendscope, rmscope, scope

"append `id` to the scope"
appendscope(ω::Ω, id) where {Ω <: AbstractΩ} =
  appendscope(traits(Ω), ω, id)
appendscope(::trait(Scope), ω, id) =
  updatetag(ω, Val{:scope}, append(id, ω.tags.scope))
appendscope(traits, ω, id) =
  tag(ω, (scope = id,))

rmscope(ω::T) where T = rmscope(ω, traithastag(T, Val{:scope}))
rmscope(ω, ::HasTag{:scope}) = rmtag(ω, Val{:scope})
rmscope(ω, _) = ω

""""
`scope(ω)`

Current scope
"""
scope(ω) = ω.tags.scope
@pre scope(ω) = hastag(ω, :scope) "Tag should have socpe"

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect