text:
@testset "Clebsch-Gordan" begin

using PyCall, Test, SHIPs.SphericalHarmonics
using SHIPs: eval_basis
using SHIPs.SphericalHarmonics: index_y
using SHIPs.Rotations: ClebschGordan, clebschgordan

sympy = pyimport("sympy")
spin = pyimport("sympy.physics.quantum.spin")

pycg(j1, m1, j2, m2, j3, m3, T=Float64) =
      spin.CG(j1, m1, j2, m2, j3, m3).doit().evalf().__float__()

cg = ClebschGordan()

ntest = 0
while ntest <= 200
   j1 = rand(0:10)
   j2 = rand(0:10)
   J = rand(abs(j1-j2):min(10,(j1+j2)))
   m1 = rand(-j1:j1)
   for m2 = -j2:j2
      M = m1+m2
      if abs(M) <= J
         ntest += 1
         print_tf(@test clebschgordan(j1,m1,j2,m2,J, M) ≈ pycg(j1,m1, j2,m2, J, M))
         print_tf(@test clebschgordan(j1,m1,j2,m2,J, M) ≈ cg(j1,m1, j2,m2, J, M))
      end
   end
end
println()

@info("Checking the SphH expansion in terms of CG coeffs")
# expansion coefficients of a product of two spherical harmonics in terms a
# single spherical harmonic
# see e.g. https://en.wikipedia.org/wiki/Clebsch–Gordan_coefficients
# this is the magic formula that we need
for ntest = 1:100
   # two random Ylm  ...
   l1, l2 = rand(1:10), rand(1:10)
   m1, m2 = rand(-l1:l1), rand(-l2:l2)
   # ... evaluated at random spherical coordinates
   θ = rand() * π
   φ = (rand()-0.5) * 2*π
   R = SVector( cos(φ)*sin(θ), sin(φ)*sin(θ), cos(θ) )
   # evaluate all relevant Ylms (up to l1 + l2)
   Ylm = eval_basis(SHBasis(l1+l2), R)
   # evaluate the product p = Y_l1_m1 * Y_l2_m2
   p = Ylm[index_y(l1,  m1)] * Ylm[index_y(l2,m2)]
   # and its expansion in terms of CG coeffs
   p2 = 0.0
   M = m1 + m2  # all other coeffs are zero

   for L = abs(M):(l1+l2)
      p2 += sqrt( (2*l1+1)*(2*l2+1) / (4 * π * (2*L+1)) ) *
            clebschgordan(l1,  0, l2,  0, L, 0) *
            clebschgordan(l1, m1, l2, m2, L, M) *
            Ylm[index_y(L, M)]
   end
   print_tf((@test (p ≈ p2) || (abs(p-p2) < 1e-15)))
end
println()

end # @testset

max_stars_repo_name: cortner/SHIPs

max_issues_repo_name: cortner/SHIPs

max_forks_repo_name: cortner/SHIPs


text:
# Author: Youngjun Kim, youngjun@stanford.edu
# Date: 06/05/2014

module CollisionAvoidanceSystem

export
    AbstractCollisionAvoidanceSystem,

    initialize,
    update,

    testThreat,
    selectRA,

    SimpleTCAS,
    SimpleTCASInput,
    SimpleTCASResolutionAdvisory,

    CoordSimpleTCAS,

    ACASX_CCAS,
    ACASX_ADD,
    ACASXInput,
    ACASXOutput

import CommonInterfaces: initialize, update

using AbstractCollisionAvoidanceSystemImpl

using SimpleTCASImpl

using CoordSimpleTCASImpl

using ACASXCommonImpl

pkgs = keys(Pkg.installed())
if "CCAS" in pkgs && "CASInterface" in pkgs
    @eval using ACASX_CCAS_Impl
end

using ACASX_ADD_Impl

end

ext: jl

max_stars_repo_name: sisl/RLES-SISLES

max_issues_repo_name: sisl/RLES-SISLES

max_forks_repo_name: sisl/RLES-SISLES


text:
let
    #######################################
    # First test the subproblem.
    srand(42)
    n = 5
    H = rand(n, n)
    H = H' * H + 4 * eye(n)
    H_eig = eigfact(H)
    U = H_eig[:vectors]

    gr = zeros(n)
    gr[1] = 1.
    s = zeros(Float64, n)

    true_s = -H \ gr
    s_norm2 = dot(true_s, true_s)
    true_m = dot(true_s, gr) + 0.5 * dot(true_s, H * true_s)

    # An interior solution
    delta = sqrt(s_norm2) + 1.0
    m, interior, lambda, hard_case, reached_solution =
        Optim.solve_tr_subproblem!(gr, H, delta, s)
    @assert interior
    @assert !hard_case
    @assert reached_solution
    @assert abs(m - true_m) < 1e-12
    @assert norm(s - true_s) < 1e-12
    @assert abs(lambda) < 1e-12

    # A boundary solution
    delta = 0.5 * sqrt(s_norm2)
    m, interior, lambda, hard_case, reached_solution =
        Optim.solve_tr_subproblem!(gr, H, delta, s)
    @assert !interior
    @assert !hard_case
    @assert reached_solution
    @assert m > true_m
    @assert abs(norm(s) - delta) < 1e-12
    @assert lambda > 0

    # A "hard case" where the gradient is orthogonal to the lowest eigenvector

    # Test the checking
    hard_case, lambda_1_multiplicity =
        Optim.check_hard_case_candidate([-1., 2., 3.], [0., 1., 1.])
    @assert hard_case
    @assert lambda_1_multiplicity == 1

    hard_case, lambda_1_multiplicity =
        Optim.check_hard_case_candidate([-1., -1., 3.], [0., 0., 1.])
    @assert hard_case
    @assert lambda_1_multiplicity == 2

    hard_case, lambda_1_multiplicity =
        Optim.check_hard_case_candidate([-1., -1., -1.], [0., 0., 0.])
    @assert hard_case
    @assert lambda_1_multiplicity == 3

    hard_case, lambda_1_multiplicity =
        Optim.check_hard_case_candidate([1., 2., 3.], [0., 1., 1.])
    @assert !hard_case

    hard_case, lambda_1_multiplicity =
        Optim.check_hard_case_candidate([-1., -1., -1.], [0., 0., 1.])
    @assert !hard_case

    hard_case, lambda_1_multiplicity =
        Optim.check_hard_case_candidate([-1., 2., 3.], [1., 1., 1.])
    @assert !hard_case

    # Now check an actual had case problem
    L = zeros(Float64, n) + 0.1
    L[1] = -1.
    H = U * diagm(L) * U'
    H = 0.5 * (H' + H)
    @assert issymmetric(H)
    gr = U[:,2][:]
    @assert abs(dot(gr, U[:,1][:])) < 1e-12
    true_s = -H \ gr
    s_norm2 = dot(true_s, true_s)
    true_m = dot(true_s, gr) + 0.5 * dot(true_s, H * true_s)

    delta = 0.5 * sqrt(s_norm2)
    m, interior, lambda, hard_case, reached_solution =
        Optim.solve_tr_subproblem!(gr, H, delta, s)
    @assert !interior
    @assert hard_case
    @assert reached_solution
    @assert abs(lambda + L[1]) < 1e-4
    @assert abs(norm(s) - delta) < 1e-12

    #######################################
    # Next, test on actual optimization problems.

    function f(x::Vector)
        (x[1] - 5.0)^4
    end

    function g!(x::Vector, storage::Vector)
        storage[1] = 4.0 * (x[1] - 5.0)^3
    end

    function h!(x::Vector, storage::Matrix)
        storage[1, 1] = 12.0 * (x[1] - 5.0)^2
    end

    d = TwiceDifferentiableFunction(f, g!, h!)

    results = Optim.optimize(d, [0.0], method=NewtonTrustRegion())
    @assert length(results.trace) == 0
    @assert results.g_converged
    @assert norm(results.minimum - [5.0]) < 0.01

    eta = 0.9

    function f_2(x::Vector)
        0.5 * (x[1]^2 + eta * x[2]^2)
    end

    function g!_2(x::Vector, storage::Vector)
        storage[1] = x[1]
        storage[2] = eta * x[2]
    end

    function h!_2(x::Vector, storage::Matrix)
        storage[1, 1] = 1.0
        storage[1, 2] = 0.0
        storage[2, 1] = 0.0
        storage[2, 2] = eta
    end

    d = TwiceDifferentiableFunction(f_2, g!_2, h!_2)

    results = Optim.optimize(d, Float64[127, 921], method=NewtonTrustRegion())
    @assert results.g_converged
    @assert norm(results.minimum - [0.0, 0.0]) < 0.01

    # Test Optim.newton for all twice differentiable functions in
    # Optim.UnconstrainedProblems.examples
    for (name, prob) in Optim.UnconstrainedProblems.examples
    	if prob.istwicedifferentiable
    		ddf = DifferentiableFunction(prob.f, prob.g!)
    		res = Optim.optimize(ddf, prob.initial_x, NewtonTrustRegion(), OptimizationOptions(autodiff = true))
    		@assert norm(Optim.minimizer(res) - prob.solutions) < 1e-2
    		res = Optim.optimize(ddf.f, prob.initial_x, NewtonTrustRegion(), OptimizationOptions(autodiff = true))
    		@assert norm(Optim.minimizer(res) - prob.solutions) < 1e-2
            res = Optim.optimize(ddf.f, ddf.g!, prob.initial_x, NewtonTrustRegion(), OptimizationOptions(autodiff = true))
    		@assert norm(Optim.minimizer(res) - prob.solutions) < 1e-2
    	end
    end
end


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

using Markdown
using InteractiveUtils

# ╔═╡ dcab9a18-5d58-4632-a325-d16af9c79344
using Pkg, DrWatson

# ╔═╡ c95f3f02-f179-4f81-9c53-6e6b87f25613
begin
	@quickactivate "NumericalMethodsforEngineers"
	using NumericalMethodsforEngineers, Symbolics, SymbolicUtils
	using StaticArrays
end

# ╔═╡ c9036632-adb9-11eb-33d1-c5193729e72a
function lagrangianpolynomial1(x::Vector{Float64},
		y::Vector{Float64}, xi::Vector{Float64})

	np = length(x)
 	yi = Vector{Float64}(undef, length(xi))
 	m = sort(hcat(x, y), dims=1)
	for l in 1:length(xi)
		for i in 1:np
			term = 1.0
			for j in 1:np
				j !== i && (term *= (xi[l] - m[j, 1])/(m[i, 1] - m[j, 1]))
			end
			yi[l] += term * m[i, 2]
		end
	end
	columns=[]
	push!(columns, m[:,1])
	push!(columns, m[:,2])
	df = DataFrame(columns, [:x, :y])

	columns=[]
	push!(columns, xi)
	push!(columns, yi)
	dfxi = DataFrame(columns, [:xi, :yi])

	(df, dfxi)
end

# ╔═╡ 298bbd55-d7b6-4229-8848-f0ea11f6448e

# ╔═╡ b8e7c2c7-5a72-4d62-9c2d-0b6cfa3e8fbb
@variables x y;

# ╔═╡ 057c029c-bd62-4ecb-9dc8-ed0f35684fd0
to_compute = [x^2 + y, y^2 + x];

# ╔═╡ d8e6cfe2-4cab-476d-8730-ec23470ac1b7
f_expr = build_function(to_compute, [x, y]);

# ╔═╡ dcb5841b-9ab1-4efe-928a-a69a02dffc95
begin
	myf = eval(f_expr[1])
	myf(SA[2.0, 3.0])
end

# ╔═╡ 2f360a8b-4ab6-4617-9481-1a10eab71579
Base.remove_linenums!(f_expr[2])

# ╔═╡ e5f5c803-d739-460f-bdf2-c9f5b55786a3
begin
	myf! = eval(f_expr[2])
	out = zeros(2)
	myf!(out, [2.0, 3.0])
	out
end

# ╔═╡ add02f10-37f9-49b1-80e4-205758be2605
f1 = build_function(to_compute, [x, y], expression=Val{false});

# ╔═╡ 2c817c8d-d108-432c-98b4-f849ec2c9c71
f1[1]([2.0, 3.0])

# ╔═╡ Cell order:
# ╠═dcab9a18-5d58-4632-a325-d16af9c79344
# ╠═c95f3f02-f179-4f81-9c53-6e6b87f25613
# ╠═c9036632-adb9-11eb-33d1-c5193729e72a
# ╠═298bbd55-d7b6-4229-8848-f0ea11f6448e
# ╠═b8e7c2c7-5a72-4d62-9c2d-0b6cfa3e8fbb
# ╠═057c029c-bd62-4ecb-9dc8-ed0f35684fd0
# ╠═d8e6cfe2-4cab-476d-8730-ec23470ac1b7
# ╠═dcb5841b-9ab1-4efe-928a-a69a02dffc95
# ╠═2f360a8b-4ab6-4617-9481-1a10eab71579
# ╠═e5f5c803-d739-460f-bdf2-c9f5b55786a3
# ╠═add02f10-37f9-49b1-80e4-205758be2605
# ╠═2c817c8d-d108-432c-98b4-f849ec2c9c71


text:
# deprecated functions

using Base: @deprecate

@deprecate add!(a::Matrix, b::AbstractPDMat) pdadd!(a, b)

@deprecate add_scal!(a::Matrix, b::AbstractPDMat, c::Real) pdadd!(a, b, c)

@deprecate add_scal(a::Matrix, b::AbstractPDMat, c::Real) pdadd(a, b, c)

@deprecate full(x::AbstractPDMat) Matrix(x)

@deprecate CholType Cholesky

@deprecate ScalMat(d::Int, x::Real, inv_x::Real) ScalMat(d, x)
@deprecate PDiagMat(v::AbstractVector, inv_v::AbstractVector) PDiagMat(v)


text:
using KernelAbstractions
using CUDA
using Parameters
using ..Architectures: devi, default_architecture

"""
Given the hitran data and necessary parameters, calculate an absorption cross-section at the
given pressure, temperature, and grid of wavenumbers (or wavelengths)
"""
function compute_absorption_cross_section(
                # Required
                model::HitranModel,          # Model to use in this cross section calculation
                                             # (Calculation from Hitran data vs. using Interpolator)
                # Wavelength [nm] or wavenumber [cm-1] grid (modify using wavelength_flag)
                grid::Union{AbstractRange{<:Real}, AbstractArray},  # Can be range OR array
                pressure::Real,              # actual pressure [hPa]
                temperature::Real;           # actual temperature [K]
                # Optionals
                wavelength_flag::Bool=false  # Use wavelength in nm (true) or wavenumber cm-1 units (false)
                )

    @unpack hitran, broadening, wing_cutoff, vmr, CEF, architecture = model

    # Notify user of wavelength grid
    if (wavelength_flag)
        @info """
        Note: Cross-section reported to wavelength grid
        (will internally convert to wavenumber grid for calculations)
        """  maxlog = 5
    end

    # Convert T to float type (ex. if Int)
    temperature = AbstractFloat(temperature)

    # Store results here to return
    result = array_type(architecture)(zeros(eltype(temperature), length(grid)))
    fill!(result, 0);

    # Calculate the minimum and maximum grid bounds, including the wing cutoff
    grid_max = maximum(grid) + wing_cutoff
    grid_min = minimum(grid) - wing_cutoff

    # Convert to wavenumber from [nm] space if necessary
    grid = wavelength_flag ? reverse(nm_per_m ./ grid) : grid
    grid_min, grid_max = wavelength_flag ? (nm_per_m /grid_max, nm_per_m/grid_min) : (grid_min, grid_max)

    # Interpolators from grid bounds to index values
    grid_idx_interp_low  = LinearInterpolation(grid, 1:1:length(grid), extrapolation_bc=1)
    grid_idx_interp_high = LinearInterpolation(grid, 1:1:length(grid), extrapolation_bc=length(grid))

    # Temporary storage array for output of qoft!. Compiler/speed issues when returning value in qoft
    rate = zeros(eltype(temperature), 1)

    # Declare the device being used
    device = devi(architecture)

    grid = array_type(architecture)(grid)

    # Loop through all transition lines:
    for j in eachindex(hitran.Sᵢ)

        # Test that this ν lies within the grid
        # (Aside: I ♥ chained comparisons)
        if grid_min < hitran.νᵢ[j] < grid_max

            # Apply pressure shift
            ν   = hitran.νᵢ[j] + pressure / p_ref * hitran.δ_air[j]

            # Compute Lorentzian HWHM
            γ_l = (hitran.γ_air[j] *
                  (1 - vmr) * pressure / p_ref + hitran.γ_self[j] *
                  vmr * pressure / p_ref) *
                  (t_ref / temperature)^hitran.n_air[j]

            # Compute Doppler HWHM
            γ_d = ((cSqrt2Ln2 / cc_) * sqrt(cBolts_ / cMassMol) * sqrt(temperature) *
            hitran.νᵢ[j] / sqrt(mol_weight(hitran.mol[j], hitran.iso[j])))

            # Ratio of widths
            y = sqrt(cLn2) * γ_l / γ_d

            # Apply line intensity temperature corrections
            S = hitran.Sᵢ[j]
            if hitran.E″[j] != -1
                qoft!(hitran.mol[j], hitran.iso[j], temperature, t_ref, rate)
                S = S * rate[1] *
                        exp(c₂ * hitran.E″[j] * (1 / t_ref - 1 / temperature)) *
                        (1 - exp(-c₂ * hitran.νᵢ[j] / temperature)) / (1 - exp(-c₂ * hitran.νᵢ[j] / t_ref));

            end

            # Calculate index range that this ν impacts
            ind_start = Int64(round(grid_idx_interp_low(ν - wing_cutoff)))
            ind_stop  = Int64(round(grid_idx_interp_high(ν + wing_cutoff)))

            # Create views from the result and grid arrays
            result_view   = view(result, ind_start:ind_stop);
            grid_view     = view(grid, ind_start:ind_stop);

            # Kernel for performing the lineshape calculation
            kernel! = line_shape!(device)

            # Run the event on the kernel
            # That this, this function adds to each element in result, the contribution from this transition
            event = kernel!(result_view, array_type(architecture)(grid_view), ν, γ_d, γ_l, y, S, broadening, CEF, ndrange=length(grid_view))
            wait(device, event)
            synchronize_if_gpu()
        end
    end

    # Return the resulting lineshape
    return (wavelength_flag ? reverse(result) : result)
end

"""
    $(FUNCTIONNAME)(model::InterpolationModel, grid::AbstractRange{<:Real}, pressure::Real, temperature::Real; wavelength_flag::Bool=false)

Given an Interpolation Model, return the interpolated absorption cross-section at the given pressure,
temperature, and grid of wavenumbers (or wavelengths)

"""
function compute_absorption_cross_section(
    # Required
    model::InterpolationModel,      # Model to use in this cross section calculation
                                    # (Calculation from Interpolator vs. Hitran Data)
    # Wavelength [nm] or wavenumber [cm-1] grid
    grid::Union{AbstractRange{<:Real}, AbstractArray},  # Can be range OR array
    pressure::Real,                 # actual pressure [hPa]
    temperature::Real;              # actual temperature [K]
    # Optionals
    wavelength_flag::Bool=false,    # Use wavelength in nm (true) or wavenumber cm-1 units (false)
    )

    # Convert to wavenumber from [nm] space if necessary
    grid = wavelength_flag ? reverse(nm_per_m ./ collect(grid)) : collect(grid)

    # Scale the interpolation to match the model grids
    sitp = scale(model.itp, model.p_grid, model.t_grid, model.ν_grid)

    # Perform the interpolation and return the resulting grid
    return sitp(pressure, temperature, grid)
end

#####
##### Lineshape functions that are called by absorption_cross_section
#####

@kernel function line_shape!(A, @Const(grid), ν, γ_d, γ_l, y, S, ::Doppler, CEF)
    FT = eltype(ν)
    I = @index(Global, Linear)
    @inbounds A[I] += FT(S) * FT(cSqrtLn2divSqrtPi) * exp(-FT(cLn2) * ((FT(grid[I]) - FT(ν)) / FT(γ_d))^2) / FT(γ_d)
end

@kernel function line_shape!(A, @Const(grid), ν, γ_d, γ_l, y, S, ::Lorentz, CEF)
    FT = eltype(ν)
    I = @index(Global, Linear)
    @inbounds A[I] += FT(S) * FT(γ_l) / (FT(pi) * (FT(γ_l)^2 + (FT(grid[I]) - FT(ν))^2))
end

@kernel function line_shape!(A, @Const(grid), ν, γ_d, γ_l, y, S, ::Voigt, CEF)
    FT = eltype(ν)
    I = @index(Global, Linear)
    @inbounds A[I] += FT(S) * FT(cSqrtLn2divSqrtPi) / FT(γ_d) * real(w(CEF, FT(cSqrtLn2) / FT(γ_d) * (FT(grid[I]) - FT(ν)) + im * FT(y)))
end

@kernel function line_shape32!(A, @Const(grid), ν, γ_d, γ_l, y, S, broadening, CEF)
    line_shape!(A, grid, Float32(ν), Float32(γ_d), Float32(γ_l), Float32(y), Float32(S), broadening, CEF)
end


text:
##### Beginning of file

module Run # Begin submodule Snapshots.Run

__precompile__(true)

import ArgParse
import Dates
import HTTP
import Pkg
import TimeZones

import ..Types
import ..Utils
import ..Common

import ..delayederror
import ..process_delayed_error_list

function run_snapshots!!(
        ;
        src_provider,
        dst_provider,
        git_user_name,
        git_user_email,
        include_branches,
        exclude_branches,
        task::String = "all",
        is_dry_run::Bool = false,
        do_not_push_to_these_destinations::Vector{String} =
            String[],
        do_not_try_url_list::Vector{String} =
            String[],
        try_but_allow_failures_url_list::Vector{String} =
            String[],
        time_zone::Dates.TimeZone =
            Dates.TimeZone("America/New_York"),
        )::Nothing
    @info("Running Snapshots.Run.run_snapshots!!")

    all_src_organization_repos::Vector{String} = src_provider(:list_all_repos)()

    unique!(all_src_organization_repos)
    sort!(all_src_organization_repos)
    @debug(
        string(
            "All repos in the source organization ",
            "($(length(all_src_organization_repos))):",
            )
        )
    for i = 1:length(all_src_organization_repos)
        @debug(
            "$(i). $(all_src_organization_repos[i])"
            )
    end

    if task == "all"
        task_src_repos = all_src_organization_repos
    elseif Types._is_interval(task)
        task_interval::Types.AbstractInterval = Types._construct_interval(
            task
            )
        @info(string("Using interval for stage 2: "),task_interval,)
        task_src_repos = Common._names_that_fall_in_interval(
            all_src_organization_repos,
            task_interval,
            )
    else
        @warn("not a valid task: ", task,)
        delayederror("\"$(task)\" is not a valid task")
    end

    unique!(task_src_repos)
    sort!(task_src_repos)
    @debug(
        string(
            "Source repos in the task interval ",
            "($(length(task_src_repos))):",
            )
        )
    for i = 1:length(task_src_repos)
        @debug(
            "$(i). $(task_src_repos[i])"
            )
    end

    n = length(task_src_repos)
    for i = 1:n
        repo_name = strip(task_src_repos[i])
        @debug(
            string(
                "Processing \"$(repo_name)\" ",
                "(repo $(i) of $(n))",
                )
            )
        if false
        else
            if repo_name in do_not_push_to_these_destinations
            else
                Common._snapshot_repo!!(
                    repo_name;
                    src_provider = src_provider,
                    dst_provider = dst_provider,
                    include_branches = include_branches,
                    exclude_branches = exclude_branches,
                    git_user_name = git_user_name,
                    git_user_email = git_user_email,
                    time_zone = time_zone,
                    is_dry_run = is_dry_run,
                    )
            end
        end
    end

    @info(
        string(
            "SUCCESS: run_snapshots completed ",
            "successfully :) Good-bye!",
            )
        )

    process_delayed_error_list()
    return nothing
end

end # End submodule Snapshots.Run

##### End of file


text:
using OhMyArtifacts, Pkg, Dates
using Scratch
using Test
include("utils.jl")

# Set to true for verbose Pkg output
const verbose = false
global const pkgio = verbose ? stderr : (VERSION < v"1.6.0-DEV.254" ? mktemp()[2] : devnull)

@testset "OhMyArtifacts.jl Basics" begin
    temp_pkg_dir() do project_dir
        artifacts_toml = @my_artifacts_toml!()
        @test isdir(OhMyArtifacts.get_scratch_dir())
        @test isdir(OhMyArtifacts.get_artifacts_dir())
        @test isfile(artifacts_toml)
        @test startswith(artifacts_toml, OhMyArtifacts.get_scratch_dir())
        @test isempty(OhMyArtifacts.load_my_artifacts_toml(artifacts_toml))

        hash_a = create_my_artifact() do artifact_dir
            cp(@__FILE__, joinpath(artifact_dir, "a"))
        end
        hash_b = create_my_artifact() do artifact_dir
            cp(joinpath(@__DIR__, "utils.jl"), joinpath(artifact_dir, "a"))
        end
        @test isfile(joinpath(OhMyArtifacts.get_artifacts_dir(), string(hash_a)))
        @test isfile(joinpath(OhMyArtifacts.get_artifacts_dir(), string(hash_b)))
        @test isfile(my_artifact_path(hash_a))
        @test isfile(my_artifact_path(hash_b))
        @test my_artifact_exists(hash_a)
        @test my_artifact_exists(hash_b)

        bind_my_artifact!(artifacts_toml, "runtestfile", hash_a)
        bind_my_artifact!(artifacts_toml, "utils.jl", hash_b)
        @test length(OhMyArtifacts.load_my_artifacts_toml(artifacts_toml)) == 2
        @test my_artifact_hash("runtestfile", artifacts_toml) == hash_a
        @test my_artifact_hash("utils.jl", artifacts_toml) == hash_b

        usagefile = OhMyArtifacts.usages_toml()
        @test isfile(usagefile)
        usage = OhMyArtifacts.parse_toml(usagefile)
        @test length(usage) == 2
        @test length(usage[my_artifact_path(hash_a)]) == 1
        @test length(usage[my_artifact_path(hash_b)]) == 1
        @test now() - usage[my_artifact_path(hash_a)][artifacts_toml]["runtestfile"] < Day(1)
        @test now() - usage[my_artifact_path(hash_b)][artifacts_toml]["utils.jl"] < Day(1)

        OhMyArtifacts.unbind_my_artifact!(artifacts_toml, "utils.jl")
        @test length(OhMyArtifacts.load_my_artifacts_toml(artifacts_toml)) == 1
        @test my_artifact_exists(hash_a)
        @test my_artifact_exists(hash_b)

        orphanfile = OhMyArtifacts.orphanages_toml()
        @test isfile(orphanfile)
        orphan = OhMyArtifacts.parse_toml(orphanfile)
        @test isempty(orphan)

        OhMyArtifacts.find_orphanages()
        orphan = OhMyArtifacts.parse_toml(orphanfile)
        @test length(orphan) == 1
        @test now() - orphan[my_artifact_path(hash_b)] < Day(1)
        @test my_artifact_exists(hash_a)
        @test my_artifact_exists(hash_b)

        OhMyArtifacts.find_orphanages(; collect_delay=Hour(0))
        orphan = OhMyArtifacts.parse_toml(orphanfile)
        @test isempty(orphan)
        @test my_artifact_exists(hash_a)
        @test !my_artifact_exists(hash_b)

        url = "https://raw.githubusercontent.com/chengchingwen/OhMyArtifacts.jl/master/README.md"
        hash_c = download_my_artifact!(url, "readme", artifacts_toml)
        @test my_artifact_hash("readme", artifacts_toml) == hash_c
        @test my_artifact_exists(hash_c)
        @test isfile(my_artifact_path(hash_c))

        another_toml = my_artifacts_toml!(Base.UUID("cf8be1f4-309d-442e-839d-29d2a0af6cb7"))
        @test isfile(another_toml)
        bind_my_artifact!(artifacts_toml, "runtests.jl", hash_a)
        bind_my_artifact!(another_toml, "runtests.jl", hash_a)
        usage = OhMyArtifacts.parse_toml(usagefile)
        @test length(usage) == 2
        @test length(usage[my_artifact_path(hash_a)]) == 2
        @test length(usage[my_artifact_path(hash_a)][artifacts_toml]) == 2
        @test length(usage[my_artifact_path(hash_c)]) == 1
    end
end


text:
"""
    SpinBasis(n)

Basis for spin-n particles.

The basis can be created for arbitrary spinnumbers by using a rational number,
e.g. `SpinBasis(3//2)`. The Pauli operators are defined for all possible
spin numbers.
"""
struct SpinBasis{S,T} <: Basis
    shape::Vector{T}
    spinnumber::Rational{T}
    function SpinBasis{S}(spinnumber::Rational{T}) where {S,T<:Int}
        @assert isa(S, Rational{Int})
        n = numerator(spinnumber)
        d = denominator(spinnumber)
        @assert d==2 || d==1
        @assert n > 0
        N = numerator(spinnumber*2 + 1)
        new{spinnumber,T}([N], spinnumber)
    end
end
SpinBasis(spinnumber::Rational{Int}) = SpinBasis{spinnumber}(spinnumber)
SpinBasis(spinnumber::Int) = SpinBasis(convert(Rational{Int}, spinnumber))

==(b1::SpinBasis, b2::SpinBasis) = b1.spinnumber==b2.spinnumber

"""
    sigmax(b::SpinBasis)

Pauli ``σ_x`` operator for the given Spin basis.
"""
function sigmax(b::SpinBasis)
    N = length(b)
    diag = ComplexF64[complex(sqrt(real((b.spinnumber + 1)*2*a - a*(a+1)))) for a=1:N-1]
    data = spdiagm(1 => diag, -1 => diag)
    SparseOperator(b, data)
end

"""
    sigmay(b::SpinBasis)

Pauli ``σ_y`` operator for the given Spin basis.
"""
function sigmay(b::SpinBasis)
    N = length(b)
    diag = ComplexF64[1im*complex(sqrt(real((b.spinnumber + 1)*2*a - a*(a+1)))) for a=1:N-1]
    data = spdiagm(-1 => diag, 1 => -diag)
    SparseOperator(b, data)
end

"""
    sigmaz(b::SpinBasis)

Pauli ``σ_z`` operator for the given Spin basis.
"""
function sigmaz(b::SpinBasis)
    N = length(b)
    diag = ComplexF64[complex(2*m) for m=b.spinnumber:-1:-b.spinnumber]
    data = spdiagm(0 => diag)
    SparseOperator(b, data)
end

"""
    sigmap(b::SpinBasis)

Raising operator ``σ_+`` for the given Spin basis.
"""
function sigmap(b::SpinBasis)
    N = length(b)
    S = (b.spinnumber + 1)*b.spinnumber
    diag = ComplexF64[complex(sqrt(float(S - m*(m+1)))) for m=b.spinnumber-1:-1:-b.spinnumber]
    data = spdiagm(1 => diag)
    SparseOperator(b, data)
end

"""
    sigmam(b::SpinBasis)

Lowering operator ``σ_-`` for the given Spin basis.
"""
function sigmam(b::SpinBasis)
    N = length(b)
    S = (b.spinnumber + 1)*b.spinnumber
    diag = [complex(sqrt(float(S - m*(m-1)))) for m=b.spinnumber:-1:-b.spinnumber+1]
    data = spdiagm(-1 => diag)
    SparseOperator(b, data)
end

"""
    spinup(b::SpinBasis)

Spin up state for the given Spin basis.
"""
spinup(b::SpinBasis) = basisstate(b, 1)

"""
    spindown(b::SpinBasis)

Spin down state for the given Spin basis.
"""
spindown(b::SpinBasis) = basisstate(b, b.shape[1])


text:
import Radio
const rad = Radio

( N, β )        = rad.kaiserord( 0.001, 0.2*π )
window          = rad.kaiser( N, β )
impulseResponse = rad.firdes( 0.5, window )
p               = rad.freqz( impulseResponse )

display( p )
file( p, "Kaiser.png" )


text:
using GeometricIntegrators
using GeometricProblems.KuboOscillator
using Test

@testset "$(rpad("Kubo Oscillator",80))" begin

    sde   = kubo_oscillator_sde_1()
    psde  = kubo_oscillator_psde_1()
    spsde = kubo_oscillator_spsde_1()

    sde_equs = get_functions(sde)
    @test_nowarn sde_equs[:v](sde.t₀, sde.q₀[begin], zero(sde.q₀[begin]))
    @test_nowarn sde_equs[:B](sde.t₀, sde.q₀[begin], zeros(eltype(sde.q₀[begin]), sde.d, sde.m))
    @test_nowarn sde_equs[:B](sde.t₀, sde.q₀[begin], zeros(eltype(sde.q₀[begin]), sde.d, sde.m), 1)

    psde_equs = get_functions(psde)
    @test_nowarn psde_equs[:v](psde.t₀, psde.q₀[begin], psde.p₀[begin], zero(psde.q₀[begin]))
    @test_nowarn psde_equs[:f](psde.t₀, psde.q₀[begin], psde.p₀[begin], zero(psde.p₀[begin]))
    @test_nowarn psde_equs[:B](psde.t₀, psde.q₀[begin], psde.p₀[begin], zero(psde.q₀[begin]))
    @test_nowarn psde_equs[:G](psde.t₀, psde.q₀[begin], psde.p₀[begin], zero(psde.p₀[begin]))
    @test_nowarn psde_equs[:B](psde.t₀, psde.q₀[begin], psde.p₀[begin], zeros(eltype(psde.q₀[begin]), psde.d, psde.m))
    @test_nowarn psde_equs[:G](psde.t₀, psde.q₀[begin], psde.p₀[begin], zeros(eltype(psde.p₀[begin]), psde.d, psde.m))

    spsde_equs = get_functions(spsde)
    @test_nowarn spsde_equs[:v ](spsde.t₀, spsde.q₀[begin], spsde.p₀[begin], zero(spsde.q₀[begin]))
    @test_nowarn spsde_equs[:f1](spsde.t₀, spsde.q₀[begin], spsde.p₀[begin], zero(spsde.p₀[begin]))
    @test_nowarn spsde_equs[:f2](spsde.t₀, spsde.q₀[begin], spsde.p₀[begin], zero(spsde.p₀[begin]))
    @test_nowarn spsde_equs[:B ](spsde.t₀, spsde.q₀[begin], spsde.p₀[begin], zero(spsde.q₀[begin]))
    @test_nowarn spsde_equs[:G1](spsde.t₀, spsde.q₀[begin], spsde.p₀[begin], zero(spsde.p₀[begin]))
    @test_nowarn spsde_equs[:G2](spsde.t₀, spsde.q₀[begin], spsde.p₀[begin], zero(spsde.p₀[begin]))
    @test_nowarn spsde_equs[:B ](spsde.t₀, spsde.q₀[begin], spsde.p₀[begin], zeros(eltype(spsde.q₀[begin]), spsde.d, spsde.m))
    @test_nowarn spsde_equs[:G1](spsde.t₀, spsde.q₀[begin], spsde.p₀[begin], zeros(eltype(spsde.p₀[begin]), spsde.d, spsde.m))
    @test_nowarn spsde_equs[:G2](spsde.t₀, spsde.q₀[begin], spsde.p₀[begin], zeros(eltype(spsde.p₀[begin]), spsde.d, spsde.m))

end


text: println("1 + 2 = $(1 + 2)")

max_stars_repo_name: bluey-crypto/Hacktoberfest

max_issues_repo_name: bluey-crypto/Hacktoberfest

max_forks_repo_name: bluey-crypto/Hacktoberfest


text:
@show ARGS
@show ARGS[1]
@show ARGS[2]


text:
#################################################################################
#Type Alias for long type signatures
const MinMax = NamedTuple{(:min, :max), NTuple{2, Float64}}
const NamedMinMax = Tuple{String, MinMax}
const UpDown = NamedTuple{(:up, :down), NTuple{2, Float64}}
const InOut = NamedTuple{(:in, :out), NTuple{2, Float64}}

# Type Alias From other Packages
const PM = PowerModels
const PSY = PowerSystems
const PSI = PowerSimulations
const IS = InfrastructureSystems
const MOI = MathOptInterface
const MOIU = MathOptInterface.Utilities
const PJ = ParameterJuMP
const MOPFM = MathOptFormat.MOF.Model
const TS = TimeSeries

#Type Alias for JuMP and PJ containers
const JuMPExpressionMatrix = Matrix{<:JuMP.AbstractJuMPScalar}
const PGAE{V} = PJ.ParametrizedGenericAffExpr{Float64, V} where V<:JuMP.AbstractVariableRef
const GAE{V} = JuMP.GenericAffExpr{Float64, V} where V<:JuMP.AbstractVariableRef
const JuMPAffineExpressionArray = Matrix{GAE{V}} where V<:JuMP.AbstractVariableRef
const JuMPConstraintArray = JuMP.Containers.DenseAxisArray{JuMP.ConstraintRef}
const JuMPParamArray = JuMP.Containers.DenseAxisArray{PJ.ParameterRef}
const DSDA = Dict{Symbol, JuMP.Containers.DenseAxisArray}


text:
using FastBEAST
using StaticArrays
using LinearAlgebra

# Do a 3D test with Laplace kernel

function OneoverRkernel(testpoint::SVector{3,T}, sourcepoint::SVector{3,T}) where T
    if isapprox(testpoint, sourcepoint, rtol=eps()*1e-4)
        return 0.0
    else
        return 1.0 / (norm(testpoint - sourcepoint))
    end
end

function assembler(kernel, testpoints, sourcepoints)
    kernelmatrix = zeros(promote_type(eltype(testpoints[1]),eltype(sourcepoints[1])),
                length(testpoints), length(sourcepoints))

    for i = 1:length(testpoints)
        for j = 1:length(sourcepoints)
            kernelmatrix[i,j] = kernel(testpoints[i], sourcepoints[j])
        end
    end
    return kernelmatrix
end

function assembler(kernel, matrix, testpoints, sourcepoints)
    for i = 1:length(testpoints)
        for j = 1:length(sourcepoints)
            matrix[i,j] = kernel(testpoints[i], sourcepoints[j])
        end
    end
end

##

N =  1000
NT = N

spoints = [@SVector rand(3) for i = 1:N]
tpoints = 0.1*[@SVector rand(3) for i = 1:NT] + [1.0*SVector(3.5, 3.5, 3.5) for i = 1:NT]

@views OneoverRkernelassembler(matrix, tdata, sdata) = assembler(OneoverRkernel, matrix, tpoints[tdata], spoints[sdata])
stree = create_tree(spoints, BoxTreeOptions(nmin=50))
ttree = create_tree(tpoints, BoxTreeOptions(nmin=50))
kmat = assembler(OneoverRkernel, tpoints, spoints)
hmat = HMatrix(OneoverRkernelassembler, ttree, stree, compressor=:naive, T=Float64)

@test estimate_reldifference(hmat,kmat) ≈ 0 atol=1e-4
@test compressionrate(hmat)*100 ≈ 99 atol=1

##
N = 4000
NT = N

spoints = [@SVector rand(3) for i = 1:N]
##
@views OneoverRkernelassembler(matrix, tdata, sdata) = assembler(OneoverRkernel, matrix, spoints[tdata], spoints[sdata])
stree = create_tree(spoints, BoxTreeOptions(nmin=400))
kmat = assembler(OneoverRkernel, spoints, spoints)
hmat = HMatrix(OneoverRkernelassembler, stree, stree, compressor=:aca, T=Float64)

@test estimate_reldifference(hmat,kmat) ≈ 0 atol=1e-4
@test compressionrate(hmat)*100 ≈ 54 atol=1

max_stars_repo_name: sbadrian/FastBEAST

max_issues_repo_name: sbadrian/FastBEAST

max_forks_repo_name: sbadrian/FastBEAST


text:
module QuadDIRECT

using LinearAlgebra
using StaticArrays, PositiveFactorizations

export Box, CountedFunction, LoggedFunction
export leaves, value, numevals, analyze, analyze!, minimize
export splitprint, splitprint_colored, treeprint

include("types.jl")
include("util.jl")
include("quadratic_model.jl")
include("algorithm.jl")

# Convenience methods
function analyze(f, x0::AbstractVector{<:Real}, lower::AbstractVector{<:Real}, upper::AbstractVector{<:Real}; kwargs...)
    inds = LinearIndices(x0)
    LinearIndices(lower) == LinearIndices(upper) == inds || error("lengths must match")
    splits = similar(x0, Vector{Float64})
    for i in inds
        lo, hi, xi = lower[i], upper[i], x0[i]
        @assert(lo<=xi<=hi)
        lo = max(xi-1, (xi+lo)/2)
        hi = min(xi+1, (xi+hi)/2)
        xi = xi == lo ? (xi+hi)/2 :
             xi == hi ? (xi+lo)/2 : xi
        @assert(lo < xi < hi)
        @assert isfinite(lo) && isfinite(xi) && isfinite(hi)
        splits[i] = [lo, xi, hi]
    end
    return analyze(f, splits, lower, upper; kwargs...)
end

function analyze(f, x0::AbstractVector{<:Real}; kwargs...)
    return analyze(f, x0, fill(-Inf, length(x0)), fill(Inf, length(x0)); kwargs...)
end

function analyze(f, lower::AbstractVector{<:Real}, upper::AbstractVector{<:Real}; kwargs...)
    x0 = (lower .+ upper)./2
    for i in eachindex(x0)
        if !isfinite(x0[i])
            x0[i] = clamp(0, lower[i], upper[i])
        end
    end
    return analyze(f, x0, lower, upper; kwargs...)
end

end # module


text:
mutable struct TD; ptr
    function TD(a::CuArray)
        d = Cptr[0]
        @cuda(cudnn,cudnnCreateTensorDescriptor,(Ptr{Cptr},),d)
        n = ndims(a)
        sz = [Cint(size(a,n-i+1)) for i=1:n]
        st = [Cint(stride(a,n-i+1)) for i=1:n]
        @cuda(cudnn,cudnnSetTensorNdDescriptor,
              (Cptr,UInt32,Cint,Ptr{Cint},Ptr{Cint}),
              d[1], DT(a), n, sz, st)
        td = new(d[1])
        finalizer(td, x->@cuda(cudnn,cudnnDestroyTensorDescriptor,(Cptr,),x.ptr))
        return td
    end
end

mutable struct FD; ptr
    function FD(a::CuArray)
        d = Cptr[0]
        @cuda(cudnn,cudnnCreateFilterDescriptor,(Ptr{Cptr},),d)
        n = ndims(a)
        sz = [Cint(size(a,n-i+1)) for i=1:n]
        cudnnVersion = cudnn_version()
        if cudnnVersion >= 5000
            @cuda(cudnn,cudnnSetFilterNdDescriptor,
                  (Cptr,UInt32,UInt32,Cint,Ptr{Cint}),
                  d[1], DT(a), 0,     n,   sz)
        elseif cudnnVersion >= 4000
            @cuda(cudnn,cudnnSetFilterNdDescriptor_v4,
                  (Cptr,UInt32,UInt32,Cint,Ptr{Cint}),
                  d[1], DT(a), 0,     n,   sz)
        else
            @cuda(cudnn,cudnnSetFilterNdDescriptor,
                  (Cptr,UInt32,Cint,Ptr{Cint}),
                  d[1], DT(a),    n,   sz)
        end
        fd = new(d[1])
        finalizer(fd, x->@cuda(cudnn,cudnnDestroyFilterDescriptor,(Cptr,),x.ptr))
        return fd
    end
end

mutable struct CD; ptr
    function CD(w::CuArray,x::CuArray; padding=0, stride=1, upscale=1, mode=0)
        d = Cptr[0]
        @cuda(cudnn,cudnnCreateConvolutionDescriptor,(Ptr{Cptr},),d)
        nd = ndims(x)-2
        cudnnVersion = cudnn_version()
        if cudnnVersion >= 4000
            @cuda(cudnn,cudnnSetConvolutionNdDescriptor,
                  (Cptr,Cint,Ptr{Cint},Ptr{Cint},Ptr{Cint},UInt32,UInt32),
                  d[1],nd,cdsize(padding,nd),cdsize(stride,nd),cdsize(upscale,nd),mode,DT(x))
        elseif cudnnVersion > 3000 # does not work when cudnnVersion==3000
            @cuda(cudnn,cudnnSetConvolutionNdDescriptor_v3,
                  (Cptr,Cint,Ptr{Cint},Ptr{Cint},Ptr{Cint},UInt32,UInt32),
                  d[1],nd,cdsize(padding,nd),cdsize(stride,nd),cdsize(upscale,nd),mode,DT(x))
        else
            @cuda(cudnn,cudnnSetConvolutionNdDescriptor,
                  (Cptr,Cint,Ptr{Cint},Ptr{Cint},Ptr{Cint},UInt32),
                  d[1],nd,cdsize(padding,nd),cdsize(stride,nd),cdsize(upscale,nd),mode)
        end
        cd = new(d[1])
        finalizer(cd, x->@cuda(cudnn,cudnnDestroyConvolutionDescriptor,(Cptr,),x.ptr))
        return cd
    end
end

mutable struct PD; ptr
    function PD(x::CuArray; window=2, padding=0, stride=window, mode=0, maxpoolingNanOpt=0)
        d = Cptr[0]
        @cuda(cudnn,cudnnCreatePoolingDescriptor,(Ptr{Cptr},),d)
        nd = ndims(x)-2
        cudnnVersion = cudnn_version()
        if cudnnVersion >= 5000
            @cuda(cudnn,cudnnSetPoolingNdDescriptor,
                  (Cptr,UInt32,UInt32,Cint,Ptr{Cint},Ptr{Cint},Ptr{Cint}),
                  d[1],mode,maxpoolingNanOpt,nd,cdsize(window,nd),cdsize(padding,nd),cdsize(stride,nd))
        elseif cudnnVersion >= 4000
            @cuda(cudnn,cudnnSetPoolingNdDescriptor_v4,
                  (Cptr,UInt32,UInt32,Cint,Ptr{Cint},Ptr{Cint},Ptr{Cint}),
                  d[1],mode,maxpoolingNanOpt,nd,cdsize(window,nd),cdsize(padding,nd),cdsize(stride,nd))
        else
            @cuda(cudnn,cudnnSetPoolingNdDescriptor,
                  (Cptr,UInt32,Cint,Ptr{Cint},Ptr{Cint},Ptr{Cint}),
                  d[1],mode,nd,cdsize(window,nd),cdsize(padding,nd),cdsize(stride,nd))
        end
        pd = new(d[1])
        finalizer(pd, x->@cuda(cudnn,cudnnDestroyPoolingDescriptor,(Cptr,),x.ptr))
        return pd
    end
end

# mutable struct RNN_D; ptr
#     function RNN_D(h_size::Int, num_layers::Int; handle=cudnnhandle())
#         d = Cptr[0]
#         @cuda(cudnn, cudnnCreateRNNDescriptor,(Ptr{Cptr},),d)

#         @cuda(cudnn, cudnnSetRNNDescriptor,(Ptr{Cptr},),d)

#         rnn_d = new(d[1])
#         finalizer(rnn_d, x->@cuda(cudnn,cudnnDestroyRNNDescriptor,(Cptr,),x.ptr))
#         return rnn_d
#     end
# end

import Base: unsafe_convert
unsafe_convert(::Type{Cptr}, td::TD)=td.ptr
unsafe_convert(::Type{Cptr}, fd::FD)=fd.ptr
unsafe_convert(::Type{Cptr}, cd::CD)=cd.ptr
unsafe_convert(::Type{Cptr}, pd::PD)=pd.ptr


text:
# Try to recreate the example from Wuestefeld & Bokelmann, BSSA, 2007,
# Figure 1, and use this to both profile the splitting calculation
# overall, and also verify the rotation-correlation analysis.

using Revise, Seis, SeisSplit, Plots
import Juno, DelimitedFiles

# Synthetic parameters
noise = √(1/15)
spol = 60
freq = 1/8
phi_true = 0
dt_true = 1.3
dt_max = 3

# Create synthetics using programs from
#   https://github.com/anowacki/seismo-fortran
run(`create_wave -n $noise -f $freq -spol $spol`)
run(`splitwave $phi_true $dt_true`)

# Filter and save
t = read_sac("wave.BH[ENZ]", ".", echo=false)
bandpass!.(t, 0.02, 1, twopass=true)
write_sac.(t, t.meta.file)

# Calculate splitting using SHEBA for reference
open("sheba.in", "w") do f
    println(f, """
        # sheba.in
        wave
        BHE
        BHN
        BHZ
        1
        1 1
        $dt_max
        0
        0
        """)
end
withenv(()->run(`sheba_exec`), "SHEBA_WRITE_XC"=>"1")
sheba_phi, sheba_dphi, sheba_dt, sheba_ddt =
    parse.(Float64, split(readlines("wave_sheba.result")[2])[11:14])
sheba_Q = parse(Float64, split(readlines("wave_sheba.stats")[2])[5])
sheba_XC = DelimitedFiles.readdlm("sheba.xc")
sheba_phi_XC, sheba_dt_XC = parse.(Float64,
    split(readlines("sheba.xc_best")[2])[[1,3]])
sheba_phis_XC = -90:90
sheba_dts_XC = range(0, step=2t[1].delta, length=size(sheba_XC, 2))

# Calculate splitting using SeisSplit
e, n = Seis.read_sac("wave.BH[EN]", ".", echo=false)
@time s = splitting(n, e, e.picks.A.time, e.picks.F.time, dt_max=dt_max)

# Uncomment for profiling display with Juno
# Juno.@profiler begin
#     a = 0.0
#     for _ in 1:20
#         t, phi, xc, phi_best, dt_best = SeisSplit.rotation_correlation(n, e,
#             e.picks.A.time, e.picks.F.time)
#         a += t[1]
#     end
#     a
# end

p = plot(layout=(2,1))

# Xcorr heatmap, λ₂ contours, and best results
heatmap!(p[1], s.xcorr_dt, s.xcorr_phi, s.xcorr_map', c=:RdBu,
    xlabel="\\delta t / s", ylabel="\\phi / °",
    framestyle=:box, xlim=(0,dt_max), ylim=(-90,90))
contour!(p[1], s.dt, s.phi, s.lam2.*maximum(s.xcorr_map)./maximum(s.lam2))
scatter!(p[1], [s.xcorr_dt_best], [s.xcorr_phi_best], label="XC", markersize=15)
scatter!(p[1], [sheba_dt_XC], [sheba_phi_XC], label="SHEBA XC", markersize=12)
scatter!(p[1], [s.dt_best], [s.phi_best], label="EV", markersize=10)
scatter!(p[1], [sheba_dt], [sheba_phi], label="SHEBA EV", markersize=6)
scatter!(p[1], [dt_true], [phi_true], label="True", markersize=4)

# Xcorr surfaces, all normalised
XC_sheba_norm = sheba_XC ./ maximum(abs, sheba_XC)
XC_seissplit_norm = s.xcorr_map./maximum(abs, s.xcorr_map)
XC_splitlab_norm = DelimitedFiles.readdlm("splitlab.xc")
XC_splitlab_norm ./= maximum(abs, XC_splitlab_norm)
XC_splitlab_norm = reverse(XC_splitlab_norm, dims=1)
splitlab_phis = -90:89
splitlab_dts = range(0, dt_max, length=size(XC_splitlab_norm, 2))
iphi, idt = Tuple(argmin(XC_splitlab_norm))
splitlab_phi = splitlab_phis[iphi]
splitlab_dt = splitlab_dts[idt]

contour!(p[2], sheba_dts_XC, sheba_phis_XC, XC_sheba_norm,
    l=:red, levels=-1:0.2:1, label="SHEBA", xlim=(0,dt_max))
contour!(p[2], s.xcorr_dt, s.xcorr_phi, XC_seissplit_norm',
    l=:black, levels=-1:0.2:1, label="SeisSplit")
contour!(p[2], splitlab_dts, splitlab_phis, XC_splitlab_norm,
    l=:blue, levels=-1:0.2:1, label="Splitlab")
scatter!(p[2], [sheba_dt s.xcorr_dt_best splitlab_dt],
    [sheba_phi s.xcorr_phi_best splitlab_phi], markercolor=[:red :black :blue],
    markersize=[10 8 6],
    label=["SHEBA" "SeisSplit" "Splitlab"])

@show s.xcorr_phi_best, s.xcorr_dt_best
@show sheba_phi_XC, sheba_dt_XC
@show SeisSplit.quality(s), sheba_Q

p # show plot

max_stars_repo_name: anowacki/SeisSplit

max_issues_repo_name: anowacki/SeisSplit

max_forks_repo_name: anowacki/SeisSplit


text:
#######################################
## Functions creating visualizations ##
#######################################

function visualize_RSUs_and_failures(OSMmap::MapData,
                                    Start::Vector{Rect},
                                    End::Vector{Rect},
                                    Agents::Vector{Agent},
                                    failedENU::Array{Array{ENU,1},1},
                                    RSUs::Array{RSU,1},
                                    range::Float64,
                                    outfile::String)
    RSU_ENU = getfield.(RSUs,:ENU)
    flm = pyimport("folium")
    matplotlib_cm = pyimport("matplotlib.cm")
    matplotlib_colors = pyimport("matplotlib.colors")
    cmap = matplotlib_cm.get_cmap("prism")

    m = flm.Map()
    #Add RSUs location and their range
    locs = [OpenStreetMapX.LLA(n, OSMmap.bounds) for n in RSU_ENU]
    for k = 1:length(RSU_ENU)
        info = "RSU number: $k"
        flm.Circle(
          location=[locs[k].lat,locs[k].lon],
          popup=info,
          tooltip=info,
          radius=range,
          color="blue",
          weight=0.5,
          fill=true,
          fill_color="blue"
       ).add_to(m)
            flm.Circle(
          location=[locs[k].lat,locs[k].lon],
          popup=info,
          tooltip=info,
          radius=1,
          color="crimson",
          weight=3,
          fill=false,
          fill_color="crimson"
       ).add_to(m)
    end

    #Add subset of agents' routes
    routes = [a.route for a in Agents if a.smart]
    #Subset routes if more than 500 smart agents - HTML performance
    if length(routes) > 500 routes = routes[1:500] end
    for z = 1:length(routes)
        locs = [OpenStreetMapX.LLA(OSMmap.nodes[n], OSMmap.bounds) for n in routes[z]]
        info = "Agent $z route\n
"*
            "Length: $(length(routes[z])) nodes\n
" *
            "From: $(routes[z][1]) $(round.((locs[1].lat, locs[1].lon),digits=4))\n
" *
            "To: $(routes[z][end]) $(round.((locs[end].lat, locs[end].lon),digits=4))"
        flm.PolyLine(
            [(loc.lat, loc.lon) for loc in locs ],
            popup=info,
            tooltip=info,
            weight = 2,
            color="green"
        ).add_to(m)
    end

    #Add points of failed transmission
    failedENU = unique(collect(Iterators.flatten(failedENU)))
    flocs = [OpenStreetMapX.LLA(n, OSMmap.bounds) for n in failedENU]
    for i = 1:length(failedENU)
        info = "Failed ID: $i"
        flm.Circle(
          location=[flocs[i].lat, flocs[i].lon],
          popup=info,
          tooltip=info,
          radius=8,
          color="black",
          weight=3,
          fill=true,
            fill_color="black"
       ).add_to(m)
    end

    #Add bounds and starting/ending area
    MAP_BOUNDS = [(OSMmap.bounds.min_y,OSMmap.bounds.min_x),(OSMmap.bounds.max_y,OSMmap.bounds.max_x)]
    flm.Rectangle(MAP_BOUNDS, color="black",weight=6).add_to(m)
    for r in Start
      flm.Rectangle([r.p1,r.p2], color="blue",weight=2).add_to(m)
    end
    for r in End
      flm.Rectangle([r.p1,r.p2], color="red",weight=2).add_to(m)
    end
    m.fit_bounds(MAP_BOUNDS)
    m.save(outfile)
end

function visualize_bounds(OSMmap::MapData,
                          Start::Vector{Rect},
                          End::Vector{Rect},
                          outfile::String)
    flm = pyimport("folium")
    matplotlib_cm = pyimport("matplotlib.cm")
    matplotlib_colors = pyimport("matplotlib.colors")
    cmap = matplotlib_cm.get_cmap("prism")

    m = flm.Map()
    #Add bounds and starting/ending area
    MAP_BOUNDS = [(OSMmap.bounds.min_y,OSMmap.bounds.min_x),(OSMmap.bounds.max_y,OSMmap.bounds.max_x)]
    flm.Rectangle(MAP_BOUNDS, color="black", weight=6).add_to(m)
    for r in Start
      flm.Rectangle([r.p1,r.p2], color="blue", weight=2).add_to(m)
    end
    for r in End
      flm.Rectangle([r.p1,r.p2], color="red", weight=2).add_to(m)
    end
    m.fit_bounds(MAP_BOUNDS)
    m.save(outfile)
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

## dummy stub for https://github.com/JuliaBinaryWrappers/GMP_jll.jl
baremodule GMP_jll
using Base, Libdl
Base.Experimental.@compiler_options compile=min optimize=0 infer=false

const PATH_list = String[]
const LIBPATH_list = String[]

export libgmp, libgmpxx

# These get calculated in __init__()
const PATH = Ref("")
const LIBPATH = Ref("")
artifact_dir = ""
libgmp_handle = C_NULL
libgmp_path = ""
libgmpxx_handle = C_NULL
libgmpxx_path = ""

if Sys.iswindows()
    const libgmp = "libgmp-10.dll"
    const libgmpxx = "libgmpxx-4.dll"
elseif Sys.isapple()
    const libgmp = "@rpath/libgmp.10.dylib"
    const libgmpxx = "@rpath/libgmpxx.4.dylib"
else
    const libgmp = "libgmp.so.10"
    const libgmpxx = "libgmpxx.so.4"
end

function __init__()
    global libgmp_handle = dlopen(libgmp)
    global libgmp_path = dlpath(libgmp_handle)
    global libgmpxx_handle = dlopen(libgmpxx)
    global libgmpxx_path = dlpath(libgmpxx_handle)
    global artifact_dir = dirname(Sys.BINDIR)
    LIBPATH[] = dirname(libgmp_path)
    push!(LIBPATH_list, LIBPATH[])
end

# JLLWrappers API compatibility shims.  Note that not all of these will really make sense.
# For instance, `find_artifact_dir()` won't actually be the artifact directory, because
# there isn't one.  It instead returns the overall Julia prefix.
is_available() = true
find_artifact_dir() = artifact_dir
dev_jll() = error("stdlib JLLs cannot be dev'ed")
best_wrapper = nothing
get_libgmp_path() = libgmp_path
get_libgmpxx_path() = libgmpxx_path

end  # module GMP_jll

max_stars_repo_name: jonas-schulze/julia

max_issues_repo_name: jonas-schulze/julia

max_forks_repo_name: jonas-schulze/julia


text:
using ScikitLearn.CrossValidation: cross_val_score
using ScikitLearn
using Statistics
@sk_import naive_bayes: GaussianNB

function ec_ranking(data, target, weights, mu_vals)

	# Get maximal ReliefF weight and compute epsilon values.
	max_weight = maximum(weights)
	eps = (max_weight .- weights)./max_weight

	# Set initial tinfo value.
	tinfo = 1

	# Initialize vector of feature ranks.
	rank = Array{Int64}(undef, length(weights))

	# Initialize vector of indices in original weights array.
	index_vector = collect(1:length(weights))

	# Initialize initial rank value.
	rank_value_next = length(weights)

	# Initialize variable that holds current best tinfo value.
	best_tinfo = tinfo

	# While there are unranked features, perform evaporation.
	while length(index_vector) > 1

		### Grid search for best temperature from previous step ###

		# Initialize variables that hold current maximal CV score
		# and index of removed feature.
		max_cv_val = 0.0
		idx_removed_feature = -1

		# Perform grid search for best value of tinfo.
		for tinfo_nxt = tinfo-0.3:0.1:tinfo+0.3

			# Compute weights for next value of tinfo.
			weights_nxt = eps .- tinfo_nxt.*mu_vals

			# Rank features.
			enumerated_weights = [weights_nxt collect(1:length(weights_nxt))]'
			rank_weights = zeros(Int64, length(weights_nxt))
			s = enumerated_weights[:, sortperm(enumerated_weights[1,:], rev=false)]
			rank_weights[convert.(Int64, s[2, :])] = 1:length(weights_nxt)

			# Remove lowest ranked feature.
			msk_rm = rank_weights .!= length(weights_nxt)

			# Perform 5-fold cross validation.
			data_filt = data[:, msk_rm]

			# Compare to current maximal CV value.
			cv_val_nxt = Statistics.mean(cross_val_score(GaussianNB(), data_filt, target; cv=5))

			# If CV value greater than current maximal, save current
			# CV score, tinfo and feature weights.
			if cv_val_nxt > max_cv_val
				# Set current tinfo value as best.
				best_tinfo = tinfo_nxt
				# Get index of removed feature.
				idx_removed_feature = findall(.!msk_rm)[1]
			end
		end

		# Remove evaporated feature from data matrix,
		# and data at corresponding index from eps vector and
		# mu values vector.
		data = data[:, 1:end .!= idx_removed_feature]
		eps = eps[1:end .!= idx_removed_feature]
		mu_vals = mu_vals[1:end .!= idx_removed_feature]

		# Get index of evaporated feature in original data matrix.
		rank_index_next = index_vector[idx_removed_feature]

		# Delete removed feature from vector of indices.
		index_vector = deleteat!(index_vector, idx_removed_feature)

		# Assign rank to removed feature.
		rank[rank_index_next] = rank_value_next

		# Decrease next rank value.
		rank_value_next -= 1

	end

	# Assign final rank.
	rank[index_vector[1]] = rank_value_next

	# Return vector of ranks.
	return rank

end

max_stars_repo_name: Architecton/diplomska-naloga-koda

max_issues_repo_name: Architecton/diplomska-naloga-koda

max_forks_repo_name: Architecton/diplomska-naloga-koda


text:
module TestMutable

using UnsafeFields.Demos: Mutable, freeze
using Test

struct Values2{A,B}
    a::A
    b::B
end

#=
struct Values3{A,B,C}
    a::A
    b::B
    c::C
end
=#

@testset "Tuple{Int, Int}" begin
    mut = Mutable((111, 222))
    @test (mut[1], mut[2]) == freeze(mut) == (111, 222)
    mut[1] = 3
    mut[2] = 4
    @test (mut[1], mut[2]) == freeze(mut) == (3, 4)
end

@testset "Values2{String,Symbol}" begin
    mut = Mutable(Values2("initial string", :initial_symbol))
    @test (mut.a, mut.b) == ("initial string", :initial_symbol)
    mut.a = "second string"
    mut.b = :second_symbol
    @test (mut.a, mut.b) == ("second string", :second_symbol)
end

@testset "Values2{Union{String,Int,Nothing},Union{Int,Nothing,Float64}}" begin
    mut = Mutable(
        Values2{Union{String,Int,Nothing},Union{Int,Nothing,Float64}}(nothing, nothing),
    )
    @test (mut.a, mut.b) === (nothing, nothing)
    mut.a = "second string"
    mut.b = 222
    @test (mut.a, mut.b) === ("second string", 222)
    mut.a = 111
    mut.b = 333.0
    @test (mut.a, mut.b) === (111, 333.0)
end

@testset "Values2{Any,Union{Integer,Symbol}}" begin
    mut = Mutable(Values2{Any,Union{Integer,Symbol}}(0, 0))
    @test (mut.a, mut.b) === (0, 0)
    mut.a = "second string"
    mut.b = :second_symbol
    @test (mut.a, mut.b) == ("second string", :second_symbol)
end

end  # module


text:
n = 5
m = 2

n + n^2.0 + m + m^2.0

2 * n^2.0 + 2 * m * n

max_stars_repo_name: zacmanchester/DirectMotionPlanning

max_issues_repo_name: zacmanchester/DirectMotionPlanning

max_forks_repo_name: zacmanchester/DirectMotionPlanning


text:
rel_error(x, y) = (x - y) / y
function check_gradient(data, a = randn(length(data)))
    acopy = copy(a)
    lp, g = logdensity_and_gradient(data, a)
    all(i -> a[i] == acopy[i], eachindex(a)) || throw("Logdensity mutated inputs!?!?!?")
    for i ∈ eachindex(a)
        aᵢ = a[i]
        step = cbrt(eps(aᵢ))
        a[i] = aᵢ + step
        lp_hi = logdensity(data, a)
        a[i] = aᵢ - step
        lp_lo = logdensity(data, a)
        a[i] = aᵢ
        fd = (lp_hi - lp_lo) / (2step)
        ad = g[i]
        relative_error = rel_error(ad, fd)
        @show (i, ad, fd, relative_error)
        if abs(relative_error) > 1e-5
            fd_f = (lp_hi - lp) / step
            fd_b = (lp - lp_lo) / step
            @show rel_error.(ad, (fd_f, fd_b))
        end
    end
end


text:
module Builder

export execute

function execute(args::Dict{Symbol,Any})
  println("Building...")
end

end # module


text:
using FantasyFootball

tests = [
  "scoring",
  "season",
  "stat"
]

for t in tests
  tfile = string(t, ".jl")
  print(" * $(tfile) ... ")
  include(tfile)
  println("ok")
end

max_stars_repo_name: jnu/FantasyFootball

max_issues_repo_name: jnu/FantasyFootball

max_forks_repo_name: jnu/FantasyFootball


text:
########################## Auxiliary functions #############################

function eye(tpe::Type,n::S) where {S <: Integer}

  if tpe <: Number
      return Matrix{tpe}(I,n,n)
  end

end

eye(n::Integer) = eye(Float64,n::Integer)

function tracem(x::Array{T,2}) where {T <: Real}

# Computes the matrix-trace as defined by Gomme and Klein (2011)

# We require the number of rows to be greater than the number of columns, so
# that m is greater than one.

trans = false

if size(x,1) < size(x,2)
  x = Matrix(x')
  trans = true
end

n = size(x,2)
m = Int(size(x,1)/n)

y = zeros(m,1)  # We want this to be a 2-d array for subsequent matrix multiplication

@inbounds for i = 1:m
  @views y[i,1] = tr(x[(n*(i-1)+1):i*n,1:n])
end

if trans == true
  y = y'
end

return Matrix(y)

end

function kron_prod_times_vector(a::AbstractArray{T,2},b::AbstractArray{T,2},v::AbstractArray{T,1}) where {T <: Real}

  # This function efficiently computes kron(a,b)*v where a and b are matrices
  # and v is a vector.

  (n1,n2) = size(a)
  (n3,n4) = size(b)
  n5      = length(v)

  if n5 != n2*n4
      error("a, b, and v do not have the correct sizes")
  end

  p = vec(b*reshape(v,n4,n2)*a')

  return p

end

function vector_times_kron_prod(v::AbstractArray{T,1},a::AbstractArray{T,2},b::AbstractArray{T,2}) where {T <: Real}

  product = Matrix(kron_prod_times_vector(a',b',v')')

  return product

end

function kron_prod_times_matrix(a::AbstractArray{T,2},b::AbstractArray{T,2},v::AbstractArray{T,2}) where {T <: Real}

  # This function efficiently computes kron(a,b)*v where a, b, and v are
  # conformable matrices.

  (n1,n2) = size(a)
  (n3,n4) = size(b)
  (n5,n6) = size(v)

  if n5 != n2*n4
      error("a, b, and v do not have the correct sizes")
  end

  p = zeros(n1*n3,n6)
  for i = 1:n6
      @views p[:,i] = kron_prod_times_vector(a,b,v[:,i])
  end

  return p

end

function matrix_times_kron_prod(v::AbstractArray{T,2},a::AbstractArray{T,2},b::AbstractArray{T,2}) where {T <: Real}

  product = Matrix(kron_prod_times_matrix(a',b',v')')

  return product

end

function kron_prod_times_vector(a::AbstractArray{T,2},b::AbstractArray{T,2},c::AbstractArray{T,2},v::AbstractArray{T,1}) where {T <: Real}

  # This function efficiently computes kron(a,kron(b,c))*v where a, b, and c
  # are matrices and v is a vector.  Included for completeness, but not
  # actually used.

  (n1,n2) = size(a)
  (n3,n4) = size(b)
  (n5,n6) = size(c)
  n7      = length(v)

  if n7 != n2*n4*n6
      error("a, b, c, and v do not have the correct sizes")
  end

  v_tilda = reshape(v,n4*n6,n2)*a'
  p = vec(kron_prod_times_matrix(b,c,v_tilda))

  return p

end

function vector_times_kron_prod(v::AbstractArray{T,1},a::AbstractArray{T,2},b::AbstractArray{T,2},c::AbstractArray{T,2}) where {T <: Real}

  product = Matrix(kron_prod_times_vector(a',b',c',v')')

  return product

end

function kron_prod_times_matrix(a::AbstractArray{T,2},b::AbstractArray{T,2},c::AbstractArray{T,2},v::AbstractArray{T,2}) where {T <: Real}

  # This function efficiently computes kron(a,kron(b,c))*v where a, b, c and
  # v are comformable matrices.

  (n1,n2) = size(a)
  (n3,n4) = size(b)
  (n5,n6) = size(c)
  (n7,n8) = size(v)

  if n7 != n2*n4*n6
      error("a, b, c, and v do not have the correct sizes")
  end

  p = zeros(n1*n3*n5,n8)
  for i = 1:n8
      @views v_tilda = reshape(v[:,i],n4*n6,n2)*a'
      p[:,i] .= vec(kron_prod_times_matrix(b,c,v_tilda))
  end

  return p

end

function matrix_times_kron_prod(v::AbstractArray{T,2},a::AbstractArray{T,2},b::AbstractArray{T,2},c::AbstractArray{T,2}) where {T <: Real}

  product = Matrix(kron_prod_times_matrix(a',b',c',v')')

  return product

end

function dsylvester(a::AbstractArray{T,2}, b::AbstractArray{T,2}, c::Union{AbstractArray{T,1},AbstractArray{T,2}}) where {T <: Real}

  #= Uses the Hessenberg-Schur method to find the bounded solution of the
     discrete Sylvester equation:

           X + A*X*B = C

     Based on Golub, Nash, and Van Loan (1979).

  =#

  n = size(a, 1)
  m = size(b, 1)
  x = zeros(size(c))

  (s,u) = schur(Matrix(b'))
  (v,t) = hessenberg(a)

  c = v'*c*u

  j = m
  while j > 0
      j1 = j
      if j == 1
          block_size = 1
      elseif isequal(s[j,j-1],0.0) == false
          block_size = 2
          j -= 1
      else;
          block_size = 1
      end
      @views ajj = kron(s[j:j1,j:j1],t) + I
      @views rhs = vec(c[:,j:j1])
      if j1 < m
          @views rhs2 = t*(x[:,(j+1):m]*s[j:j1,(j+1):m]')
          rhs -= vec(rhs2)
      end
      w = ajj\rhs
      @views x[:,j] = w[1:n]
      if block_size == 2
          @views x[:,j1] = w[(n+1):2*n]
      end
      j -= 1
  end

  x = v*x*u'

  return x

end

function trm(x::AbstractArray{T,2}) where {T <: Real}

  # Computes the matrix trace as defined by Binning (2013).  Used for
  # computing the second-order terms, hss, gss.

  (n1,n2) = size(x)
  k = Int(round(sqrt(n2)))

  y = zeros(n1)
  for i = 1:k
      @views y += x[:,i+(i-1)*k]
  end

  return y

end

function trm2(x::AbstractArray{T,2}) where {T <: Real}

  # Computes the matrix trace as defined by Binning (2013).  Used for
  # computing the third-order terms, hssx, gssx.

  (n1,n2) = size(x)
  k = Int(round(n2^(1//3)))

  y = zeros(n1,k)
  for j = 1:k
      for i = 1:k
          @views y[:,j] += x[:,(j-1)+i+(i-1)*k^2]
      end
  end

  return y

end

function create_omega3(n::S) where {S <: Integer}

  # Creates the combination matrix for a third-order perturbation as defined in Levintal (2017).

  # This function is a simplified version of the create_OMEGA function originally written in
  # Matlab by Oren Levintal for his paper "Fifth Order Perturbation Solution to DSGE Models"
  # published in the Journal of Economic Dynamics and Control, 2017.  Permission to translate
  # this function into Julia and release it within the SolveDSGE module was granted by Oren
  # Levintal on February 5, 2020.

  ind = [1:n^3;]
  M = reshape(ind,1,n,n,n)
  Ix = eye(S,n^3)
  omega3 = (reshape(Ix[:,PermutedDimsArray(M,[1,4,2,3])],n^3,n^3)
         + reshape(Ix[:,PermutedDimsArray(M,[1,2,4,3])],n^3,n^3)
         + reshape(Ix[:,PermutedDimsArray(M,[1,2,3,4])],n^3,n^3))

  return omega3

end

function kron_prod_times_vector(A::Union{Array{Array{T,2},1},Array{Array{Complex{T},2},1}},x::Union{Array{T,1},Array{Complex{T},1}},n::Array{S,1},p::S) where {T <: Real, S <: Integer}

  # Computes y = (A[p] * A[p-1] * ... * A[1] )*x

  N = prod(n[1:p])
  z = copy(x)
  for i = 1:p
      z = (A[i]*reshape(z,n[i],Int(N/n[i])))'
  end
  y = reshape(Matrix(z),N)

  return y

end

function kron_prod_times_matrix(A::Union{Array{Array{T,2},1},Array{Array{Complex{T},2},1}},x::Union{Array{T,2},Array{Complex{T},2}},n::Array{S,1},p::S) where {T <: Real, S <: Integer}

  # Computes y = (A[p] * A[p-1] * ... * A[1] )*x

  N = prod(n[1:p])
  y = Array{Complex{T}}(undef,N,size(x,2))
  for j = 1:size(x,2)
      @views z = x[:,j]
      for i = 1:p
          z = (A[i]*reshape(z,n[i],Int(N/n[i])))'
      end
      y[:,j] .= reshape(Matrix(z),N)
  end

  return y

end

function KPShiftSolve(TT::Union{Array{Array{T,2},1},Array{Array{Complex{T},2},1}},n::Array{S,1},c::Union{Array{T,1},Array{Complex{T},1}},lambda::T,alpha::Union{T,Complex{T}}) where {T <: Real, S <: Integer}

  p = length(n)
  N = prod(n)

  c = copy(c)

  TT[p] = alpha*TT[p]
  if p == 1
      y = (TT[1] + lambda*Matrix{T}(I,n[1],n[1]))\c
  else
      y = Array{Complex{T}}(undef,N)
      mp = Int(N/n[p])
      for i = n[p]:-1:1
          idx = ((i-1)*mp+1):(i*mp)
          y[idx] = KPShiftSolve(TT[1:(p-1)],n[1:(p-1)],c[idx],lambda,TT[p][i,i])
          z = kron_prod_times_vector(TT[1:p],y[idx],n,p-1)
          for j = 1:(i-1)
              jdx = ((j-1)*mp+1):(j*mp)
              c[jdx] = c[jdx] - TT[p][j,i]*z
          end
      end
  end

  return y

end

function martin_van_loan(a::Array{T,2},b::Array{T,2},c::Array{T,2},d::Array{T,2},k::S) where {T <: Real, S <: Integer}

  #= Uses a recursive Schur method to find the bounded solution of the Sylvester
     equation:

     AX + B*X*(Kron^(k)C) = D

     Based on Martin and Van Loan (2006).  This is a simplified implementation
     of their algorithm, but it captures most of the gains over Golub, Nash,
     and van Loan (1979).

  =#

  a = copy(a)
  b = a\copy(b)
  c = copy(c)
  d = a\copy(d)

  (v,s) = hessenberg(b)       # v*s*v' = b
  (t,q) = schur(complex(c'))  # q*t*q' = c'

  v = Matrix(v)

  p = k + 2
  TT = Array{typeof(t)}(undef,p)
  TT[1] = s
  for i = 2:p
      TT[i] = conj(t)
  end

  Q     = fill(q,k+1)
  inv_Q = fill(Matrix(q'),k+1)

  n = fill(size(c,1),p)
  n[1] = size(b,1)
  N = prod(n)

  lambda = 1.0

  e = vec(v'*kron_prod_times_matrix(inv_Q,Matrix(d'),n[2:end],k+1)')

  y = reshape(KPShiftSolve(TT,n,e,lambda,1.0),size(d))
  x = real(Matrix(v*kron_prod_times_matrix(Q,Matrix(y'),n[2:end],k+1)'))

  return x

end

function dlyap(a::Array{T,2}, b::Array{T,2}) where {T <: Real}

  n = size(a,1)
  x = zeros(n,n)
  j = n

  (s,u) = schur(a)
  b = u'b*u

  while j > 0
      j1 = j
      if j == 1
          block = 1
      elseif !isequal(s[j,j-1],0.0)
          block = 2
          j -= 1
      else
          block = 1
      end
      @views lhs = kron(s[j:j1,j:j1],s) - I # I = eye(block*n)
      @views rhs = vec(b[:,j:j1])
      if j1 < n
          @views rhs2 = s*(x[:,(j1+1):n]*s[j:j1,(j1+1):n]')
          rhs += vec(rhs2)
      end
      w = -lhs\rhs
      @views x[:,j] = w[1:n]
      if block == 2
          @views x[:,j1] = w[(n+1):block*n]
      end
      j -= 1
  end

  x = u*x*u'

  return x

end

function ind2sub(i::S,dims::Tuple{S,Vararg{S}}) where {S <: Integer}

  if i < 1 || i > prod(dims)
      error("index is out of bounds.")
  end

  subs = CartesianIndices(dims)[i]

  return subs

end

function compute_variances(soln::FirstOrderSolutionStoch)

  hx = soln.hx
  k = soln.k
  gx = soln.gx
  sigma = soln.sigma

  var_states = dlyap(hx,k*sigma*k')
  var_jumps  = gx*var_states*gx'

return var_states,var_jumps

end

function _compute_chebyshev_integrals(eps_nodes::Array{T,1},eps_weights::Array{T,1},nodes::Array{T,1},order::S,rho::T,sigma::T) where {T <: AbstractFloat, S <: Integer}

  # Case where shocks are AR(1) and innovations are independent

  if rho == 0.0
    rho = eps()
  end

  selected_nodes = copy(nodes)
  if length(nodes) > 2
    selected_number = ceil(Int,length(nodes)^(1/2))
    if isodd(length(nodes)) && iseven(selected_number)
      selected_number -= 1
    elseif iseven(length(nodes)) && isodd(selected_number)
      selected_number -= 1
    end
    start_index = Int((length(nodes)-selected_number)/2)
    selected_nodes = nodes[start_index+1:start_index+selected_number]
  end

  terms_num  = Array{T}(undef,length(eps_nodes))
  integrals  = Array{T,2}(undef,order+1,length(selected_nodes))
  integrals2 = Array{T}(undef,order+1)
  for i = 1:(order+1)
    integrals2[i] = sum(exp.(sqrt(2)*sigma*(i-1)*eps_nodes).*eps_weights)*pi^(-1/2)
    for j = 1:length(selected_nodes)
      terms_num     .= rho*nodes[j] .+ sqrt(2)*sigma*eps_nodes
      terms_den      = rho*nodes[j]
      terms_num     .= chebyshev_polynomial(i,terms_num)[:,i]
      terms_den      = chebyshev_polynomial(i,terms_den)[i]
      integrals[i,j] = sum((terms_num/terms_den).*eps_weights)*pi^(-1/2)
     end
  end

  nodetoosmall = abs.(selected_nodes) .< sqrt(eps())
  if sum(nodetoosmall) > 0
    if length(selected_nodes) == 1
      integrals[:,1] .= integrals2
    else
      for i = 1:length(selected_nodes)
        if nodetoosmall[i] == 1
          if i == 1
            integrals[:,i] .= integrals[:,i+1]
        elseif i == length(selected_nodes)
            integrals[:,i] .= integrals[:,i-1]
          else
            integrals[:,i] .= (integrals[:,i-1]+integrals[:,i+1])/2
          end
        end
      end
    end
  end

  return reshape(sum(integrals,dims=2)/length(selected_nodes),order+1)

end

function _compute_chebyshev_integrals(eps_nodes::Array{T,1},eps_weights::Array{T,1},nodes::Array{Array{T,1},1},order::Union{S,Array{S,1}},Ρ::Array{T,2},k::Array{T,2}) where {T <: AbstractFloat, S <: Integer}

N = size(k,2)
for i = 1:N
  if Ρ[i,i] == 0.0
    Ρ[i,i] = eps()
  end
end

selected_nodes = similar(nodes)
for i = 1:length(nodes)
  if length(nodes[i]) > 2
  selected_number = ceil(Int,sqrt(length(nodes[i])))
    if isodd(length(nodes[i])) && iseven(selected_number)
      selected_number -= 1
    elseif iseven(length(nodes[i])) && isodd(selected_number)
      selected_number -= 1
    end
    start_index = Int((length(nodes[i])-selected_number)/2)
    selected_nodes[i] = nodes[i][start_index+1:start_index+selected_number]
  end
end

if typeof(order) == S
  ord = fill(order,N)
else
  ord = copy(order)
end

terms_num  = Array{T}(undef,length(eps_nodes))
integrals  = Array{T}(undef,((ord.+1)...,length.(selected_nodes)...))
integrals2 = Array{T}(undef,(ord.+1)...)

order_prod = prod(ord.+1)
nodes_prod = prod(length.(selected_nodes))
eps_prod   = length(eps_nodes)^N

for i = 1:order_prod
  ii = ind2sub(i,Tuple(ord.+1))
  int2 = 0.0
  for j = 1:eps_prod
    jj = ind2sub(j,Tuple(fill(length(eps_nodes),N)))
    eps_w    = eps_weights[jj[1]]
    eps_node = eps_nodes[jj[1]]
    for k = 2:N
      eps_w    *= eps_weights[jj[k]]
      eps_node = [eps_node;eps_nodes[jj[k]]]
    end
    int2 += exp(2^(N/2)*sqrt(det(k*k'))*collect(ii.-1)'*(k*k')*eps_node)*eps_w*pi^(-N/2)
  end
  integrals2[ii...] = int2 # exp(collect(ii.-1)'*(k*k')*collect(ii.-1)/2) # This is the analytic expression
end

for i = 1:order_prod
  ii = ind2sub(i,Tuple(ord.+1))
  for m = 1:nodes_prod
    mm = ind2sub(m,Tuple(length.(selected_nodes)))
    int2 = 0.0
    node = selected_nodes[1][mm[1]]
    for k = 2:N
      node = [node;selected_nodes[k][mm[k]]]
    end
    for j = 1:eps_prod
      jj = ind2sub(j,Tuple(fill(length(eps_nodes),N)))
      eps_w    = eps_weights[jj[1]]
      eps_node = eps_nodes[jj[1]]
      for k = 2:N
        eps_w    *= eps_weights[jj[k]]
        eps_node = [eps_node;eps_nodes[jj[k]]]
      end
      terms_num = Ρ*node + sqrt(2)*k*eps_node
      terms_den = Ρ*node
      num       = chebyshev_polynomial(ii[1],terms_num[1])[ii[1]]
      den       = chebyshev_polynomial(ii[1],terms_den[1])[ii[1]]
      for k = 2:N
        num = [num;chebyshev_polynomial(ii[k],terms_num[k])[ii[k]]]
        den = [den;chebyshev_polynomial(ii[k],terms_den[k])[ii[k]]]
      end
      int2 += (prod(num)/prod(den))*eps_w*pi^(-N/2)
    end
    integrals[ii...,mm...] = int2
  end
end

for i in eachindex(integrals)
  ii = ind2sub(i,Tuple(size(integrals)))
  if abs(integrals[i]) > 2.0 || isnan(integrals[i])
    if length(integrals) == 1
      integrals[i] = integrals2[i]
    else
      integrals[i] = integrals2[CartesianIndex(Tuple(ii)[1:N])]
    end
  end
end

for i = N:-1:1
  integrals = sum(integrals,dims = (N+i))/length(selected_nodes[i])
end

return reshape(integrals,Tuple(ord.+1))

end

function compute_chebyshev_integrals(eps_nodes::Array{T,1},eps_weights::Array{T,1},nodes::Array{Array{T,1},1},order::Union{S,Array{S,1}},Ρ::Array{T,2},k::Array{T,2}) where {T <: AbstractFloat, S <: Integer}

if !isdiag(Ρ.>sqrt(eps()))
  error("The autoregression matrix for the shocks must be diagonal")
end

ns = size(k,2)

if typeof(order) == S
  ord = fill(order,ns)
else
  ord = order[1:ns]
end

if !isdiag(abs.(k).>sqrt(eps())) # Correlated innovations
  integrals = _compute_chebyshev_integrals(eps_nodes,eps_weights,nodes[1:ns],ord,Ρ,k)
  return integrals
else
  integrals = Array{Array{T,1},1}(undef,ns)
  for i = 1:ns
    integrals[i] = _compute_chebyshev_integrals(eps_nodes,eps_weights,nodes[i],ord[i],Ρ[i,i],k[i,i])
  end
  return integrals
end

end

function scale_chebyshev_weights!(weights::Array{Array{T,N},1},scaled_weights::Array{Array{T,N},1},integrals::Array{Array{T,1},1},j_approx::Union{S,Array{S,1}},ns::S) where {T <: AbstractFloat, N, S <: Integer}

  for i = 1:length(j_approx)
      for j = 1:ns
          index = [1:ndims(weights[i]);]
          index[1],index[j] = index[j],index[1]
          scaled_weights[i] .= permutedims(integrals[j].*permutedims(weights[i],index),index)
      end
  end

end

function scale_chebyshev_weights!(weights::Array{Array{T,N},1},scaled_weights::Array{Array{T,N},1},integrals::Array{T,N2},j_approx::Union{S,Array{S,1}},ns::S) where {T <: AbstractFloat, N, N2, S <: Integer}

  for i = 1:length(j_approx)
      for j = 1:ns
          scaled_weights[i] .= integrals.*weights[i]
      end
  end

end

function _compute_smolyak_integrals(eps_nodes::Array{T,1},eps_weights::Array{T,1},nodes::Array{T,1},order::S,rho::T,sigma::T) where {T <: AbstractFloat, S <: Integer}

  # Case where shocks are AR(1) and innovations are independent

  if rho == 0.0
    rho = eps()
  end

  terms_num  = Array{T}(undef,length(eps_nodes))
  integrals  = Array{T,2}(undef,order+1,length(nodes))
  integrals2 = Array{T}(undef,order+1)
  for i = 1:(order+1)
    integrals2[i] = sum(exp.(sqrt(2)*sigma*(i-1)*eps_nodes).*eps_weights)*pi^(-1/2)
    for j = 1:length(nodes)
      terms_num     .= rho*nodes[j] .+ sqrt(2)*sigma*eps_nodes
      terms_den      = rho*nodes[j]
      terms_num     .= chebyshev_polynomial(i,terms_num)[:,i]
      terms_den      = chebyshev_polynomial(i,terms_den)[i]
      integrals[i,j] = sum((terms_num/terms_den).*eps_weights)*pi^(-1/2)
     end
  end

  nodetoosmall = abs.(nodes) .< sqrt(eps())
  if sum(nodetoosmall) > 0
    if length(nodes) == 1
      integrals[:,1] .= integrals2
    else
      for i = 1:length(nodes)
        if nodetoosmall[i] == 1
          if i == 1
            integrals[:,i] .= integrals[:,i+1]
        elseif i == length(nodes)
            integrals[:,i] .= integrals[:,i-1]
          else
            integrals[:,i] .= (integrals[:,i-1]+integrals[:,i+1])/2
          end
        end
      end
    end
  end

  return reshape(sum(integrals,dims=2)/length(nodes),order+1)

end

function compute_smolyak_integrals(eps_nodes,eps_weights,nx,order,grid,RHO,k)

integrals = ones(nx,order+1)
for j = 1:size(k,2)
    nodes = unique(grid[:,j])
    integrals[j,:] .= _compute_smolyak_integrals(eps_nodes,eps_weights,nodes,order,RHO[j,j],k[j,j])
end

#  for j = 1:size(sigma,2)
#      for i = 1:(order+1)
#          integrals[j,i] = sum(exp.(sqrt(2)*k[j,j]*(i-1)*eps_nodes).*eps_weights)*pi^(-1/2)
#      end
#  end

return integrals

end

function weight_scale_factors(eps_nodes,eps_weights,multi_index,nx,grid,RHO,sigma)

unique_multi_index = sort(unique(multi_index))
unique_orders = SmolyakApprox.m_i(unique_multi_index).-1

# Here we construct the base integrals

base_integrals = Array{Array{Float64,2}}(undef,length(unique_orders))
for i = 1:length(unique_orders)
  base_integrals[i] = compute_smolyak_integrals(eps_nodes,eps_weights,nx,unique_orders[i],grid,RHO,sigma)
end

# Compute the unique polynomial terms from the base polynomials

unique_base_integrals = Array{Array{Float64,2}}(undef,length(unique_orders))
for i = length(unique_orders):-1:2
  unique_base_integrals[i] = base_integrals[i][:,size(base_integrals[i-1],2)+1:end]
end
unique_base_integrals[1] = base_integrals[1]

# Construct the first row of the interplation matrix

new_integrals = unique_base_integrals[multi_index[1,1]][1,:]
for i = 2:size(multi_index,2)
  new_integrals = kron(new_integrals,unique_base_integrals[multi_index[1,i]][i,:])
end

weight_scale_factor = copy(new_integrals)

# Iterate over nodes, doing the above three steps at each iteration

for j = 2:size(multi_index,1)

  new_integrals = unique_base_integrals[multi_index[j,1]][1,:]
  for i = 2:size(multi_index,2)
    new_integrals = kron(new_integrals,unique_base_integrals[multi_index[j,i]][i,:])
  end
  weight_scale_factor = [weight_scale_factor; new_integrals]

end

return weight_scale_factor

end

function scale_smolyak_weights(weights,weight_scale_factor)

  scaled_weights = weights.*weight_scale_factor

end

function compute_piecewise_linear_integrals(eps_nodes,eps_weights,sigma)

integral = 1.0
#integral = sum(exp.(sqrt(2)*sigma*eps_nodes).*eps_weights)*pi^(-1/2)

return integral

end

function mylen(len::T,a::Array{T,N},b::Array{T,N}) where {T <: AbstractFloat, N}

  for i in eachindex(a)
      len = max(len,abs(a[i]-b[i]))
  end

  return len

end


text:
function plot_summary(summary::RobustnessSummary, X; nbins = 8, image_shape = (16, 8),
    reshape_im = x -> reshape(x, 16, 8)', threshold = 0.9, masked_opacity = 0.2)
    """ Plots parts of the robustness summary
        Args
            - summary: robustness summary
            - X: training examples (each column is an example)
            --------------
            - nbins: number of bins for histograms
            - reshape_im: function to reshape vectorized images to images
            - threshold: highlight pixels above this percentile
            - masked_opacity: opacity for non-highlighted pixels
    """
    default(titlefont = 30, guidefont = 24, tickfont = 20)

    xmax = 1.05 * max(max(maximum(summary.errors_nominal), maximum(summary.errors_adversarial_true)),
    maximum(summary.mae_adversarial_pred))

    # Histograms
    bins = 0:(xmax / nbins):xmax
    h1 = fit(Histogram, summary.errors_nominal, bins)
    h2 = fit(Histogram, summary.errors_adversarial_true, bins)
    h3 = fit(Histogram, summary.errors_adversarial_pred, bins)

    ymax = maximum([maximum(h1.weights), maximum(h2.weights), maximum(h3.weights)])

    p1 = plot(h1, title = "Nominal Absolute Errors", xlims = (0.0, xmax),
    ylims = (0.0, ymax), legend = false, color = :teal, xlabel = "Absolute Error")
    p2 = plot(h2, title = "Adversarial Absolute Errors", xlims = (0.0, xmax),
    ylims = (0.0, ymax), legend = false, color = :teal, xlabel = "Absolute Error")
    p3 = plot(h3, title = "Adversarial Perturbations", xlims = (0.0, xmax),
    ylims = (0.0, ymax), legend = false, color = :teal, xlabel = "Absolute Error")

    # Explanations
    p4 = plot_explanations(summary, X, reshape_im = reshape_im,
    threshold = threshold, masked_opacity = masked_opacity)

    l = @layout [a b c ; d]
    return plot(p1, p2, p3, p4, layout = l, size = (2200, 1600), margin = 10mm)
end

function plot_nominal_errors(summary::RobustnessSummary; xmax = nothing)
    if xmax === nothing
        xmax = 1.05 * maximum(summary.errors_nominal)
    end

    p = Plots.histogram(summary.errors_nominal, title = "Nominal Absolute Errors",
    xlims = (0.0, xmax), legend = false, xlabel = "Absolute Error", size = (1000, 1000))

    return p
end

function plot_adversarial_errors(summary::RobustnessSummary; xmax = nothing)
    if xmax === nothing
        xmax = 1.05 * maximum(summary.errors_adversarial_true)
    end

    p = Plots.histogram(summary.errors_adversarial_true,
    title = "Adversarial Absolute Errors", xlims = (0.0, xmax), legend = false,
    xlabel = "Absolute Error", size = (1000, 1000))

    return p
end

function plot_adversarial_perturbations(summary::RobustnessSummary; xmax = nothing)
    if xmax === nothing
        xmax = 1.05 * maximum(summary.errors_adversarial_pred)
    end

    p = Plots.histogram(summary.errors_adversarial_pred,
    title = "Adversarial Perturbations", xlims = (0.0, xmax), legend = false,
    xlabel = "Output Perturbation", size = (1000, 1000))

    return p
end

function plot_explanations(summary::RobustnessSummary, X; N = 5,
    reshape_im = x -> reshape(x, 16, 8), threshold = 0.9, masked_opacity = 0.2)

    X_samples = X[:, summary.samples]
    sample_im = reshape_im(X_samples[:, 1])
    m, n = size(sample_im)

    total_im = ones(m * 3 + 2, n * N + (N - 1))

    for i = 1:N
        im = reshape_im(X_samples[:, i])
        total_im[1:m, n * (i-1) + 1 + (i - 1):n * i + (i - 1)] = im

        attrs_sal = summary.saliency_maps[:, i]
        masked_attrs_sal = get_mask(abs.(attrs_sal), X_samples[:, i], threshold = threshold,
                            masked_opacity = masked_opacity)
        masked_im_sal = reshape_im(masked_attrs_sal)
        total_im[m + 2:2 * m + 1, n * (i-1) + 1 + (i - 1):n * i + (i - 1)] = masked_im_sal

        attrs_ig = summary.ig_maps[:, i]
        masked_attrs_ig = get_mask(abs.(attrs_ig), X_samples[:, i], threshold = threshold,
                            masked_opacity = masked_opacity)
        masked_im_ig = reshape_im(masked_attrs_ig)
        total_im[2 * m + 3:3 * m + 2, n * (i-1) + 1 + (i - 1):n * i + (i - 1)] = masked_im_ig
    end

    p = plot(Gray.(total_im), axis = [],
        title = "Explanations (Original, Saliency Map, Integrated Gradients)")

    return p
end


text:
image_1 = img
image_2 = imgw
filter_half_size = fhs
window_size = window_size

## function
image_size = size(image_1)
pixel_count = length(image_1)
point_count = size(points, 2)
filter_size = 2 * filter_half_size + 1

# Prepare filters:
# Calculate separable filters from basis:
sigma = (filter_half_size + 2) / 4
centered_inds = centered(-filter_half_size:filter_half_size)
gaus = gaussian(sigma, filter_size)
gaus_der = gaus .* centered_inds .* (-1)/sigma^2
gaus_der_flip = reverse(gaus_der, dims=1)

# prepare 2D filter basis for later use
basis = zeros(filter_size, filter_size, filter_num)

# dif = (points of forward image - points of backward image) for each filter
dif = zeros(image_size..., filter_num)

if filter_num == 3
    # temporary place to store filtered images
    tmp_filtered_1 = zeros(image_size)
    tmp_filtered_2 = zeros(image_size)
    tmp_filtered_1 = fill(NaN, image_size)
    tmp_filtered_2 = fill(NaN, image_size)

    # basis 1 - (gaus, gaus) both forward and backward
    kernf_1 = kernelfactors((gaus, gaus))
    basis[:, :, 1] = broadcast(*, kernf_1...)
    filt_onebyone!(tmp_filtered_1, image_1, kernf_1, filter_half_size, points)
    filt_onebyone!(tmp_filtered_2, image_2, kernf_1, filter_half_size, points)
    dif[:, :, 1] = tmp_filtered_2 - tmp_filtered_1

    # basis 2 - (gaus, gaus_der_flip) as forward, (gaus, gaus_der) as backward
    kernf_2f = kernelfactors((gaus, gaus_der_flip))
    basis[:, :, 2] = broadcast(*, kernf_2f...)
    filt_onebyone!(tmp_filtered_1, image_1, kernf_2f, filter_half_size, points)
    kernf_2b = kernelfactors((gaus, gaus_der))
    filt_onebyone!(tmp_filtered_2, image_2, kernf_2b, filter_half_size, points)
    dif[:, :, 2] = tmp_filtered_2 - tmp_filtered_1

    # basis 3 - (gaus_der_flip, gaus) as forward, (gaus_der, gaus) as backward
    kernf_3f = kernelfactors((gaus_der_flip, gaus))
    basis[:, :, 3] = broadcast(*, kernf_3f...)
    filt_onebyone!(tmp_filtered_1, image_1, kernf_3f, filter_half_size, points)
    kernf_3b = kernelfactors((gaus_der, gaus))
    filt_onebyone!(tmp_filtered_2, image_2, kernf_3b, filter_half_size, points)
    dif[:, :, 3] = tmp_filtered_2 - tmp_filtered_1
end

dif = reshape(dif, (:, filter_num))

# prepare matrices for linear system of equations
A = zeros(filter_num-1, filter_num-1, pixel_count)
b = zeros(filter_num-1, pixel_count)
A = fill(NaN, filter_num-1, filter_num-1, pixel_count)
b = fill(NaN, filter_num-1, pixel_count)

for k in 1:filter_num-1
    for l in k:filter_num-1
        @views window_sum_around_points!(A[k, l, :], dif[:, k+1] .* dif[:, l+1], image_size, window_size, points)
        A[l, k, :] = A[k, l, :]
    end
    @views window_sum_around_points!(b[k, :], dif[:, k+1] .* dif[:, 1] .* (-1), image_size, window_size, points)
end

println(length(filter(!isnan, real.(A[1, 1, :]))))
imgshow(reshape(A[1, 1, :], image_size...))
println(filter(!isnan, real.(A[1, 1, :])))
println(filter(!isnan, real.(A[2, 1, :])))
println(filter(!isnan, real.(A[2, 2, :])))
println(filter(!isnan, real.(A[1, 2, :])))

# Perform Gauss elimination on all pixels in parallel:
# coeffs will be of shape: pixel_count, filter_num-1
@views coeffs = multi_mat_div_at_points(A, b, points, image_size)
# adding ones so that all base filters have their coefficients even the first one
all_coeffs = [ones(pixel_count) coeffs]

# make a border mask. Is 0 if its in a border of filter_half_size size.
window_half_size = Int64.((window_size .- 1) ./ 2)
border_mask = parent(padarray(ones((image_size .- (2 .* window_half_size))...),
                Fill(NaN, window_half_size, window_half_size)))
@views all_coeffs = all_coeffs .* reshape(border_mask, (:, 1))

k = (-filter_half_size:filter_half_size)

# Get the displacement vector field from the filters
u1_top = zeros(image_size);
u1_bot = zeros(image_size);
u2_top = zeros(image_size);
u2_bot = zeros(image_size);

for n in 1:filter_num
    @views u1_top[:] = u1_top[:] .- sum(transpose(basis[:, :, n]) * k) .* all_coeffs[:, n];
    @views u1_bot[:] = u1_bot[:] .+ sum(basis[:, :, n]) .* all_coeffs[:, n];

    @views u2_top[:] = u2_top[:] .- sum(basis[:, :, n] * k) .* all_coeffs[:, n];
    @views u2_bot[:] = u2_bot[:] .+ sum(basis[:, :, n]) .* all_coeffs[:, n];
end

u_est = 2 .* ((im .* u1_top ./ u1_bot) .+ (u2_top ./ u2_bot));

# dont use estimations whose displacement is larger than the filter_half_size
displacement = real(u_est).^2 + imag(u_est).^2
displacement_mask = displacement .> filter_half_size^2
u_est[displacement_mask] .= NaN .+ NaN .* 1im;

any(.!isnan.(real(u_est)))

showflow(u_est)


text:
using ApproxFun, Plots

x=Fun()
u0=0.0x
N=u->[u(-1.),u(1.),0.001u''+6*(1-x^2)*u'+u^2-1.]
u=newton(N,u0)

x=Fun()
u0=0.0x
N=u->[u(-1.),u(1.),0.1u''+sin(u)-1.]
u=newton(N,u0)

##
# Heaven and hell for PI
##

M = 10.0
x = Fun(0..M)
N = u -> [u(0), u(M)-sqrt(M), u'' - u^2 + x]
u0 = 0.0x
@time heaven = newton(N, u0)

u0 = heaven-10x * exp(-x)
@time hell = newton(N, u0)
plot(heaven; label="heaven", linewidth=1.5, legend=:bottomright, title="Painleve I: Heaven and Hell")
plot!(hell; label="hell", linewidth=1.5, linestyle=:dash)
savefig("heavenhell.pdf")

scatter(abs.(heaven.coefficients); yscale=:log10,  label="heaven", markersize=1.5, title="Painleve I: Heaven and Hell Coefficients")
scatter!(abs.(hell.coefficients); yscale=:log10,  label="hell", markersize=1.5, marker=:square)
savefig("heavenhellcoefficients.pdf")

max_stars_repo_path: ODEs/Newton iteration for Nonlinear BVPs.jl

max_stars_repo_name: Gyoshi/ApproxFunExamples

max_issues_repo_name: Gyoshi/ApproxFunExamples

max_forks_repo_name: Gyoshi/ApproxFunExamples


text:
using CKANQuery
using Test

@testset "CKANQuery.jl" begin
    # Write your tests here.
end


text:
using Test

include("frames.jl")


text:
@testset "Unions" begin
    @testset "Low-level" begin
        test(StructDefinition{false}, union_by_name, :VkClearColorValue, :(
            struct _ClearColorValue <: VulkanStruct{false}
                vks::VkClearColorValue
            end
        ))

        test(StructDefinition{false}, union_by_name, :VkClearValue, :(
            struct _ClearValue <: VulkanStruct{false}
                vks::VkClearValue
            end
        ))
    end

    @testset "High-level" begin
        test(StructDefinition{true}, union_by_name, :VkClearColorValue, :(
            struct ClearColorValue <: HighLevelStruct
                data::VkClearColorValue
            end
        ))

        test(StructDefinition{true}, union_by_name, :VkClearValue, :(
            struct ClearValue <: HighLevelStruct
                data::VkClearValue
            end
        ))
    end
end


text:
function H5TWRI()

end


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

using Markdown
using InteractiveUtils

# ╔═╡ 8dedd8b2-0f0a-4523-8503-3c8405b3ef64
using PlutoUI

# ╔═╡ cad5438e-df75-11eb-1882-e384f5a8554c
md"""
# Chapter 6
#### Nicetak
 $(import Dates; Dates.format(Dates.today(), Dates.DateFormat("U d, Y")))
"""

# ╔═╡ c1e55426-349f-4419-aacc-756b9db75ce5
md"""
## 6.1
どの２要素も互いに相異なる$N$要素からなる整数列$a_0, a_1, \dots, a_{N-1}$が与えられます.
$i = 0, 1, \dots, N-1$に対して, $a_i$が全体の中で何番目に小さい値であるかを$O(N \log N)$で求めるアルゴリズムを設計してください.
"""

# ╔═╡ bb6c9fc1-1087-41c2-9685-38734bda7abe
md"""
**答**
Juliaの標準ソートアルゴリズム`sort`が$O(N \log N)$で動作し, 各$a_i$がソート列のどの位置に挿入されるかは`searchsortedfirst`によって$O(\log N)$で調べられ(二分探索より), それを$N$回繰り返せばよい.
"""

# ╔═╡ bfc7a793-e646-42bf-9f93-246ffe0f10f1
function f₁(as)
	as_sorted = sort(as)
	return [searchsortedfirst(as_sorted, a) for a in as]
end

# ╔═╡ cdda7a2e-a79b-4585-b39a-c7efc39af10e
md"""
なお, Juliaには標準でこの座標圧縮を行う関数`sortperm`が実装されている.
"""

# ╔═╡ db96c894-8e17-46c4-867b-ce8628c80939
with_terminal() do
	as = [12, 43, 7, 15, 9]
	println("Given as = $as,\n")
	@show f₁(as)
	@show sortperm(as)
end

# ╔═╡ 0b3ae17a-cbed-4b12-b9fb-1ff125ae7774
md"""
 $N = 100$に対して, $M = 10000$回の試行を行う計算時間を比べると, `sortperm`がやや早い.
"""

# ╔═╡ 56f632a0-e7c8-46da-a9b2-34bb96f053a4
begin
	function time₁(A, f)
		M = size(A)[2]

		for i in 1:M
			f(A[:, i])
		end
	end

	N, M = 100, 10000
	A = rand(1:N, N, M)

	t_f₁ = @timed time₁(A, f₁)
	t_sp = @timed time₁(A, sortperm)

	with_terminal() do
		println("Time f₁: $(t_f₁.time)")
		println("Time sortperm: $(t_sp.time)")
	end
end

# ╔═╡ Cell order:
# ╟─cad5438e-df75-11eb-1882-e384f5a8554c
# ╟─8dedd8b2-0f0a-4523-8503-3c8405b3ef64
# ╟─c1e55426-349f-4419-aacc-756b9db75ce5
# ╟─bb6c9fc1-1087-41c2-9685-38734bda7abe
# ╠═bfc7a793-e646-42bf-9f93-246ffe0f10f1
# ╟─cdda7a2e-a79b-4585-b39a-c7efc39af10e
# ╟─db96c894-8e17-46c4-867b-ce8628c80939
# ╟─0b3ae17a-cbed-4b12-b9fb-1ff125ae7774
# ╟─56f632a0-e7c8-46da-a9b2-34bb96f053a4


text:
x = [0.0, 1.0, 2.0]                        # function knots
u = x.^2                                   # function values in knots
s = [2.0]                                  # function first derivative knots
v = [4.0]                                  # function first derivative values
t = [0.0, 1.0]                             # function second derivative knots
w = [2.0 ,2.0]                             # function second derivative values
d1_u = v                                   # exact function first derivative values in knots
d2_u = w                                   # exact function second derivative values in knots

@testset "Bessel potentials space reproducing kernels" begin
    @testset "RK_H1 kernel" begin
        interpolate(x, u, RK_H1(0.1))  # create spline
        cond = get_cond()              # get estimation of the Gram matrix condition number
        @test cond ≈ 100.0

        σ = evaluate(x)                # evaluate spline in knots
        @test σ ≈ u                    # compare with exact function values in knots

        # Check that we get close when evaluating near the knots
        p = x .+ 1e-3*randn(size(x))   # evaluation points near the knots
        f = p.^2                       # exact function values in evaluation points
        σ = evaluate(p)                # evaluate spline in evaluation points
        # compare spline values with exact function values in evaluation point
        @test all(isapprox.(σ, f, atol = 0.05))
    end

    @testset "RK_H2 kernel" begin
        interpolate(x, u, RK_H2(0.1))  # create spline
        cond = get_cond()              # get estimation of the gram matrix condition number
        @test cond ≈ 1.0e5

        σ = evaluate(x)                # evaluate spline in knots
        @test σ ≈ u                    # compare with exact function values in knots

        # Check that we get close when evaluating near the knots
        p = x .+ 1e-3*randn(size(x))   # evaluation points near the knots
        f = p.^2                       # exact function values in evaluation points
        σ = evaluate(p)                # evaluate spline in evaluation points
        # compare spline values with exact function values in evaluation point
        @test all(isapprox.(σ, f, atol = 1e-2))

        ###
        interpolate(x, u, s, v, RK_H2(0.1)) # create spline by function and
                                            # first derivative values in knots
        cond = get_cond()              # get estimation of the gram matrix condition number
        @test cond ≈ 1.0e5

        σ = evaluate(x)                # evaluate spline in knots
        @test σ ≈ u                    # compare with exact function values in knots

        # Check that we get close when evaluating near the knots
        p = x .+ 1e-3*randn(size(x))   # evaluation points near the knots
        f = p.^2                       # exact function values in evaluation points
        σ = evaluate(p)                # evaluate spline in evaluation points
        # compare spline values with exact function values in evaluation point
        @test all(isapprox.(σ, f, atol = 1e-2))

        d1_σ = evaluate(s, 1)             # evaluate spline first derivative in knots
        @test d1_σ ≈ d1_u                 # compare with exact function first derivative values in knots

        # Check that we get close when evaluating near the knots
        ps = s .+ 1e-3*randn(size(s))   # evaluation points near the knots
        d1_f = ps.*2                    # exact function first derivative values in evaluation points
        d1_σ = evaluate(ps, 1)             # evaluate spline first derivative in evaluation points
        # compare spline first derivative values with exact function first derivative values in evaluation point
        @test all(isapprox.(d1_σ, d1_f, atol = 1e-2))
    end

    @testset "RK_H3 kernel" begin
        interpolate(x, u, RK_H3(0.1))  # create spline
        cond = get_cond()              # get estimation of the gram matrix condition number
        @test cond ≈ 1.0e7

        σ = evaluate(x)                # evaluate spline in knots
        @test σ ≈ u                    # compare with exact function values in knots

        # Check that we get close when evaluating near the knots
        p = x .+ 1e-3*randn(size(x))   # evaluation points near the knots
        f = p.^2                       # exact function values in evaluation points
        σ = evaluate(p)                # evaluate spline in evaluation points
        # compare spline values with exact function values in evaluation point
        @test all(isapprox.(σ, f, atol = 1e-2))

        ###
        interpolate(x, u, s, v, RK_H3(0.1))     # create spline
        cond = get_cond()              # get estimation of the gram matrix condition number
        @test cond ≈ 1.0e8

        σ = evaluate(x)                # evaluate spline in knots
        @test σ ≈ u                    # compare with exact function values in knots

        # Check that we get close when evaluating near the knots
        p = x .+ 1e-3*randn(size(x))   # evaluation points near the knots
        f = p.^2                       # exact function values in evaluation points
        σ = evaluate(p)                # evaluate spline in evaluation points
        # compare spline values with exact function values in evaluation point
        @test all(isapprox.(σ, f, atol = 1e-2))

        d1_σ = evaluate(s, 1)             # evaluate spline first derivative in knots
        @test d1_σ ≈ d1_u                 # compare with exact function first derivative values in knots

        # Check that we get close when evaluating near the knots
        ps = s .+ 1e-3*randn(size(s))   # evaluation points near the knots
        d1_f = ps.*2                    # exact function first derivative values in evaluation points
        d1_σ = evaluate(ps, 1)             # evaluate spline first derivative in evaluation points
        # compare spline first derivative values with exact function first derivative values in evaluation point
        @test all(isapprox.(d1_σ, d1_f, atol = 1e-2))

        ###
        interpolate(x, u, t, w, RK_H3(0.1), 2)     # create spline
        cond = get_cond()              # get estimation of the gram matrix condition number
        @test cond ≈ 1.0e8

        σ = evaluate(x)                # evaluate spline in knots
        @test σ ≈ u                    # compare with exact function values in knots

        # Check that we get close when evaluating near the knots
        p = x .+ 1e-3*randn(size(x))   # evaluation points near the knots
        f = p.^2                       # exact function values in evaluation points
        σ = evaluate(p)                # evaluate spline in evaluation points
        # compare spline values with exact function values in evaluation point
        @test all(isapprox.(σ, f, atol = 1e-2))

        d1_σ = evaluate(s, 1)          # evaluate spline first derivative in knots
        @test d1_σ ≈ d1_u atol = 0.05  # compare with exact function first derivative values in knots

        # Check that we get close when evaluating near the knots
        ps = s .+ 1e-3*randn(size(s))   # evaluation points near the knots
        d1_f = ps.*2                    # exact function first derivative values in evaluation points
        d1_σ = evaluate(ps, 1)          # evaluate spline first derivative in evaluation points
        # compare spline first derivative values with exact function first derivative values in evaluation point
        @test all(isapprox.(d1_σ, d1_f, atol = 0.05))

        d2_σ = evaluate(t, 2)           # evaluate spline second derivative in knots
        @test d2_σ ≈ d2_u               # compare with exact function second derivative values in knots

        # Check that we get close when evaluating near the knots
        pt = t .+ 1e-3*randn(size(t))   # evaluation points near the knots
        d2_f = 2.0                      # exact function second derivative values in evaluation points
        d2_σ = evaluate(pt, 2)          # evaluate spline second derivative in evaluation points
        # compare spline second derivative values with exact function second derivative values in evaluation point
        @test all(isapprox.(d2_σ, d2_f, atol = 1e-2))

        ###
        interpolate(x, u, s, v, t, w, RK_H3(0.1)) # create spline
        cond = get_cond()              # get estimation of the gram matrix condition number
        @test cond ≈ 1.0e9

        σ = evaluate(x)                # evaluate spline in knots
        @test σ ≈ u                    # compare with exact function values in knots

        # Check that we get close when evaluating near the knots
        p = x .+ 1e-3*randn(size(x))   # evaluation points near the knots
        f = p.^2                       # exact function values in evaluation points
        σ = evaluate(p)                # evaluate spline in evaluation points
        # compare spline values with exact function values in evaluation point
        @test all(isapprox.(σ, f, atol = 1e-2))

        d1_σ = evaluate(s, 1)          # evaluate spline first derivative in knots
        @test d1_σ ≈ d1_u              # compare with exact function first derivative values in knots

        # Check that we get close when evaluating near the knots
        ps = s .+ 1e-3*randn(size(s))  # evaluation points near the knots
        d1_f = ps.*2                   # exact function first derivative values in evaluation points
        d1_σ = evaluate(ps, 1)         # evaluate spline first derivative in evaluation points
        # compare spline first derivative values with exact function first derivative values in evaluation point
        @test all(isapprox.(d1_σ, d1_f, atol = 1e-2))

        d2_σ = evaluate(t, 2)          # evaluate spline second derivative in knots
        @test d2_σ ≈ d2_u              # compare with exact function second derivative values in knots

        # Check that we get close when evaluating near the knots
        pt = t .+ 1e-3*randn(size(t))   # evaluation points near the knots
        d2_f = 2.0                      # exact function second derivative values in evaluation points
        d2_σ = evaluate(pt, 2)          # evaluate spline second derivative in evaluation points
        # compare spline second derivative values with exact function second derivative values in evaluation point
        @test all(isapprox.(d2_σ, d2_f, atol = 1e-2))
    end

end


text:
using Plots
using BenchmarkTools

function vfsolvex(vnew, kgrid, tolerance, imax, σ)
    β = 0.9

    v = vnew .+ 2*tolerance
    cartesianindex = Array{CartesianIndex{2}, length(v)}
    i = 1

    c = kgrid .- kgrid'
    c[c .< 0] .= 0
    u = (c .^ (1-σ) .- 1) ./ (1-σ)
    u[kgrid .- kgrid' .< 0] .= -Inf

    while maximum(abs.(v - vnew)) > tolerance && i <= imax
        v = vnew;

        (vnew, cartesianindex) = findmax(u .+ β*v', dims = 2);
        i += 1;
    end
    kprimeindex = getindex.(cartesianindex, 2)
    return (v = vnew, kprime = kgrid[kprimeindex], kprimeindex = kprimeindex)
end

kupper = 5
klower = 0.01
n = 1000
kgrid = collect(range(klower, stop = kupper, length = n))
(v, kprime, kprimeindex) = vfsolvex(zeros(n), kgrid, 0.001, 1000, 1.5);

scatter(kgrid, v, label = "v")
scatter(kgrid, kprime, label = "k'")

function policy(x)
    e=getindex(findall(kgrid .== x),1)
    return kprime[e]
end

T = 100
kpath = zeros(T)
k = kgrid[1000]
for i in 1:T
    k = policy(k)
    kpath[i] = k
end
time = collect(range(1, T, length = T))
scatter(time, kpath, label = "kpath")

max_stars_repo_path: Dynamic Programming/DM2.jl

max_issues_repo_path: Dynamic Programming/DM2.jl

max_forks_repo_path: Dynamic Programming/DM2.jl