text:
using ArgParse
using Formatting
using JSON,TOML
using DataStructures

include("runCmd.jl")

s = ArgParseSettings()
@add_arg_table s begin

    "--casename"
        help = "Prefix of the file."
        arg_type = String
        required = true

    "--input-dir"
        help = "The folder where case folders are contained"
        arg_type = String
        required = true

    "--output-dir"
        help = "Output folder"
        arg_type = String
        required = true

    "--year-rng"
        help = "Range of years being diagnosed"
        arg_type = Int64
        nargs = 2
        required = true

    "--domain-file"
        help = "Domain file of ocean model"
        arg_type = String
        required = true

    "--overwrite"
        action = :store_true

    "--remap-file-nn"
        help = "Remap weighting file"
        arg_type = String
        required = true

    "--pop2"
        arg_type = Bool
        default = false
end

parsed = DataStructures.OrderedDict(parse_args(ARGS, s))

JSON.print(parsed, 4)

for y = parsed["year-rng"][1]:parsed["year-rng"][2], m=1:12

    date_str = format("{:04d}-{:02d}", y, m)

    if ! parsed["pop2"]
        old_file  = joinpath(parsed["input-dir"],  "$(parsed["casename"]).EMOM.h0.mon.$(date_str).nc")
        Nx = "Nx"
        Ny = "Ny"
        Nz = "Nz"
    else
        old_file  = joinpath(parsed["input-dir"],  "$(parsed["casename"]).pop.h.$(date_str).nc")
        Nx = "nlon"
        Ny = "nlat"
        Nz = "z_t"
    end

    new_file1 = joinpath(parsed["output-dir"], "$(parsed["casename"]).EMOM_extra1.$(date_str).nc")
    new_file2 = joinpath(parsed["output-dir"], "$(parsed["casename"]).EMOM_extra1_rg.$(date_str).nc")

    if !isfile(new_file1) || !isfile(new_file2) || parsed["overwrite"]

        println("Generating file: $(new_file1)")

#        pleaseRun(`ncap2 -h -O -v
#                    -s "SST=array(0.0,0.0,/\$time,\$$Ny,\$$Nx/);  SST(0,:,:)=TEMP(0, 0, :, :);"
#                    -s '*tmp=HMXL*1.0;'
#                    -s "HMXL=array(0.0,0.0,/\$time,\$$Ny,\$$Nx/); HMXL(0,:,:)=tmp;"
#                    $old_file $new_file1
#        `; igs=false)

        pleaseRun(`ncap2 -h -O -v
                    -s "SST=array(0.0,0.0,/\$time,\$$Ny,\$$Nx/);  SST(0,:,:)=TEMP(0, 0, :, :);"
                    $old_file $new_file1
        `; igs=false)

        if parsed["pop2"]
            pleaseRun(`ncrename -d nlat,Ny -d nlon,Nx -d z_t,Nz $new_file1`; igs=true)
        end
        println("Generating file: $(new_file2)")
        pleaseRun(`ncremap -R '--rgr lat_nm_in=Ny --rgr lon_nm_in=Nx' -m $(parsed["remap-file-nn"]) $new_file1 $new_file2`;igs=false)

    end

    # Special: stratification
    new_file3 = joinpath(parsed["output-dir"], "$(parsed["casename"]).EMOM_extra2.$(date_str).nc")
    new_file4 = joinpath(parsed["output-dir"], "$(parsed["casename"]).EMOM_extra2_rg.$(date_str).nc")
    if !isfile(new_file3) || parsed["overwrite"]

        println("Generating file: $(new_file3)")

        if !parsed["pop2"]
            dz_cmd = "*dz = dz_cT(:, 0, 0);"
        else
            dz_cmd = "*dz=array(0.0,0.0,/\$$Nz/); for (i=0;i<\$z_t.size;i++) { dz(i) = (z_w_bot(i) - z_w_top(i))/100.0; };"
        end
        pleaseRun(`ncap2 -h -O -v
                    -s 'defdim("Nz1", 5); defdim("Nz2", 28);'
                    -s "$dz_cmd"
                    -s '*dz1=dz(0:4).ttl(); *dz2=dz(5:32).ttl();'
                    -s "*TEMPdz[\$time,\$$Nz,\$$Ny,\$$Nx]=TEMP*dz; TEMPdz.set_miss(TEMP@_FillValue);"
                    -s "*_T1=array(0.0,0.0,/\$time,\$Nz1,\$$Ny,\$$Nx/);_T1.set_miss(TEMP@_FillValue);"
                    -s "*_T2=array(0.0,0.0,/\$time,\$Nz2,\$$Ny,\$$Nx/);_T2.set_miss(TEMP@_FillValue);"
                    -s '_T1(:, :, :, :)=TEMPdz(:, 0:4, :, :);'
                    -s '_T2(:, :, :, :)=TEMPdz(:, 5:32, :, :);'
                    -s 'T1=_T1.ttl($Nz1)/dz1;'
                    -s 'T2=_T2.ttl($Nz2)/dz2;'
                    -s 'STRAT=T1-T2;'
                    $old_file $new_file3
        `; igs=false)

        println("Generating file: $(new_file4)")
        pleaseRun(`ncremap -R '--rgr lat_nm_in=Ny --rgr lon_nm_in=Nx' -m $(parsed["remap-file-nn"]) $new_file3 $new_file4`;igs=false)

    end

    # Ocean heat content
    new_file5 = joinpath(parsed["output-dir"], "$(parsed["casename"]).EMOM_extra3.$(date_str).nc")
    if !isfile(new_file5) || parsed["overwrite"]

        println("Generating file: $(new_file5)")

        if !parsed["pop2"]
            pleaseRun(`ncwa -a N1 -v TEMP,SALT,area_sT $old_file $new_file5`)
            pleaseRun(`ncwa -O -w area_sT -a Nx,Ny -v TEMP,SALT $new_file5 $new_file5`)
        else
            pleaseRun(`ncwa -w TAREA -a nlon,nlat -v TEMP,SALT $old_file $new_file5`)
        end

    end

    # Ocean heat transport
    new_file6 = joinpath(parsed["output-dir"], "$(parsed["casename"]).EMOM_extra4_OHT.$(date_str).nc")
    if !isfile(new_file6) || parsed["overwrite"]

        println("Generating file: $(new_file6)")

        if !parsed["pop2"]
            pleaseRun(`ncwa -a N1 -v ADVT,ADVS,WKRSTT,WKRSTS,Q_LOST,TEMP,SALT,dz_cT $old_file $new_file6`)
            pleaseRun(`ncwa -O -y ttl -w dz_cT -a Nz -v ADVT,ADVS,WKRSTT,WKRSTS,Q_LOST,TEMP,SALT $new_file6 $new_file6`)
#            pleaseRun(`ncap2 -O -s 'SHF=ADVT+WKRSTT;' $new_file6 $new_file6`)
        else
            pleaseRun(`ncwa -O -y ttl -w dz -a z_t -v TEMP,SALT $old_file $new_file6`)
            pleaseRun(`ncks -A -v SHF,N_HEAT $old_file $new_file6`)
        end

    end

end

max_stars_repo_name: meteorologytoday/CESM-diagnostic

max_issues_repo_name: meteorologytoday/CESM-diagnostic

max_forks_repo_name: meteorologytoday/CESM-diagnostic


text:
function solve_find_invalid(input, n)
	preamble = input[1:n]
	sums = zeros(Int64, n, n)
	for i in 1:n
		for j in i + 1:n
			sums[i,j] = preamble[i] + preamble[j]
		end
	end
	ret_i = n + 1
	while !isnothing(findfirst(x -> x == input[ret_i], sums))
		i = ret_i % n != 0 ? ret_i % n : n
		preamble[i] = input[ret_i]
		for j in filter(x -> x != i, 1:n)
			sums[i, j] = preamble[i] + preamble[j]
		end
		ret_i += 1
	end
	input[ret_i]
end

function solve_find_weakness(input, target)
	i = 2
	arr = input[1:2]
	while sum(arr) != target
		if sum(arr) < target
			i += 1
			push!(arr, input[i])
		else
			popfirst!(arr)
		end
	end
	minimum(arr) + maximum(arr)
end

@assert solve_find_invalid(parse.(Int64, readlines("input/test/test_day9.txt")), 5) == 127 ["test case 1 failed"]
@assert solve_find_weakness(parse.(Int64, readlines("input/test/test_day9.txt")), 127) == 62 ["test case 2 failed"]

invalid_number = solve_find_invalid(parse.(Int64, readlines("input/day9.txt")), 25)
println("First star: " * string(invalid_number))
println("Second star: " * string(solve_find_weakness(parse.(Int64, readlines("input/day9.txt")), invalid_number)))

ext: jl


text:
@kwargs function select!(
    op::SelectUnion,
    C::GBVecOrMat,
    A::GBVecMatOrTrans;
    thunk::Union{GBScalar, Nothing} = nothing
)
    thunk === nothing && (thunk = C_NULL)
    A, desc, _ = _handletranspose(A, desc)
    if A isa GBVector && C isa GBVector
        libgb.GxB_Vector_select(C, mask, accum, op, A, thunk, desc)
    elseif A isa GBMatrix && C isa GBMatrix
        libgb.GxB_Matrix_select(C, mask, accum, op, A, thunk, desc)
    end
end

@kwargs function select(
    op::SelectUnion,
    A::GBVecMatOrTrans;
    thunk::Union{GBScalar, Nothing} = nothing
)
    C = similar(A)
    select!(op, C, A; thunk, accum, mask, desc)
    return C
end


text:
using TensorOperations

"""
	h0::Real                # constant energy shift
	h1::Array{Float64,2}    # one electron integrals
	h2::Array{Float64,4}    # two electron integrals (chemist's notation)

Type to hold a second quantized Hamiltonian coefficients in memory
"""
struct InCoreInts
    h0::Float64
    h1::Array{Float64,2}
    h2::Array{Float64,4}
end

# @with_kw struct ElectronicProblem
# 	#=
# 	Structure to hold all problem specific parameters
# 	=#
# 	no::Int = 0
# 	na::Int = 0  # number of alpha
# 	nb::Int = 0  # number of beta
# 	basis::String = "sto-3g"
# 	dima::Int = calc_nchk(no,na)
# 	dimb::Int = calc_nchk(no,nb)
# 	#dim::Int = dima*dimb
# 	#converged::Bool = false
# 	#restarted::Bool = false
# 	#iteration::Int = 0
# 	#algorithm::String = "direct"    #  options: direct/davidson
# 	#n_roots::Int = 1
# end

"""
    id::Integer             #index of atom in the molecule
    symbol::String          #Atomic ID (E.g. H, He, ...)
    xyz::Array{Float64,1}   #list of XYZ coordinates

Simply an Atom
"""
struct Atom
    #=
    Type defining an atom
    =#
    id::Integer
    symbol::String
    xyz::Array{Float64,1}
end

"""
    charge::Integer         #overall charge on molecule
    multiplicity::Integer   #2S+1
    atoms::Vector{Atom}     #Vector of `Atoms`
    basis::String           #Basis set
Molecule essentially as a Vector of atoms, number of electrons and basis set
"""
struct Molecule
	charge::Integer
	multiplicity::Integer
	atoms::Array{Atom,1}
    basis::String
end

"""
	orbital_rotation!(ints::InCoreInts, U)

Transform electronic integrals, by U
i.e.,
```math
h_{pq} = U_{rp}h_{rs}U_{sq}
```
```math
(pq|rs) = (tu|vw)U_{tp}U_{uq}U_{vr}U_{ws}
```
"""
function orbital_rotation!(ints::InCoreInts, U)
	@tensor begin
		ints.h1[p,q] = U[r,p]*U[s,q]*ints.h1[r,s]
		ints.h2[p,q,r,s] = U[t,p]*U[u,q]*U[v,r]*U[w,s]*ints.h2[t,u,v,w]
	end
end

@doc raw"""
	orbital_rotation(ints::InCoreInts, U)

Transform electronic integrals, by U
i.e.,
```math
h_{pq} = U_{rp}h_{rs}U_{sq}
```
```math
(pq|rs) = (tu|vw)U_{tp}U_{uq}U_{vr}U_{ws}
```
"""
function orbital_rotation(ints::InCoreInts, U)
    @tensor begin
        h1[p,q] := U[r,p]*U[s,q]*ints.h1[r,s]
        # h2[p,q,r,s] := U[t,p]*U[u,q]*U[v,r]*U[w,s]*ints.h2[t,u,v,w]
        h2[p,q,r,s] := U[t,p]*ints.h2[t,q,r,s]
        h2[p,q,r,s] := U[t,q]*h2[p,t,r,s]
        h2[p,q,r,s] := U[t,r]*h2[p,q,t,s]
        h2[p,q,r,s] := U[t,s]*h2[p,q,r,t]
    end
    return InCoreInts(ints.h0,h1,h2)
end

"""
    subset(ints::InCoreInts, list; rmd1a=nothing, rdm1b=nothing)

Extract a subset of integrals acting on orbitals in list, returned as `InCoreInts` type
- `ints::InCoreInts`: Integrals for full system
- `list`: list of orbital indices in subset
- `rdm1a`: 1RDM for embedding α density to make CASCI hamiltonian
- `rdm1b`: 1RDM for embedding β density to make CASCI hamiltonian
"""
function subset(ints::InCoreInts, list, rdm1a=nothing, rdm1b=nothing)
    ints_i = InCoreInts(ints.h0, view(ints.h1,list,list), view(ints.h2,list,list,list,list))
    if rdm1b != nothing
        if rdm1a == nothing
            throw(Exception)
        end
    end

    if rdm1a != nothing
        if rdm1b == nothing
            throw(Exception)
        end
        da = deepcopy(rdm1a)
        db = deepcopy(rdm1b)
        da[:,list] .= 0
        db[:,list] .= 0
        da[list,:] .= 0
        db[list,:] .= 0
        viirs = ints.h2[list, list,:,:]
        viqri = ints.h2[list, :, :, list]
        fa = zeros(length(list),length(list))
        fb = copy(fa)
        @tensor begin
            ints_i.h1[p,q] += viirs[p,q,r,s] * (da+db)[r,s]
            ints_i.h1[p,s] -= .5*viqri[p,q,r,s] * da[q,r]
            ints_i.h1[p,s] -= .5*viqri[p,q,r,s] * db[q,r]
        end
    end
    return ints_i
end

"""
	compute_energy(h0, h1, h2, rdm1, rdm2)

Given an energy shift `h0`, 1e integrals `h1`, and 2e ints `h2`
along with a 1rdm and 2rdm on the same space, return the energy
"""
function compute_energy(h0, h1, h2, rdm1, rdm2)
    e = h0
    e += sum(h1 .* rdm1)
    e += .5*sum(h2 .* rdm2)
    # @tensor begin
    # 	e  += .5 * (ints.h2[p,q,r,s] * rdm2[p,q,r,s])
    # end
    return e
end
"""
	compute_energy(ints::InCoreInts, rdm1, rdm2)

Return energy defined by `rdm1` and `rdm2`
"""
function compute_energy(ints::InCoreInts, rdm1, rdm2)
    e = ints.h0
    e += sum(ints.h1 .* rdm1)
    e += .5*sum(ints.h2 .* rdm2)
    # @tensor begin
    # 	e  += .5 * (ints.h2[p,q,r,s] * rdm2[p,q,r,s])
    # end
    return e
end

max_stars_repo_name: nmayhall-vt/FermiCG

max_issues_repo_name: nmayhall-vt/FermiCG

max_forks_repo_name: nmayhall-vt/FermiCG


text:
"""
    fit_mvnormals(θs, ∇logpθs; history_length=5)

Fit a multivariate-normal distribution to each point on the trajectory `θs`.

Given `θs` with gradients `∇logpθs`, construct LBFGS inverse Hessian approximations with
the provided `history_length`. The inverse Hessians approximate a covariance. The
covariances and corresponding means that define multivariate normal approximations per
point are returned.
"""
function fit_mvnormals(θs, ∇logpθs; history_length, kwargs...)
    if length(θs) < history_length
        error("Insufficient history in trajectory.")
    end
    Σs = lbfgs_inverse_hessians(θs, ∇logpθs; kwargs...)
    l = length(Σs)
    μs = @views muladd.(Σs, ∇logpθs[1:l], θs[1:l])
    return Distributions.MvNormal.(μs, Σs)
end

# faster than computing `logpdf` and `rand` independently
function rand_and_logpdf(rng, dist::Distributions.MvNormal, ndraws)
    μ = dist.μ
    Σ = dist.Σ
    N = length(μ)

    # draw points
    u = Random.randn!(rng, similar(μ, N, ndraws))
    unormsq = vec(sum(abs2, u; dims=1))
    x = PDMats.unwhiten!(u, Σ, u)
    x .+= μ

    # compute log density at each point
    logpx = (muladd(N, log2π, logdet(Σ)) .+ unormsq) ./ -2

    return x, logpx
end


text:
using Test
using CrypticCrosswords

@testset "Derivations" begin
    clue = "aerial worker anne on the way up"
    pattern = r"^.{7}$"
    solutions, state = CC.solve(clue, pattern=pattern, min_grammar_score=2e-5)
    sol = first(solutions)
    derivations = derive!(state, sol)
    @test length(derivations) == 1
    d = first(derivations)
    @test d.output == "antenna"
    @test d.derivation.constituents[1].output == "aerial"
    @test d.derivation.constituents[2].output == "antenna"
    @test d.derivation.constituents[2].constituents[1].output == "ant"
    @test d.derivation.constituents[2].constituents[1].constituents[1].output == "ant"
    @test d.derivation.constituents[2].constituents[1].constituents[1].constituents[1].output == "worker"
    @test d.derivation.constituents[2].constituents[2].output == "enna"
    @test d.derivation.constituents[2].constituents[2].constituents[1].output == "enna"
    @test d.derivation.constituents[2].constituents[2].constituents[1].constituents[1].output == "anne"
end


text:
function trimolecular_1path_diff!(du,u,p,t,dens,cons)
    # unpack rates and constants
    nᵣ,nᵣₓ,n₃ = u
    kₓ,k₋ₓ,k₃,k₋₃ = p
    k₁,k₋₁ = cons
    mᵣ,mₗ,mₓ,A = dens

    # model
    du[1] = dnᵣ  = A*k₁*mᵣ*mₗ - k₋₁*nᵣ - A*kₓ*mₓ*nᵣ + k₋ₓ*nᵣₓ
    du[2] = dnᵣₓ = A*kₓ*mₓ*nᵣ - k₋ₓ*nᵣₓ - k₃*nᵣₓ + k₋₃*n₃
    du[3] = dn₃  = k₃*nᵣₓ - k₋₃*n₃

end


text:
# This file is a part of MGVI.jl, licensed under the MIT License (MIT).

function fisher_information_and_jac(
    fwd_model::Function, ξ::AbstractVector,
    jac_method::Type{JF}
) where JF<:AbstractJacobianFunc
    ℐ_λ = fisher_information(fwd_model(ξ))
    dλ_dξ= jac_method(flat_params ∘ fwd_model)(ξ)
    ℐ_λ, dλ_dξ
end

const rs_default_options=NamedTuple()
const optim_default_options = Optim.Options()
const optim_default_solver = LBFGS()

function _create_residual_sampler(f::Function, center_p::Vector;
                                  residual_sampler::Type{RS}=ImplicitResidualSampler,
                                  jacobian_func::Type{JF}=FwdDerJacobianFunc,
                                  residual_sampler_options::NamedTuple
                                 ) where RS <: AbstractResidualSampler where JF <: AbstractJacobianFunc
    ℐ_λ, dλ_dξ = fisher_information_and_jac(f, center_p, jacobian_func)
    residual_sampler(ℐ_λ, dλ_dξ; residual_sampler_options...)
end

function posterior_loglike(model, p, data)
    logpdf(model(p), data) - dot(p, p)/2
end

function mgvi_kl(f::Function, data, residual_samples::AbstractMatrix{<:Real}, center_p::AbstractVector{<:Real})
    kl_component(p::AbstractVector) = -posterior_loglike(f, p, data)
    kl_comp_both(rs::AbstractVector) = kl_component(center_p + rs) + kl_component(center_p - rs)
    res = 0
    for rs in eachcol(residual_samples)
        res += kl_comp_both(rs)
    end
    res/size(residual_samples, 2)/2
end

function _gradient_for_optim(kl::Function)
    (res::AbstractVector, x::AbstractVector) -> begin
        res .= Zygote.gradient(kl, x)[1]
    end
end

"""
    mgvi_kl_optimize_step(rng, forward_model, data, init_param_point;
                          jacobian_func=jacobian_func,
                          residual_sampler=residual_sampler,
                          [num_residuals=3,]
                          [residual_sampler_options=NamedTuple(),]
                          [optim_solver=LBFGS(),]
                          [optim_options=Optim.Options()])

Performs one MGVI iteration.

The posterior distribution is approximated with a multivariate normal distribution.
The covariance is approximated with the inverse Fisher information valuated at
`init_param_point`. Samples are drawn according to this covariance, which are then
used to estimate and minimize the KL divergence between the true posterior and the
approximation.

# Arguments

* `rng::AbstractRNG`: instance of the random number generator
* `forward_model::Function`: turns model parameters into an instance of Distribution
* `data::AbstractVector`: data on which model's pdf is evaluated
* `init_param_point::Vector`: initial estimate of the model parameters
* `jacobian_func::Type{<:AbstractJacobianFunc}`: method to calculate the Jacobian
  of the forward model
* `residual_sampler::Type{<:AbstractResidualSampler}`: method to draw samples from
  the approximation
* `num_residuals::Integer = 3`: number of samples used to estimate the KL divergence
* `residual_sampler_options::NamedTuple = NamedTuple()`: further options to pass to the
  residual sampler. Important is `cg_params` that is passed to the CG solver used inside of ImplicitResidualSampler
* `optim_solver::Optim.AbstractOptimizer = LBFGS()`: optimizer used for minimizing KL divergence
* `optim_options::Optim.Options = Optim.Options()`: options to pass to the KL optimizer

# Example

```julia
using Random, Distributions, MGVI

model(x::AbstractVector) = Normal(x[1], 0.2)
true_param = [2.0]
data = rand(model(true_param), 1)[1]
init_param = [1.3]

res = mgvi_kl_optimize_step(Random.GlobalRNG, model, data, init_param;
                            jacobian_func=FwdRevADJacobianFunc,
                            residual_sampler=ImplicitResidualSampler,
                            num_residuals=5,
                            residual_sampler_options=(;cg_params=(;maxiter=10, verbose=true))),
                            optim_solver=LBFGS(;m=5),
                            optim_options=Optim.Options(iterations=7, show_trace=true))

next_param_point = res.result

optim_optimized_object = res.optimized
Optim.summary(optim_optimized_object)

samples_from_est_covariance = res.samples
```

# See also

* Residual samplers: [`AbstractResidualSampler`](@ref), [`ImplicitResidualSampler`](@ref), [`FullResidualSampler`](@ref)
* Jacobian functions: [`AbstractJacobianFunc`](@ref), [`FwdRevADJacobianFunc`](@ref), [`FwdDerJacobianFunc`](@ref)

"""
function mgvi_kl_optimize_step(rng::AbstractRNG,
                               forward_model::Function, data, init_param_point::AbstractVector;
                               jacobian_func::Type{JF},
                               residual_sampler::Type{RS},
                               num_residuals::Integer=3,
                               residual_sampler_options::NamedTuple=rs_default_options,
                               optim_solver::Union{Optim.AbstractOptimizer, NewtonCG}=optim_default_solver,
                               optim_options::Union{Optim.Options, Nothing}=optim_default_options
                              ) where RS <: AbstractResidualSampler where JF <: AbstractJacobianFunc
    est_res_sampler = _create_residual_sampler(forward_model, init_param_point;
                                               residual_sampler=residual_sampler,
                                               jacobian_func=jacobian_func,
                                               residual_sampler_options=residual_sampler_options)
    residual_samples = rand(rng, est_res_sampler, num_residuals)
    kl(params::AbstractVector) = mgvi_kl(forward_model, data, residual_samples, params)
    ∇kl! =  _gradient_for_optim(kl)
    Σ⁻¹(ξ) = _inv_cov_est(forward_model, ξ, jacobian_func)
    Σ̅⁻¹(ξ) = mean(Σ⁻¹.(collect.(eachcol(ξ .+ residual_samples))))
    res = _optimize(
        kl, ∇kl!, Σ̅⁻¹, init_param_point, optim_solver, optim_options)
    updated_p = res.minimizer

    (result=updated_p, optimized=res, samples=hcat(updated_p .+ residual_samples, updated_p .- residual_samples))
end

export mgvi_kl_optimize_step

function _inv_cov_est(fwd_model::Function, ξ::AbstractVector, jac_method::Type{JF}) where JF <: AbstractJacobianFunc
    ℐ_λ, dλ_dξ = fisher_information_and_jac(fwd_model, ξ, jac_method)
    dλ_dξ' * ℐ_λ * dλ_dξ + I
end


text:
#-----------------------------------------------------------------------# DataValueArray

function mmwrite(io::AbstractSerializer, xs::DataValueArray)
    Serialization.serialize_type(io, MMSer{DataValueArray})

    mmwrite(io, BitArray(xs.isna))
    mmwrite(io, xs.values)
end

function mmread(::Type{DataValueArray}, io, mmap)
    isnull = deserialize(io)
    vals = deserialize(io)
    DataValueArray(vals, isnull)
end

#-----------------------------------------------------------------------# PooledArray
function mmwrite(io::AbstractSerializer, xs::PooledArray)
    Serialization.serialize_type(io, MMSer{PooledArray})

    mmwrite(io, xs.refs)
    mmwrite(io, xs.invpool)
end

function mmread(::Type{PooledArray}, io, mmap)
    refs = deserialize(io)
    invpool = deserialize(io)
    PooledArray(PooledArrays.RefArray(refs), invpool)
end

#-----------------------------------------------------------------------# Columns

function mmwrite(io::AbstractSerializer, xs::Columns)
    Serialization.serialize_type(io, MMSer{Columns})

    if eltype(xs) <: NamedTuple
        fnames = fieldnames(eltype(xs))
    else
        fnames = length(columns(xs))
    end

    serialize(io, fnames)
    for c in columns(xs)
        mmwrite(io, c)
    end
end

function mmread(::Type{Columns}, io, mmap)
    fnames = deserialize(io)
    if isa(fnames, Int)
        cols = [deserialize(io) for i=1:fnames]
        Columns(Tuple(cols))
    else
        cols = [deserialize(io) for i=1:length(fnames)]
        Columns(Tuple(cols); names=fnames)
    end
end

#-----------------------------------------------------------------------# NDSparse
function mmwrite(io::AbstractSerializer, xs::NDSparse)
    Serialization.serialize_type(io, MMSer{NDSparse})

    flush!(xs)
    mmwrite(io, xs.index)
    mmwrite(io, xs.data)
end

function mmread(::Type{NDSparse}, io, mmap)
    idx = deserialize(io)
    data = deserialize(io)
    NDSparse(idx, data, presorted=true, copy=false)
end

#-----------------------------------------------------------------------# IndexedTable
function mmwrite(io::AbstractSerializer, xs::IndexedTable)
    Serialization.serialize_type(io, MMSer{IndexedTable})

    #flush!(xs)
    mmwrite(io, rows(xs))
    mmwrite(io, xs.pkey)
    mmwrite(io, xs.perms)
end

function mmread(::Type{IndexedTable}, io, mmap)
    data = deserialize(io)
    pkey = deserialize(io)
    perms = deserialize(io)
    table(data, pkey=pkey, perms=perms, presorted=true, copy=false)
end


text:
module Synthetic

export gen_blob

using Celeste, Celeste.Model
import Celeste: Infer, DeterministicVI

import WCS
import Distributions

include("../src/bivariate_normals.jl")

# Generate synthetic data.

function wrapped_poisson(rate::Float64)
    0 < rate ? float(rand(Distributions.Poisson(rate))) : 0.
end

function get_patch(the_mean::SVector{2,Float64}, H::Int, W::Int)
    radius = 50
    hm = round(Int, the_mean[1])
    wm = round(Int, the_mean[2])
    w11 = max(1, wm - radius):min(W, wm + radius)
    h11 = max(1, hm - radius):min(H, hm + radius)
    return (w11, h11)
end

function write_gaussian(the_mean, the_cov, intensity, pixels;
                        expectation=false)
    the_precision = inv(the_cov)
    c = sqrt(det(the_precision)) / 2pi

    H, W = size(pixels)
    w_range, h_range = get_patch(the_mean, H, W)

    for w in w_range, h in h_range
        y = @SVector [the_mean[1] - h, the_mean[2] - w] # Maybe not hard code Float64
        ypy = dot(y,  the_precision * y)
        pdf_hw = c * exp(-0.5 * ypy)
        pixel_rate = intensity * pdf_hw
        pixels[h, w] += expectation ? pixel_rate : wrapped_poisson(pixel_rate)
    end

    pixels
end

function write_star(img0::Image, ce::CatalogEntry, pixels::Matrix{Float64};
                    expectation=false)
    iota = median(img0.iota_vec)
    for k in 1:length(img0.psf)
        the_mean = WCS.world_to_pix(img0.wcs, ce.pos) + img0.psf[k].xiBar
        the_cov = img0.psf[k].tauBar
        intensity = ce.star_fluxes[img0.b] * iota * img0.psf[k].alphaBar
        write_gaussian(the_mean, the_cov, intensity, pixels,
            expectation = expectation)
    end
end

function write_galaxy(img0::Image, ce::CatalogEntry, pixels::Matrix{Float64};
                      expectation=false)
    iota = median(img0.iota_vec)
    e_devs = [ce.gal_frac_dev, 1 - ce.gal_frac_dev]
    XiXi = DeterministicVI.get_bvn_cov(ce.gal_ab, ce.gal_angle, ce.gal_scale)

    for i in 1:2
        for gproto in galaxy_prototypes[i]
            for k in 1:length(img0.psf)
                the_mean = WCS.world_to_pix(img0.wcs, ce.pos) +
                           img0.psf[k].xiBar
                the_cov = img0.psf[k].tauBar + gproto.nuBar * XiXi
                intensity = ce.gal_fluxes[img0.b] * iota *
                    img0.psf[k].alphaBar * e_devs[i] * gproto.etaBar
                write_gaussian(the_mean, the_cov, intensity, pixels,
                    expectation=expectation)
            end
        end
    end
end

function gen_image(img0::Image, n_bodies::Vector{CatalogEntry}; expectation=false)
    epsilon = img0.epsilon_mat[1]
    iota = img0.iota_vec[1]

    if expectation
        pixels = [epsilon * iota for h=1:img0.H, w=1:img0.W]
    else
        pixels = reshape(float(rand(Distributions.Poisson(epsilon * iota),
                         img0.H * img0.W)), img0.H, img0.W)
    end

    for body in n_bodies
        body.is_star ? write_star(img0, body, pixels) : write_galaxy(img0, body, pixels)
    end

    epsilon_mat = fill(epsilon, img0.H, img0.W)
    iota_vec = fill(iota, img0.H)

    return Image(img0.H, img0.W, pixels, img0.b, img0.wcs,
                 img0.psf,
                 img0.run_num, img0.camcol_num, img0.field_num,
                 epsilon_mat, iota_vec,
                 img0.raw_psf_comp)
end

"""
Generate a simulated blob based on a vector of catalog entries using
identity world coordinates.
"""
function gen_blob(blob0::Vector{Image}, n_bodies::Vector{CatalogEntry}; expectation=false)
    [gen_image(blob0[b], n_bodies, expectation=expectation) for b in 1:5]
end

#######################################

const pp = Model.load_prior()

function sample_fluxes(i::Int, r_s)
#    r_s = rand(Distributions.Normal(pp.r_mean[i], pp.r_var[i]))
    k_s = rand(Distributions.Categorical(pp.k[i]))
    c_s = rand(Distributions.MvNormal(pp.c[i][:, k_s], pp.c[i][:, :, k_s]))

    l_s = Array(Float64, 5)
    l_s[3] = r_s
    l_s[4] = l_s[3] * exp(c_s[3])
    l_s[5] = l_s[4] * exp(c_s[4])
    l_s[2] = l_s[3] / exp(c_s[2])
    l_s[1] = l_s[2] / exp(c_s[1])
    l_s
end

function synthetic_body(ce::CatalogEntry)
    ce2 = deepcopy(ce)
#    ce2.is_star = rand(Distributions.Bernoulli(pp.a[1]))
    ce2.star_fluxes[:] = sample_fluxes(1, ce.star_fluxes[3])
    ce2.gal_fluxes[:] = sample_fluxes(2, ce.gal_fluxes[3])
    ce2
end

function synthetic_bodies(n_bodies::Vector{CatalogEntry})
    CatalogEntry[synthetic_body(ce) for ce in n_bodies]
end

end


text:
##########
###Photons
##########

###Sequence of uniform random numbers
photonsequence(r::AbstractVector,len::Int,timescale =0.001) = GeneralizedMetropolisHastings.data(:array,timescale*(1:len),rand(r,len))

###Load an existing photon sequence
function photonsequence(filename::AbstractString,timescale =0.001f0)
    photons = trunc(Int,JLD.load(filename,"photons"))
    GeneralizedMetropolisHastings.data(:array,timescale*(1:length(photons)),photons)
end

##########
###Current
##########

###Load an existing light-induced current sequence
function lightinducedcurrent(filename::AbstractString,timescale =0.001f0)
    lic = JLD.load(filename,"lightinducedcurrent")
    GeneralizedMetropolisHastings.data(:array,timescale*(1:length(lic)),lic)
end


text:
function foo()
    lines = readlines(open("data/6.txt", "r"))
    total_distance = 10000
    points = zeros(Int, length(lines), 2)
    maxx = 0
    for (i, line) in enumerate(lines)
        point = parse.(Int, split(line, ','))
        if any(point .> maxx)
            maxx = maximum(point)
        end
        points[i, :] = point
    end

    stop = maxx + 10
    start = -10
    area = 0
    for i in start:stop
        for j in start:stop
            dist = abs.(points[:, 1] .- i) .+ abs.(points[:, 2] .- j)
            area += sum(dist) < total_distance
        end
    end
    return area
end

println(foo())


text:
#=

=#
function sortededgesbyweight(G::DictGraph)
    @assert !G.directed "the graph must be undirected"
    𝑄 = PriorityQueue{Tuple{Int,Int}, Int}()
    for k in keys(G.relation)
        for link in G.relation[k]
            local t = (k, link.u)
            (!haskey(𝑄, t) && !haskey(𝑄, reverse(t))) && push!(𝑄, t => link.w)
        end
    end
    return 𝑄

end

function mstkruskal(G::DictGraph)
    T = DictGraph()
    Q = sortededgesbyweight(G)
    DS = IntDisjointSets(length(G.relation))

    while !isempty(Q)
        e = dequeue_pair!(Q)
        if !in_same_set(DS, e.first[1] + 1, e.first[2] + 1)
            union!(DS, e.first[1] + 1, e.first[2] + 1)
            addedge!(T, e.first[1], e.first[2]; weight = e.second)
        end
    end

    return T
end


text:
module NLSolversTests

using Gridap.Algebra
using Gridap.Algebra: NonlinearOperatorMock

op = NonlinearOperatorMock()
nls = NLSolver(show_trace=false,method=:newton)

x0 = zero_initial_guess(op)
x = [1.0, 3.0]
test_nonlinear_solver(nls,op,x0,x)

x0 = [2.1,2.9]
x = [2.0, 3.0]
test_nonlinear_solver(nls,op,x0,x)

x0 = zero_initial_guess(op)
cache = solve!(x0,nls,op)

_ = cache.result

end # module


text:
#solve eq x=cos(x)
#Reference: http://nbviewer.jupyter.org/gist/genkuroki/ce7f62d10365b9a25c82aad499c886c2

function bisection(f, x₀::T, x₁::T; ε=eps(T), maxiter=10^8) where T<:AbstractFloat
    y₀, y₁ = f(x₀), f(x₁)
    iter = 0
    y₀ == 0 && return x₀, iter
    y₁ == 0 && return x₁, iter
    y₀ * y₁ > 0 && return T(NaN), -1
    for iter in 1:maxiter
        x₂ = (x₀ + x₁)/2
        y₂ = f(x₂)
        y₂ == 0 && return x₂, iter
        if y₀ * y₂ > 0
            x₀, y₀ = x₂, y₂
        else
            x₁, y₁ = x₂, y₂
        end
        abs(x₁ - x₀) < ε && return x₂, iter
    end
    return T(NaN), iter
end

function main()
    f(x)=x-cos(x)
    x₀, x₁ = 0.0, 3.14
    @show α, iter = bisection(f, x₀, x₁)
    @show f(α)
end

main()

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
module Modes

using Unitful: ustrip, @u_str
using LinearAlgebra: norm, ⋅, ×
using StaticArrays: SVector
using ..Constants: twopi, η₀, dbmin
using ..Layers: Layer, TEorTM, TE, TM
using ..Sheets: RWGSheet, find_unique_periods
using ..GSMs: Gblock
using ..Rings: Ring
using ..PGF: mysqrt
using ..Log: @logfile

const MNmax_default = 6

"""
    choose_layer_modes!(layers, sheets, junc, gbl, k0max, dbmin)

    #
Set up the arrays of modal indices `M`, `N`, and `P` for all layers not included
in `Gblock`s given the wavenumber at the maximum operating frequency `k0max`,
the list of layers and sheets `strata`, the list of `Gblock`s `gbl`, and the
desired minimum attenuation `dbmin` for each neglected mode to encounter
when passing through the layer.  Also, allocate the arrays `β`, `γ`, `Y`,
`c`, and `tvec` for each of the excluded layers.  Note that the modes are not
necessarily defined consistently in each layer because the periodicity of the
multiple FSS/PSS sheets may not all be identical.  Therefore, we store the β₁
and β₂ reciprocal lattice vectors for each layer in the `Layer` type and use
these values to define the periodicity for the modes in a given layer.

## Arguments

- `layers`: A vector of Layer objects. It is assumed that for each `Layer` not included in a `Gblock`,
the permeability and permittivity have been correctly initialized.  On exit, these same excluded layers
will have fields `β₁` and `β₂` appropriately set, and will have the arrays `β`, `γ`, `Y`,
`c`, and `tvec` allocated for the  appropriate number of modes.

- `sheets`: A vector of Sheet objects.

- `junc`: An integer vector.  `junk[k]` contains the sheet number located at junction `k` or `0` if
  no sheet is present.

- `gbl`: An collection of `Gblock`s containing the definitions of the GSM block entities.

- `k0max`: Maximum free-space wavenumber to be analyzed (1/m).

- `dbmin`  Minimum attenuation any neglected mode must incur (dB > 0).
"""
function choose_layer_modes!(layers::Vector{Layer}, sheets::Vector{RWGSheet}, junc::Vector{Int},
                                                 gbl::Vector{Gblock}, k0max::Float64, dbmin::Float64=dbmin)

    nl = length(layers)
    nj = nl-1 # number of dielectric junctions
    ns = length(sheets)
    ngbl = length(gbl)

    # Make all layers unique:
    for i in 2:length(layers)
        for j in 1:i-1
            if layers[i] === layers[j]
                layers[i] = deepcopy(layers[j])
                @goto nexti
            end
        end
        @label nexti
    end

    # Possible mode indices m and n will range from -MNmax to MNmax,for a total of 2(2MNmax+1)^2
    # modes in the layer (first factor of 2 is for both TE/TM).
    mset = falses(nl)  # Indicates that modes have not been determined yet for any layers
    for g in gbl  # Don't need to set modes for layers included in a Gblock
        mset[g.rng[2:end]] .= true
    end

    pure_radome = true
    for isht in junc
      if isht ≠ 0 && sheets[isht].style ≠ "NULL"
        pure_radome = false
        break
      end
    end

    if pure_radome
        for l in layers
            l.P = [TE, TM]
            l.M = [0, 0]
            l.N = [0, 0]
            l.β₁ = SVector(twopi, 0.0)
            l.β₂ = SVector(0.0, twopi) # cell 1m in x by 1m in y.
        end
        mset .= true
        notinablock = collect(1:length(layers))::Vector{Int}

    else   # There is at least one non-NULL sheet in the structure
        last_gbl = findlast(x -> !iszero(x.j) && sheets[junc[x.j]].style ≠ "NULL", gbl)
        last_sheet = junc[gbl[last_gbl].j] # sheet index
        first_gbl = findfirst(x -> !iszero(x.j) && sheets[junc[x.j]].style ≠ "NULL", gbl)
        first_sheet = junc[gbl[first_gbl].j] # sheet index
        # Assign equivalence class numbers to the interfaces based on the
        # unit cell of the sheet located there:
        upa = find_unique_periods(junc, sheets)
        # Treat all unassigned layers that are adjacent to a gbl containing an FSS screen:
        for igbl in 1:ngbl
            g = gbl[igbl]
            ijunc = g.j # Junction # where FSS is located
            ijunc == 0 && continue
            isht = junc[ijunc]  # Sheet number at junction ijunc
            sheet = sheets[isht]
            sheet.style == "NULL" && continue
            one_meter = ustrip(Float64, sheet.units, 1u"m")
            β₁ = sheet.β₁ * one_meter # radians/meter
            β₂ = sheet.β₂ * one_meter # radians/meter
            i1 = g.rng[1] # Index of layer to left of current GBLOCK.
            i2 = 1 + g.rng[end] # Index of layer to right of current GBLOCK.
            for i in [i1, i2]
                (i==1 || i==nl) && continue  # Skip boundary layers.
                # Determine neighboring gblock to the outside:
                if i == i1
                    igbl_other = max(igbl-1, 1)
                else
                    igbl_other = min(igbl+1, ngbl)
                end
                ijunc_other = gbl[igbl_other].j
                if ijunc_other ≠ 0
                    isht_other = junc[ijunc_other]  # Sheet number at junction ijunc_other
                else
                    isht_other = 0
                end
                if !mset[i]  # This layer not yet assigned
                    layers[i].β₁ = copy(β₁)  # Save periodicity in layer i.
                    layers[i].β₂ = copy(β₂)  # Save periodicity in layer i.
                    MNmax = MNmax_default
                    # These values must be changed to 0 if there is also another GBLOCK
                    # adjacent to layer i containing a nonnull sheet with a different periodicity:
                    if ijunc_other ≠ 0 && upa[ijunc] ≠ upa[ijunc_other] &&
                       sheets[isht_other].style ≠ "NULL"
                        @logfile    """
                                    ******************* Warning ***********************
                                       Unequal unit cells in sheets $(isht) and $(isht_other)
                                       Setting #modes in dividing layer $(i) to 2
                                    ******************* Warning ***********************

                                    """
                        MNmax = 0
                    end
                    # Select the mode indices for layer i:
                    fill_mmax_pmn!(layers[i], MNmax, k0max, dbmin)
                    mset[i] = true
                else
                    # mset is true, indicating that this layer belongs to a gblock containing an FSS.
                    # IF the periodicity of that FSS is inconsistent with that of the current gblock
                    # then we can't satisfy both requirements, so abort.
                    if ijunc_other ≠ 0 && upa[ijunc] ≠ upa[ijunc_other] && length(layers[i].P) ≠ 2 &&
                        sheets[isht_other].style ≠ "NULL"
                        error("""
                        Unequal unit cells in sheets $(isht) and $(isht_other).
                               Try dividing layer $(i) into 2 halves. """)
                    end
                end
            end
        end # for

        # At this point, all layers adjacent to a Gblock containing an FSS sheet have
        # been assigned modes.  Now we assign mode index values to the remaining layers
        # by stepping inwards from the already assigned layers. Loop until all
        # layers are assigned.
        while any(.!mset[2:nl-1])
            for i in 2:nl-1  # Step over each layer
                mset[i] && continue
                if mset[i-1] && mset[i+1]
                    # Both adjacent layers have been assigned.
                    if equal_layer_periodicity(layers[i-1],layers[i+1])
                        # The two adjacent layers have the same periodicity.  Choose
                        # the accessible modes from among the two sets of adjacent ones.
                        layers[i].β₁ = copy(layers[i+1].β₁)
                        layers[i].β₂ = copy(layers[i+1].β₂)
                        layers[i].P = vcat(layers[i-1].P, layers[i+1].P)
                        layers[i].M = vcat(layers[i-1].M, layers[i+1].M)
                        layers[i].N = vcat(layers[i-1].N, layers[i+1].N)
                        select_pmn!(layers[i], k0max, dbmin)
                        mset[i] = true
                    else
                        # Both adjacent layers have been assigned, but they have different
                        # periodicities.  If they each have only two modes (TE/TM), then
                        # all is well, since these m=n=0 modes are independent of the
                        # periodicity.  If not, we must truncate the modes to only these
                        # two, and print a warning message:
                        layers[i].β₁ = copy(layers[i+1].β₁) # Assign same periodicity.
                        layers[i].β₂ = copy(layers[i+1].β₂)
                        layers[i].P = [TE, TM]
                        layers[i].M = [0, 0]
                        layers[i].N = [0, 0]
                        mset[i] = true
                        if length(layers[i-1].P) ≠ 2 || length(layers[i+1].P) ≠ 2
                            @logfile """
                            ******************* Warning ***********************
                                Setting # modes in Layer $(i) to 2 due to
                                unequal unit cells in surrounding FSS sheets
                            ******************* Warning ***********************

                            """
                        end
                    end
                elseif  mset[i-1] && !mset[i+1]
                    # Only layer to the left has been assigned.
                    i-1 == 1 && @goto nextouter # Don't rely on end layer.
                    copy_βs_pmn!(layers[i], layers[i-1])
                    select_pmn!(layers[i], k0max, dbmin)
                    mset[i] = true
                elseif !mset[i-1] && mset[i+1]
                    # Only layer to the right has been assigned.
                    i+1 == nl && @goto nextouter # Don't rely on end layer.
                    copy_βs_pmn!(layers[i], layers[i+1])
                    select_pmn!(layers[i], k0max, dbmin)
                    mset[i] = true
                end # if
            end # for
            @label nextouter
        end # while
        # Done with all interior layers.
        # Choose two modes in the end layers (the principal Floquet modes).
        # Use the periodicity of the closest FSS sheets for each end layer.
        for (l, sh) in ((layers[1], sheets[first_sheet]), (layers[nl], sheets[last_sheet]))
            one_meter = ustrip(sh.units, 1u"m")
            l.β₁ = sh.β₁ * one_meter
            l.β₂ = sh.β₂ * one_meter
            fill_mmax_pmn!(l, 0, k0max, dbmin)
        end
        mset[1] = mset[nl] = true

        # Ensure that any adjacent non-gblock layers either have the same
        # periodicity or have only dominant modes:
        inagblock = reduce(union, (g.rng[2:end] for g in gbl)) # layer indices in Gblocks
        notinablock = sort(setdiff(1:nl, inagblock))
        for i in 1:length(notinablock)-1
            notinablock[i+1] - notinablock[i] == 1 || continue
            il = notinablock[i] # layers[il] and layers[il+1] are adjacent
            equal_layer_periodicity(layers[il], layers[il+1]) && continue
            for l in layers[il:il+1]
                l.P = [TE,TM]
                l.M = [0, 0]
                l.N = [0, 0]
            end
            @logfile """
            ******************* Warning ***********************
                Setting # modes in layers $(il) and $(il+1) to 2 due
                to unequal unit cells in surrounding FSS sheets
            ******************* Warning ***********************

            """
        end
    end # if

    # Now we allocate the other Layer field vectors
    maxmodes = 2*(2*MNmax_default+1)^2
    for i in notinablock
        layer = layers[i]
        nmodes = length(layer.P)
        for field in (:β, :γ, :Y, :c, :tvec)
            setfield!(layer, field, zeros(eltype(getfield(layer,field)), nmodes))
        end
        nmodes == maxmodes &&  @logfile """
        ******************* Warning ************************
        Maximum number of modes reached. GSM results may not
        be accurate.  Consider increasing MNmax_default from
        its current value of $(MNmax_default)
        ******************* Warning ************************

        """
    end

    return nothing
end

"""
    setup_modes!(layer::Layer, k0::Real, kvec::AbstractVector)

Fill the modal layer fields.  Needed for layers not contained in a Gblock.
The arrays are assumed to have  been already allocated, and the index arrays
`M`, `N`, and `P` are asssumed to have been already initialized.

### Input arguments

- `layer`: It is assumed that the `ϵᵣ`, `μᵣ`, `β₁`, and `β₂` have been
    initialized, as have the modal index arrays `P`, `M`, and `N`.
- `k0`: Free-space wavenumber (1/meter).
- `kvec`: A real-valued 2-vector containing the x and y components of the incident
    plane wave unit vector that defines the unit cell incremental
    phase shifts.

### Outputs

There is no explicit output, but the fields of `layer` will be modified, including
`β`, `tvec`, `γ`, `c`, and `Y`.
"""
function setup_modes!(layer::Layer, k0::Real, kvec::AbstractVector)
    β₀₀ = SVector(kvec[1], kvec[2])
    β₁, β₂ = layer.β₁, layer.β₂
    area = twopi^2 / norm(β₁ × β₂)
    ksq = k0^2 * layer.ϵᵣ * layer.μᵣ
    for mode in 1:length(layer.M)
        m,n,p = layer.M[mode], layer.N[mode], layer.P[mode]
        β = β₀₀ + m*β₁ + n*β₂
        β² = β ⋅ β
        β̂ = β²*area < 1e-14 ? SVector(1.0, 0.0) : β/norm(β)
        layer.β[mode] = β
        layer.γ[mode] = γ = mysqrt(β² - ksq)
        if p == TE
            Y = γ / (im * (k0*η₀) * layer.μᵣ)
            tvec = zhatcross(β̂)
        else
          Y = im * (k0/η₀) * layer.ϵᵣ / γ
          tvec = β̂
        end
        layer.Y[mode] = Y
        layer.tvec[mode] = tvec
        layer.c[mode] = mysqrt(1 / (area * Y))
    end
    return nothing
end

function zhatcross(x::T) where T
    T([-x[2], x[1]])
end

"""
    equal_layer_periodicity(l1::Layer, l2::Layer) -> Bool

Determine if `l1` and `l2` have identical unit cells
by examining their reciprocal lattice vectors.
"""
function equal_layer_periodicity(l1::Layer, l2::Layer)
    denom = norm(l1.β₁) + norm(l1.β₂) + norm(l2.β₁) + norm(l2.β₂)
    norm(l1.β₁ - l2.β₁) / denom < 1e-5 && (norm(l1.β₂ - l2.β₂) / denom < 1e-5)
end

"""
copy_βs_pmn!(l1::Layer, l2::Layer)

Copy fields `β₁`, `β₂`, `P`, `M`, `N` from `l2` to `l1`, modifying `l1`.
"""
@inline function copy_βs_pmn!(l1::Layer, l2::Layer)
    l1.β₁ = copy(l2.β₁)
    l1.β₂ = copy(l2.β₂)
    l1.P = copy(l2.P)
    l1.M = copy(l2.M)
    l1.N = copy(l2.N)
    nothing
end

"""
    fill_mmax_pmn!(layer::Layer, MNmax::Int, k0max, dbmin)

Fill the `layer` modal index arrays `M`, `N`, and `P` given the `layer`
material parameters and reciprocal lattice vectors, the wavenumber `k0max`
at the maximum analysis frequency, and the minimum desired attenuation in dB
`dbmin` for the neglected modes. `MNmax` is the maximum ring number in the
mode lattice to consider.
"""
function fill_mmax_pmn!(layer::Layer, MNmax::Int, k0max, dbmin)
    MNmax ≥ 0 || error("MNmax must be nonnegative")
    n = 2 * (2MNmax+1)^2
    layer.M = zeros(Int, n)
    layer.N = zeros(Int, n)
    layer.P = [TE for _ in 1:n]
    mode = 1
    for ring in 0:MNmax, (m,n) in Ring(ring)
        layer.M[mode:mode+1] .= m
        layer.N[mode:mode+1] .= n
        layer.P[mode:mode+1] .= (TE,TM)
        mode += 2
    end
    select_pmn!(layer, k0max, dbmin)
    return nothing
end

"""
    select_pmn!(layer:Layer, k0max, dbmin)

Select the accessible modes from among the candidates provided
in the `layer`. A mode is considered not accessible if it undergoes more
than `dbmin` decibels of attenuation when propagating/attenuating through
the thickness of dielectric substrate `layer` at the highest frequency of
operation, for any possible angle of incidence of the plane-wave that excites
the FSS structure.

The special value `dbmin=-Inf` instructs this function to retain only
propagating modes
"""
function select_pmn!(layer::Layer, k0max, dbmin=-Inf)
    β₁, β₂ = layer.β₁, layer.β₂  # Reciprocal lattice vectors (1/m)
    kmaxsq = k0max^2  * (layer.ϵᵣ * layer.μᵣ)
    if isfinite(dbmin)
        αh_test = dbmin / 8.6859 # Convert from dB to nepers
        keep = map(1:length(layer.P)) do mode
            (layer.M[mode] == layer.N[mode] == 0) && (return true)
            # Obtain min value of ||βₘₙ|| for this mode for any angle of incidence
            β = norm(layer.M[mode]*β₁ + layer.N[mode]*β₂)
            βmin = β ≤ k0max ? 0.0 : β - k0max
            γ = mysqrt(βmin^2 - kmaxsq) # modal complex atten. constant
            αh = real(γ) * layer.width
            return  αh < αh_test
        end
    else
        keep = map(1:length(layer.P)) do mode
            βvec = layer.M[mode]*β₁ + layer.N[mode]*β₂
            β² = βvec ⋅ βvec
            γ = mysqrt(β² - kmaxsq) # modal complex atten. constant
            return  imag(γ) ≥ real(γ)
        end
    end
    layer.M = layer.M[keep]
    layer.N = layer.N[keep]
    layer.P = layer.P[keep]

    # Eliminate possible duplicate mode triples (p,m,n)
    MN = repeat(unique([layer.M  layer.N], dims=1), inner=(2,1))
    layer.M = MN[:,1]
    layer.N = MN[:,2]
    layer.P = repeat([TE,TM], outer=length(layer.M)÷2)
    return nothing
end

end # module


text:
using ExplainableAI

# NOTE: Heatmapping assumes Flux's WHCN convention (width, height, color channels, batch size).
shape = (2, 2, 3, 1)
A = reshape(collect(Float32, 1:prod(shape)), shape)

shape = (2, 2, 3, 2)
batch = reshape(collect(Float32, 1:prod(shape)), shape)

reducers = [:sum, :maxabs, :norm]
rangescales = [:extrema, :centered]
for r in reducers
    for n in rangescales
        h = heatmap(A; reduce=r, rangescale=n)
        @test_reference "references/heatmaps/reduce_$(r)_rangescale_$(n).txt" h
        @test h ≈ heatmap(A; reduce=r, rangescale=n, unpack_singleton=false)[1]

        h = heatmap(batch; reduce=r, rangescale=n)
        @test_reference "references/heatmaps/reduce_$(r)_rangescale_$(n).txt" h[1]
        @test_reference "references/heatmaps/reduce_$(r)_rangescale_$(n)2.txt" h[2]
    end
end

@test_throws ArgumentError heatmap(A, reduce=:foo)
@test_throws ErrorException heatmap(A, rangescale=:bar)

B = reshape(A, 2, 2, 3, 1, 1)
@test_throws DomainError heatmap(B)
B = reshape(A, 2, 2, 3)
@test_throws DomainError heatmap(B)

A1 = rand(3, 3, 1)
A2 = ExplainableAI._reduce(A1, :sum)
@test A1 == A2


text:
@enum TenorUnits TDays TWeeks TMonths TYears

"""
Structure for tenors, consisting of a multiplier and a unit (e.g. days, month, year)
"""
struct Tenor
    unit:: TenorUnits
    multiplier:: Int
end

const TenorUnitMapping = Dict("D" => TDays, "W" => TWeeks, "M" => TMonths, "Y"=> TYears)
const UnitTenorMapping = Dict(val => key for (key, val) in TenorUnitMapping )

"""
    Tenor(x:: AbstractString)

Constructor for creating a Tenor object from a string.
The input string must start with the multiplier as Integer followed by the
unit (as last character).

Example:

    julia> Tenor.(("1D", "3W", "1M", "10y"))
    (Tenor(TDays, 1), Tenor(TWeeks, 3), Tenor(TMonths, 1), Tenor(TYears, 10))

"""
function Tenor(x:: AbstractString)
    try
        unit = TenorUnitMapping[uppercase(x[end:end])]
    catch ex
        if ex isa KeyError
            @error "$(uppercase(x[end:end])) is not a valid tenor multiplier (occurred in tenor $x)"
        end
        throw(ex)
    end

    multiplier = parse(Int, x[1:end-1])
    # remove ambiguities
    if unit == TDays && multiplier % 7 == 0
        unit = TWeeks
        multiplier = multiplier ÷ 7
    elseif unit == TMonths && multiplier % 12 == 0
        unit = TYears
        multiplier = multiplier ÷ 12
    end
    Tenor(unit, multiplier)
end

"""
String macro for the creation of Tenor objects.

Example:

    julia> t"1W"
    Tenor(Curves.TWeeks, 1)

"""
macro t_str(s)
    esc(:(Tenor($s)))
end

const TenorInDays = Dict(TDays => 1, TWeeks => 7, TMonths => 30, TYears => 365)
# sorting of a Dict is not deterministic
const TenorSorting = sortperm(collect(keys(Curves.TenorInDays)), rev=true)
const Tenors = collect(keys(TenorInDays))[TenorSorting]
const TenorDays = collect(values(TenorInDays))[TenorSorting]

"""
    get_days(x:: Tenor):: Int

Gets the number of days corresponding to the Tenor object.
The days are calculated in a simplified way, assuming 30 days/ month and
365 days/ year.
"""
get_days(x:: Tenor):: Int = TenorInDays[x.unit]*x.multiplier

"""
    get_tenor(x:: Integer):: Tenor

Converts the given number of days to a Tenor.

The following mapping is used (exactly the reverse of `get_days`):
TDays => 1, TWeeks => 7, TMonths => 30, TYears => 365
"""
function get_tenor(x:: Integer):: Tenor
    x > 0 || error("number of days must be positive, obtained $x")
    for (i, unit_days) in enumerate(TenorDays)
        if x % unit_days == 0
            return Tenor(Tenors[i], x ÷ unit_days)
        end
    end
    error("could not convert days $x to tenor")
end

Base.Broadcast.broadcastable(q:: Tenor) = Ref(q) # treat it as a scalar in broadcasting

Base.show(io:: IO, x:: Tenor) = print(io, x.multiplier, UnitTenorMapping[x.unit])

Base.isless(x:: Tenor, y:: Tenor) = isless(get_days(x), get_days(y))

Base.isless(x:: Tenor, y:: AbstractString) = isless(x, Tenor(y))
Base.isless(x:: Tenor, y:: Integer) = isless(get_days(x), y)
Base.isless(x:: AbstractString, y:: Tenor) = isless(Tenor(x), y)
Base.isless(x:: Integer, y:: Tenor) = isless(x, get_days(y))

import Base: ==
==(x:: Tenor, y:: Tenor) = x.unit == y.unit && x.multiplier == y.multiplier
==(x:: Tenor, y:: AbstractString) = ==(x, Tenor(y))
==(x:: Tenor, y:: Integer) = ==(get_days(x), y)
==(x, y:: Tenor) = ==(y, x)


text:
using DataTypesBasic

@test map(string, Const(4)) isa Const{Int}
@test map(string, Const(4)).value == 4
@test eltype(Const(4)) == Any


text:
import HypothesisTests

function model_parameters(model::HypothesisTests.OneSampleTTest)

    parameters = Dict("testname" => HypothesisTests.testname(model),
                "n" => model.n,
                "mean" => model.xbar,
                "stderr" => model.stderr,
                "df" => model.df,
                "u" => model.μ0,
                "statistic" => model.t,
                "CI" => HypothesisTests.confint(model),
                "p" => HypothesisTests.pvalue(model))
    return parameters
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:4] <= 1)

                     @variable(m, -11.08168115564476 <= q <= 12.629895565298366)

                     add_NL_constraint(m, :(log(1 + exp(-0.8013086439641905 + -0.577380120181807*log(1 + exp(-0.15161094843442724 + 0.3901524114721284*$(x[1]) + -0.7620628362671806*$(x[2]) + 0.5531251925343885*$(x[3]) + 0.8877751701397409*$(x[4]))) + -0.9073468716494122*log(1 + exp(-0.9517329762307463 + 0.7324043290359472*$(x[1]) + -0.6011940599882695*$(x[2]) + 0.8132214951552479*$(x[3]) + 0.2016959253369457*$(x[4]))) + -0.34890000179172853*log(1 + exp(0.5700564024677015 + -0.6508434167160906*$(x[1]) + 0.3531203702734036*$(x[2]) + 0.9592602519851474*$(x[3]) + -0.2457804299770685*$(x[4]))))) + log(1 + exp(-0.6400914085711236 + 0.5999574270201746*log(1 + exp(-0.15161094843442724 + 0.3901524114721284*$(x[1]) + -0.7620628362671806*$(x[2]) + 0.5531251925343885*$(x[3]) + 0.8877751701397409*$(x[4]))) + -0.3302091884595031*log(1 + exp(-0.9517329762307463 + 0.7324043290359472*$(x[1]) + -0.6011940599882695*$(x[2]) + 0.8132214951552479*$(x[3]) + 0.2016959253369457*$(x[4]))) + 0.9262651374928534*log(1 + exp(0.5700564024677015 + -0.6508434167160906*$(x[1]) + 0.3531203702734036*$(x[2]) + 0.9592602519851474*$(x[3]) + -0.2457804299770685*$(x[4]))))) + log(1 + exp(0.7171327203675486 + 0.9873737246227612*log(1 + exp(-0.15161094843442724 + 0.3901524114721284*$(x[1]) + -0.7620628362671806*$(x[2]) + 0.5531251925343885*$(x[3]) + 0.8877751701397409*$(x[4]))) + -0.03055409136192777*log(1 + exp(-0.9517329762307463 + 0.7324043290359472*$(x[1]) + -0.6011940599882695*$(x[2]) + 0.8132214951552479*$(x[3]) + 0.2016959253369457*$(x[4]))) + 0.20254351925716962*log(1 + exp(0.5700564024677015 + -0.6508434167160906*$(x[1]) + 0.3531203702734036*$(x[2]) + 0.9592602519851474*$(x[3]) + -0.2457804299770685*$(x[4]))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
abstract UniverseEstimate

type MuSigmaUniverse{T<:MuSigmaModel} <: UniverseEstimate
    Universe::T
    Date::Date
    Data::Timematr
    # EstDetails::Estimator

    function MuSigmaUniverse(univ::T, dat::Date, data::Timematr)
        return new(univ, dat, data)
    end
end

function MuSigmaUniverse{T}(univ::T, dat::Date, data::Timematr)
    return MuSigmaUniverse{T}(univ, dat, data)
end

## check for undefined universe
##-----------------------------

function isDef(univ::UniverseEstimate)
    return isDef(univ.Universe)
end

## display
##--------

import Base.Multimedia.display
function display(mod::UniverseEstimate)
    ## display type information
    typ = typeof(mod)
    println("\ntype: $typ")

    if !isDef(mod)
        println("Model is still undefined!")
        println("   mu: NA")
        println("sigma: NA")
    else
        ## display model
        print("model: ")
        display(typeof(mod.Universe))

        ## display mus
        println("Date: $(mod.Date)")

    end

    ## display data information
    print("Data: ")
    fromDate = mod.Data.idx[1];
    toDate = mod.Data.idx[end];
    nObs, nAss = size(mod.Data)
    println("From: $fromDate, To: $toDate, Size: $nObs x $nAss")
end

##############
## estimate ##
##############

macro defineEstimateForMuSigmaModel(myType)
    esc(quote
        function estimate(t::Type{$myType}, data::Timematr, dat::Date)
            ## estimate model and create UniverseEstimate instance
            modFit = fitModel(t, data, dat)
            return MuSigmaUniverse(modFit, dat, data)
        end
    end)
end

for t = (:(SampleMoments), :(MovWinSampleMoments), :(ExpWeighted))
    eval(macroexpand(:(@defineEstimateForMuSigmaModel($t))))
end


text:
# Copyright (c) 2022 Code Komali
#
# This software is released under the MIT License.
# https://opensource.org/licenses/MIT

module Parameters

# SIGNALS RELATED
const SIGNAL_POS_X_MARGIN = 80.0
const SIGNAL_POS_Y_MARGIN = 80.0
const SIGNAL_RED_TIME = 1500
const SIGNAL_GREEN_TIME = 1200
const SIGNAL_YELLOW_TIME = 300
const SIGNAL_MS = 5

# LANE DIMENSIONS
const H_NUM_LANES = 2
const H_LANE_WIDTH = 50
const H_MEDIAN_WIDTH = 2
const V_NUM_LANES = 2
const V_LANE_WIDTH = 50
const V_MEDIAN_WIDTH = 2

# PEDESTRIAN WALKWAYS AND CROSSWALKS RELATED
const PEDESTRIAN_WALKWAY_WIDTH = 28
# Average CrossPath width 2M i.e 2 * 15.5 = 31
const CROSS_PATH_WIDTH = 31
# Instead of right at the edge of the road (near intersection)
# this margin dictates how much away the cross path is located
const CROSS_PATH_MARGIN = 10

# VEHICLE RELATED
# cars 10 km/h ≈ 2.5 m/s
# lets use milliseconds instead of seconds
# for smooth animation we use 24 frames/sec ≈ 25 f/s
# 25 frames take 1 sec (1000 msec) so 1 frame takes 1000/25 = 40 ms
# ∴ 1 unit of time in our simulation = 40 milliseconds
# if a car moves 2.5 m/1000ms then it moves (2.5/1000) * 40 = 0.1m in 40 ms
const VEHICLE_INITIAL_SPEED = 0.1
# Average Car width 1.8 m
# ∴ Car width = 1.8 * 15.5 ≈ 27.9 units (rounded to 28)
# Average Car length 4.5 m
# ∴ Car length = 4.5 * 15.5 ≈ 68.75 units (rounded to 69)
# source : https://measuringstuff.com/car-length-and-width-measured-in-feet/
const VEHICLE_LENGTH = 69
const VEHICLE_WIDTH = 28

# VEHICLE SPAWN RATE (ONCE EVERY 'N' TICKS)
const VSR_MIN = 100
const VSR_INC = 100
const VSR_MAX = 5000

# ENVIRONMENT RELATED
const EXTENT_WIDTH = 4100
const EXTENT_HEIGHT = 4100

# PLOT AND AESTHETICS RELATED
# road
const ROAD_BOUNDARY_COLOR = :red
const ROAD_BOUNDARY_LW = 2
# lane
const LANE_MARKER_COLOR = :blue
const LANE_MARKER_LS = :dash
# median
const MEDIAN_COLOR = :green
const MEDIAN_LW = 2
# cross path
const CP_COLOR = :darkorange
const CP_LW = 2
# pedestrian walk way boundary (aka platform)
const PW_COLOR = :black
# plot title and alignment
const PLOT_TITLE = "Indian Traffic Simulator"
const PLOT_TITLE_ALIGN = :center

end


text:
module GPUifyLoops

using StaticArrays
using Requires

export @setup, @loop, @synchronize
export @scratch, @shmem

@init @require CUDAnative="be33ccc6-a3ff-5ff2-a52e-74243cff1e17" begin
    using .CUDAnative
end

###
# Simple macros that help to write functions that run
# both on the CPU and GPU
###

iscpu(::Val{:GPU}) = false
iscpu(::Val{:CPU}) = true
sync(::Val{:CPU}) = nothing
sync(::Val{:GPU}) = sync_threads()

macro setup(sym)
    esc(:(local __DEVICE = Val($sym)))
end

macro synchronize()
    esc(:($sync(__DEVICE)))
end

# TODO:
# - check if __DEVICE is defined

macro loop(expr)
    if expr.head != :for
        error("Syntax error: @loop needs a for loop")
    end

    induction = expr.args[1]
    body = expr.args[2]

    if induction.head != :(=)
        error("Syntax error: @loop needs a induction variable")
    end

    rhs = induction.args[2]
    if rhs.head == :block
        @assert length(rhs.args) == 3
        # rhs[2] is a linenode
        cpuidx = rhs.args[1]
        gpuidx = rhs.args[3]

        rhs = Expr(:if, :($iscpu(__DEVICE)), cpuidx, gpuidx)
        induction.args[2] = rhs

        # use cpuidx calculation to check bounds of on GPU.
        bounds_chk = quote
            if !$iscpu(__DEVICE) && !($gpuidx in $cpuidx)
                continue
            end
        end

        pushfirst!(body.args, bounds_chk)
    end

    return esc(Expr(:for, induction, body))
end

include("scratch.jl")
include("shmem.jl")

end


text:
include("envelope.jl")
include("distortion.jl")
include("measure.jl")
include("synth.jl")
include("Instrument.jl")
using WAV
function main()
    getRand(.5)
    env = Envelope(41000,.01,.03,.2,.02,.02)
    a = applyEnvelope((getSine(1,"a4").+getSine(1,"a5").+getSine(1,"e5"))./3,env)
    a2= applyEnvelope((getSine(1,"f#5").+getSine(1,"c4").+getSine(1,"e4"))./3,env)
    b=Float64[]
    c=Float64[]
    #for i in 100.0:-inv(2):1.1 append!(b,sortDistortion(softClip(a,i),Int(ceil(i)))) end
    #for i in π/12:π/12:8π append!(b,softClip(a,abs(8sin(i))+1)) end
    for i in 32π:-π/6:π append!(b,linearDistortion(softClip(a,abs(4sin(i))+1),Int(ceil(i)))) end
    for i in 32π:-π/6:π append!(c,linearDistortion(softClip(a2,abs(4sin(i))+1),Int(ceil(i)))) end
    print(getFreq("a9"))
    print(length(b), " c: ", length(c))
    d = chebby((b.+ c)./2,.5)
    wavplay(d,41000)
 end
 function otherMain()
     env = Envelope(41000,.1,.03,.01,.012,.012)
     b = Float64[]
     c = getSine(1,880)
     e = getSine(1,880.508)
     g = getSine(1,879.982)
     a=(c.+e.+g)./3
     @time a=applyEnvelope(a,Envelope(41000,.07,.03,.7,.07,.03))
     for i in 1:1:8 append!(b,a) end
     for i in 1:inv(10):4 append!(b,chebby(a,i)) end
     for i in 1:1:40 append!(b,softClip(a,i)) end

     @time for i in 1:1:8 append!(b, cheb(a)) end

     wavplay(b,41000)

 end
 function sortTest()
    getRand(.5)
    env = Envelope(41000,.02,.02,.2,.02,.03)
    a = applyEnvelope((getSine(1,"a4").+getSine(1,"a5").+getSine(1,"e5"))./3,env)
    a2= applyEnvelope(getRand(1),env)
    b=Float64[]
    c=Float64[]

    for i in 1:1:200 append!(b,sortDistortion(a,Int(ceil(i)))) end
    for i in 1:1:200 append!(c,sortDistortion(a2,Int(ceil(i)))) end
    print(getFreq("a9"))
    print(length(b), " c: ", length(c))
    #d = chebby((b.+ c)./2,.5)
    wavplay(b,41000)
    wavwrite(b,41000,"a4.wav")
 end

 function disco()
     drumVelope(signal) = applyEnvelope(signal, Envelope(41000,.04,.03,.7,.03,.03))
     bassDrum(x) = softClip(cheb(sortDistortion(drumVelope(getSine(.5,x)),100)), 10)
     cymbal() = drumVelope(getRand(.5))./5
     cym() = append!(cymbal(),cymbal())
     sawVelope(signal) = applyEnvelope(signal, Envelope(41000,.2,.1,.7,.1,.1))
     pad(freq) = sawVelope(getSaw(1, freq))./5
     measure1 = Measure(110,41000.0,4,[])
     addToMeasure!(measure1, cym(),.5)
     addToMeasure!(measure1, cym(),1.5)
     addToMeasure!(measure1, cym(),2.5)
     addToMeasure!(measure1, cymbal(),3.5)

     addToMeasure!(measure1, bassDrum("a2"),0)
     addToMeasure!(measure1, bassDrum("a1"),1)
     addToMeasure!(measure1, bassDrum("a2"),2)
     addToMeasure!(measure1, bassDrum("a3"),3)
     rendered1 = renderMeasure(measure1)
     append!(rendered1,linearDistortion(rendered1,10))
     for i in 0:inv(4):1 addToMeasure!(measure1, pad("a4"),i) end
     for i in 1:inv(4):2 addToMeasure!(measure1, pad("c4"),i) end
     for i in 2:inv(4):3 addToMeasure!(measure1, pad("e4"),i) end
     addToMeasure!(measure1, pad("a5"),3)
     rendered2 = linearDistortion(renderMeasure(measure1),10)
     append!(rendered1,linearDistortion(rendered2,15))
     append!(rendered1,softClip(rendered2,100))
     measure2 = Measure(110,41000.0,4,[])
     addToMeasure!(measure2, drumVelope(pad("a5")),0)
     addToMeasure!(measure2, drumVelope(pad("b5")),.5)
     addToMeasure!(measure2, drumVelope(pad("c5")),1)
     addToMeasure!(measure2, drumVelope(pad("d5")),1.5)
     addToMeasure!(measure2, drumVelope(pad("e5")),2)
     addToMeasure!(measure2, drumVelope(pad("f5")),2.5)
     addToMeasure!(measure2, drumVelope(pad("g5")),3)
     addToMeasure!(measure2, drumVelope(pad("a6")),3.5)
     for i in .5:1:3 addToMeasure!(measure2, linearDistortion(cym(),14),i) end
     addToMeasure!(measure2, bassDrum("e3"),0)
     addToMeasure!(measure2, bassDrum("a3"),1)
     addToMeasure!(measure2, bassDrum("e3"),2)
     addToMeasure!(measure2, bassDrum("a3"),3)
     rendered3 = renderMeasure(measure2)
     append!(rendered1, rendered3)
     addToMeasure!(measure2, cheb(getSine(2,"a2")),0)
     addToMeasure!(measure2, cheb(getSine(2,"c2")),0)
     addToMeasure!(measure2, cheb(getSine(2,"e2")),0)
     addToMeasure!(measure2, cheb(getSine(2,"a3")),0)
     addToMeasure!(measure2, cheb(getSine(2,"c3")),0)
     append!(rendered1, renderMeasure(measure2))
     append!(rendered1, softClip(rendered1,32),30)
     wavplay(rendered1,41000)

 end
 function wavReadTest()

    files = wavread.(["loudSinJingle.wav", "percussion.wav", "summerSin.wav", "violinThing.wav"])
    filesF64 = [] # the float array vector
    for i in files
        push!(filesF64, i[1][:,2])
    end
    f = filesF64[2]

    f = softClip(f,120.1)
   wavplay(f,41000)

 end
 function bassoon(x)
    a = (getSine(1,(2.001*getFreq(x)),41000.0)+
        getSine(1,(.5003*getFreq(x)),41000.0)+
        getSine(1,(0.998*getFreq(x)),41000.0))/3

        return applyEnvelope(softClip(a,3), Envelope(41000.0,.09,.09,7,.1,.1))

end
function hat()
    drumVelope(signal) = applyEnvelope(signal, Envelope(41000,.03,.02,.5,.03,.03))
    a = Float64[]
    a = drumVelope(getRand(.5))
    append!(a, drumVelope(getRand(.5)))
    return a

end
function sawBASS(x)
    drumVelope(signal) = applyEnvelope(signal, Envelope(41000,.1,.02,.5,.03,.05))
    return cheb(drumVelope(getSaw(1, x)))

end
function noteTest()
    bpm = 130
    m = Measure(bpm, 41000, 8,[])
    addNotes(m, ["a4", "e7", "d6", "d4","e4", "e4" , "a3","a4"],1,bassoon)
    addNotes(m, ["c5", "c4", "f5", "f4","g#5", "g#4" , "c5","c4"],.98,bassoon)
    addNotes(m, ["e5", "e3", "a5", "a6","b6", "d6" , "a6","e5"],.97,bassoon)
    addNotes(m, ["e2", "e2", "a2", "a2","b2", "d2" , "a2","e2"],1,sawBASS)

    a = renderMeasure(m)
    f = Float64[]
    append!(f,a)
   for i in .5:1:7.5 addToMeasure!(m,hat(),i+i*.001) end
   append!(f,renderMeasure(m))
   addNotes(m, ["e2", "e2", "a2", "a2","b2", "d2" , "a2","e2"],1,sawBASS)
   addNotes(m, ["e2", "e2", "a2", "a2","b2", "d2" , "a2","e2"],.9,sawBASS)

   for i in 0:2:7.5 addToMeasure!(m,hat(),i) end
   addNotes(m, ["e2", "e2", "a2", "a2","b2", "d2" , "a2"],1.1,sawBASS)
   addNotes(m, ["e2", "e2", "a2", "a2","b2", "d2" , "a2"],1.2,sawBASS)
   append!(f,renderMeasure(m))

   addNotes(m, ["a4", "e7", "d6", "d4","e4", "b4" , "a3"],1.1,bassoon)
   addNotes(m, ["c5", "c4", "f5", "f4","g#5", "g#4" ],1.2,bassoon)
   addNotes(m, ["e5", "e3", "a5", "a6","b6", "d6" , "a6"],1.05,bassoon)
   addNotes(m, ["e2", "e2", "a2", "a2","b2", "b3" , "a2"],1.12,sawBASS)
   append!(f,renderMeasure(m))
   for i in .5:1:7.5 addToMeasure!(m,hat(),i+i*.003) end
   for i in 0:2:7.5 addToMeasure!(m,hat(),i*.002) end
   addNotes(m, ["a4", "e7", "d6", "d4","e4", "b4" , "a3"],1.1,bassoon)
   addNotes(m, ["c5", "c4", "f5", "f4","g#5", "g#4" ],1.2,bassoon)
   addNotes(m, ["e5", "e3", "a5", "a6","b6", "d6" , "a6"],1.05,bassoon)
   addNotes(m, ["e2", "e2", "a2", "a2","b2", "b3" , "a2"],1.12,sawBASS)
   append!(f,renderMeasure(m))
   addNotes(m, ["a4", "e7", "d6", "d4","e4", "b4" , "a3"],.6,bassoon)
   addNotes(m, ["c5", "c4", "f5", "f4","g#5", "g#4" ],.6,bassoon)
   addNotes(m, ["e5", "e3", "a5", "a6","b6", "d6" , "a6","e5", "e3", "a5", "a6","b6"],.6,bassoon)
   addNotes(m, ["e2", "e2", "a2", "a2","b2", "b3" , "a2","a3"],.6,sawBASS)
   append!(f,renderMeasure(m))
    wavplay(f, 41000)
    #wavwrite(f,41000,"sorting.wav")

end
function lfoTest()
    a = (getSine(4,"a4").+getRand(4,))/2
    cheb(x,limit)=limit*x[1]^3+(1-limit)x[1]
    a = lfo(a,.4,cheb,(.1,3.9))
    a = lfo(a,.9,cheb,(0,2))
    a = lfo(a,3,cheb,(3,4))
    wavplay(a,41000)

end
function chebTEST()
    cheb(x,limit)=limit*x[1]^3+(1-limit)x[1]
    bpm = 125
    m = Measure(bpm, 41000, 8,[])
    sawVelope(signal) = applyEnvelope(signal, Envelope(41000,.17,.1,.7,.1,.1))
    pad(freq) = softClip(lfo(sawVelope(getSine(1, freq))./5,.2,cheb,(2,3)),50)

    addNotes(m, ["b2", "b3", "d#3", "c3","c#3", "c#3" , "e3","e2"],1,pad)
    addNotes(m, ["f#4", "f#4", "d#4", "f#4","g#4", "g#4" , "b5","b5"],1,bassoon)
    a = renderMeasure(m)
    wavplay(a,41000)
end
function chebbyBrato()
    c(x,order) = cos(order*acos(x))
    a = (getSine(5,"f4")
    .+getSine(5,"a4")
    .+getSine(5,"c5")
    #.+getRand(5)/100
    .+getSine(5,"f2")
    )
    a = a/(maximum(a)+1)

    b = lfo(a,.5,c,(0,5))
    wavplay(b,41000)
end
function adaaTest()
    a = getSine(1,"c4")
    b = getSine(10,"f4")
    @time hardC(a,10)
    @time hardClipAD1(a,10)
    @time hardClipAD2(a,10)
    @time hardC(b,40)
    @time hardClipAD1(b,40)
    @time hardClipAD2(b,40)

    f(x) = (abs(x) < 1 ? x : sign(x) )
    f1(x) = (abs(x) < 1 ? x^2/2 : x*sign(x) - .5)
    f2(x) = (abs(x) < 1 ? x^3/6 : (x^2/2 + 1/6) * sign(x) - x/2)

end
function whatever(x::Vector{Float64},c)::Vector{Float64}
    b = c^2
    a = 1-b-c
    f(t) = b*t + a*(4t^3-3t) + c*(16t^5-20t^3+5t)
    return f.(x)
end
function w()
    a = getRand(1)
    @time hardClipAD2(a,1)
    @time whatever(a,.5)
    @time cheb(a)
    @time hardC(a,123)
    println()

    b = getRand(100)
    @time hardClipAD2(b,1)
    @time whatever(b,.5)
    @time cheb(b)
    @time hardC(b,123)
end
function organ(note)
    hz = getFreq(note)
    return applyEnvelope(hardClipAD2(compound(hz),200), Envelope(41000.0,.15,.09,7,.1,.04))
end
function slow_organ(note)
    hz = getFreq(note)

    return applyEnvelope(sortDistortion(softClipAD2(compound2(hz),5),50), Envelope(41000.0,.3,.19,7,.2,.15))
end
function lacrimosa()
    bpm = 140
    m = Measure(bpm, 41000, 12,[])
    r = nothing
    addNotes(m, [r, "c#5", "d5", r,"a5","a#5", r,
    "d5" , "c#5",r,"c6","a#5"],1,organ)
    addNotes(m, ["d4","f4","e4","g4"],3,slow_organ)
    addNotes(m, ["f4","a4","g4","c#4"],3,slow_organ)
    audio = renderMeasure(m)

    m2 = Measure(bpm, 41000, 12,[])
    addNotes(m2, [r, "a5", "d6", r,"a#5","g5", r,
    "e5" , "f5",r,"a5","c#5"],1,organ)
    addNotes(m2, ["d5","e4","d4","g4"],3,slow_organ)
    addNotes(m2, ["f4","g3","a3","a3"],3,slow_organ)
    append!(audio, renderMeasure(m2))
    wavplay(audio,41000)
end
#disco()
lacrimosa()
#noteTest()
#main()
#chebbyBrato()
#wavplay(a,41000)

max_stars_repo_name: justinbroce/julia-audio-stuffs

max_issues_repo_name: justinbroce/julia-audio-stuffs

max_forks_repo_name: justinbroce/julia-audio-stuffs


text:
struct RandomOutlierPQs <: ModelBasedPQs
    occ::SVDD.OCClassifier
end

function qs_score(qs::RandomOutlierPQs, x::Array{T, 2}, labels::Dict{Symbol, Array{Int, 1}})::Array{Float64, 1} where T <: Real
    prediction = SVDD.classify.(SVDD.predict(qs.occ, x))
    scores = rand(size(x, 2))
    scores[prediction .== :inlier] .= 0
    return scores
end


text:
using Base.Test

include("test_acoustic.jl")

max_stars_repo_name: WneleiGao/AcousticWaveLab

max_issues_repo_name: WneleiGao/AcousticWaveLab

max_forks_repo_name: WneleiGao/AcousticWaveLab


text:
function probe(iter::Class_iterate, kkt_solver::abstract_KKT_system_solver, pars::Class_parameters, timer::class_advanced_timer)
    kkt_associate_rhs!(kkt_solver, iter, Reduct_affine(), timer)
    compute_direction!(kkt_solver, timer)

    lb_s = lb_s_predict(iter, kkt_solver.dir, pars)

    #@show minimum(abs(iter.point.s ./ kkt_solver.dir.y)), minimum(abs(iter.point.s ./ iter.point.y))
    #if true
    #  KNITRO_guess(iter,kkt_solver.dir, pars)
    #end

    m = length(iter.point.s)
    alpha_s = simple_max_step(iter.point.s, kkt_solver.dir.s, lb_s)
    alpha_y = simple_max_step(iter.point.y, kkt_solver.dir.y, zeros(m))
    return min(alpha_s, alpha_y)
end

#=function take_step!(iter::Class_iterate, eta::Class_reduction_factors, kkt_solver::abstract_KKT_system_solver, ls_mode::Symbol, filter::Array{Class_filter,1}, pars::Class_parameters, min_step_size::Float64, timer::class_advanced_timer)
    kkt_associate_rhs!(kkt_solver, iter, eta, timer)
    compute_direction!(kkt_solver, timer)

    if false #kkt_solver.kkt_err_norm.ratio > 1e-1
        @show kkt_solver.kkt_err_norm
        return :failure, iter, Class_ls_info()
    end

    success, iter, step_info = simple_ls(iter, kkt_solver.dir, ls_mode, filter, pars, min_step_size, timer)

    return success, iter, step_info
end=#

function take_step2!( be_aggressive::Bool, iter::Class_iterate, kkt_solver::abstract_KKT_system_solver, filter::Array{Class_filter,1}, pars::Class_parameters, timer::class_advanced_timer)
    if be_aggressive
        if pars.ls.agg_gamma == :mehrotra || pars.ls.agg_gamma == :mehrotra_stb
          σ = probe(iter, kkt_solver, pars, timer)
          #gamma = (1.0 - σ)
          gamma = min(0.5, (1.0 - σ)^2)
          #gamma = (1.0 - σ) #min(0.5, )
          #gamma = min(0.5,(1.0 - σ)^1.5)
          if pars.ls.agg_gamma == :mehrotra
            reduct_factors = Class_reduction_factors(gamma, gamma, gamma)
          elseif pars.ls.agg_gamma == :mehrotra_stb
            reduct_factors = Class_reduction_factors(gamma, 0.0, gamma)
          end
        else
          if pars.ls.agg_gamma == :affine
            reduct_factors = Reduct_affine()
          elseif pars.ls.agg_gamma == :constant
            reduct_factors = Class_reduction_factors(0.2, 0.0, 0.2)
          else
            error("pars.ls.agg_gamma choice not recognized")
          end
        end

        ls_mode = :accept_aggressive
        min_step_size = pars.ls.min_step_size_agg_ratio * min(1.0, 1.0 / maximum(- get_primal_res(iter) ./ iter.point.s))
    else
        ls_mode = :accept_filter
        reduct_factors = Reduct_stable()
        min_step_size = pars.ls.min_step_size_stable
    end

    kkt_associate_rhs!(kkt_solver, iter, reduct_factors, timer)
    compute_direction!(kkt_solver, timer)

    if false #kkt_solver.kkt_err_norm.ratio > 1e-1
        @show kkt_solver.kkt_err_norm
        return :failure, iter, Class_ls_info(), reduct_factors
    end

    success, iter, step_info = simple_ls(iter, kkt_solver.factor_it, kkt_solver.dir, ls_mode, filter, pars, min_step_size, timer)
    return success, iter, step_info, reduct_factors
end


text:
println("Hello world again!")

a = 1
b = [a,2]

max_stars_repo_name: sylvaticus/IntroSPMLJuliaCourse

max_issues_repo_name: sylvaticus/IntroSPMLJuliaCourse

max_forks_repo_name: sylvaticus/IntroSPMLJuliaCourse


text:
"""
    WienerKernel{i}()

i-times integrated Wiener process kernel function.

- For i=-1, this is just the white noise covariance, see [`WhiteKernel`](@ref).
- For i= 0, this is the Wiener process covariance,
- For i= 1, this is the integrated Wiener process covariance (velocity),
- For i= 2, this is the twice-integrated Wiener process covariance (accel.),
- For i= 3, this is the thrice-integrated Wiener process covariance,

where ``κᵢ`` is given by

```math
    κ₋₁(x, y) =  δ(x, y)
    κ₀(x, y)  =  min(x, y)
```
and for ``i >= 1``,
```math
    κᵢ(x, y) = 1 / aᵢ * min(x, y)^(2i + 1) + bᵢ * min(x, y)^(i + 1) * |x - y| * rᵢ(x, y),
```
    with the coefficients ``aᵢ``, ``bᵢ`` and the residual ``rᵢ(x, y)`` defined as follows:
```math
    a₁ = 3, b₁ = 1/2, r₁(x, y) = 1,
    a₂ = 20, b₂ = 1/12, r₂(x, y) = x + y - min(x, y) / 2,
    a₃ = 252, b₃ = 1/720, r₃(x, y) = 5 * max(x, y)² + 2 * x * z + 3 * min(x, y)²

```

# References:
See the paper *Probabilistic ODE Solvers with Runge-Kutta Means* by Schober, Duvenaud and
Hennig, NIPS, 2014, for more details.

"""
struct WienerKernel{I} <: Kernel
    function WienerKernel{I}() where I
        @assert I ∈ (-1, 0, 1, 2, 3) "Invalid parameter i=$(I). Should be -1, 0, 1, 2 or 3."
        if I == -1
            return WhiteKernel()
        end
        return new{I}()
    end
end

function WienerKernel(;i::Integer=0)
    return WienerKernel{i}()
end

function (::WienerKernel{0})(x, y)
    X = sqrt(sum(abs2, x))
    Y = sqrt(sum(abs2, y))
    return min(X, Y)
end

function (::WienerKernel{1})(x, y)
    X = sqrt(sum(abs2, x))
    Y = sqrt(sum(abs2, y))
    minXY = min(X, Y)
    return 1 / 3 * minXY^3 + 1 / 2 * minXY^2 * euclidean(x, y)
end

function (::WienerKernel{2})(x, y)
    X = sqrt(sum(abs2, x))
    Y = sqrt(sum(abs2, y))
    minXY = min(X, Y)
    return 1 / 20 * minXY^5 + 1 / 12 * minXY^3 * euclidean(x, y) *
        ( X + Y - 1 / 2 * minXY )
end

function (::WienerKernel{3})(x, y)
    X = sqrt(sum(abs2, x))
    Y = sqrt(sum(abs2, y))
    minXY = min(X, Y)
    return 1 / 252 * minXY^7 + 1 / 720 * minXY^4 * euclidean(x, y) *
        ( 5 * max(X, Y)^2 + 2 * X * Y + 3 * minXY^2 )
end

Base.show(io::IO, ::WienerKernel{I}) where I = print(io, I, "-times integrated Wiener kernel")


text:
module ConstantCellFieldsTests

using Gridap

using ..FieldsMocks

l = 10
f = MockField(2,Float64)

cf = ConstantCellValue(f,l)

p1 = Point(1,2)
p2 = Point(0,3)
p3 = Point(7,4)
p4 = Point(8,1)
p = [p1,p2,p3,p4]

cp = ConstantCellValue(p,l)

r = evaluate(f,p)
rg = evaluate(∇(f),p)
v = [r for i in 1:l]
g = [rg for i in 1:l]

test_index_cell_field(cf,cp,v,g)

end # module


text:
ipa = Dict()
ipb = Dict()
ops = Dict()

# Input membership functions
ipa["small"] = TriangularMF(1, 2, 3)
ipa["large"] = TriangularMF(5, 6, 7)
ipb["small"] = TriangularMF(1, 2, 3)
ips = [ipa, ipb]

# Output membership functions
ops["small"] = TriangularMF(1, 2, 3)
ops["large"] = TriangularMF(5, 6, 7)

# Simple rule
rule1 = Rule(["small", "small"], "small")
rule2 = Rule(["large", "small"], "large")
rules = [rule1, rule2]

fis = FISMamdani(ips, ops, rules)

# Input values
ival1 = [1.2, 2.3]
ival2 = [5.7, 2.0]

@assert eval_fis(fis, ival1) == 2.0
@assert eval_fis(fis, ival2) == 6.0

service = Dict( "poor"      => TrapezoidalMF(0.0, 0.0, 2.0, 4.0),
                "good"      => TriangularMF(3., 5., 7.),
                "excellent" => TrapezoidalMF(6., 8., 10., 10.))

food = Dict("rancid"    => TrapezoidalMF(0., 0., 3., 6.),
            "delicious" => TrapezoidalMF(4., 7., 10., 10.))

inputs = [service, food]

tip  = Dict("cheap"     => TrapezoidalMF(10., 10., 20., 30.),
            "average"   => TriangularMF(20., 30., 40.),
            "generous"  => TrapezoidalMF(30., 40, 50., 50.))

rule1 = Rule(["poor", "rancid"], "cheap", "MAX")

# How to not include food input ?
rule2 = Rule(["good", ""], "average", "MAX")

rule3 = Rule(["excellent", "delicious"], "generous", "MAX")
rules = [rule1, rule2, rule3]

fis_tips = FISMamdani(inputs, tip, rules)

in_vals = [7., 8,]

@assert eval_fis(fis_tips, in_vals) == 45.0


text:
function Rect(geometry::AbstractArray{<:Point{N,T}}) where {N,T}
    return Rect{N,T}(geometry)
end

"""
Construct a HyperRectangle enclosing all points.
"""
function Rect{N1,T1}(geometry::AbstractArray{PT}) where {N1,T1,PT<:AbstractPoint}
    N2, T2 = length(PT), eltype(PT)
    @assert N1 >= N2
    vmin = Point{N2,T2}(typemax(T2))
    vmax = Point{N2,T2}(typemin(T2))
    for p in geometry
        vmin, vmax = minmax(p, vmin, vmax)
    end
    o = vmin
    w = vmax - vmin
    return if N1 > N2
        z = zero(Vec{N1 - N2,T1})
        Rect{N1,T1}(vcat(o, z), vcat(w, z))
    else
        Rect{N1,T1}(o, w)
    end
end

function Rect(primitive::GeometryPrimitive{N,T}) where {N,T}
    return Rect{N,T}(primitive)
end

function Rect{T}(primitive::GeometryPrimitive{N,T}) where {N,T}
    return Rect{N,T}(primitive)
end

function Rect{T}(a::Pyramid) where {T}
    w, h = a.width / T(2), a.length
    m = Vec{3,T}(a.middle)
    return Rect{T}(m .- Vec{3,T}(w, w, 0), m .+ Vec{3,T}(w, w, h))
end

function Rect{T}(a::Sphere) where {T}
    mini, maxi = extrema(a)
    return Rect{T}(mini, maxi .- mini)
end

Rect{T}(a) where {T} = Rect{T}(coordinates(a))
Rect{N,T}(a) where {N,T} = Rect{N,T}(coordinates(a))


text:
import MLJBase: package_name, package_url, package_license,
    is_wrapper, is_pure_julia, input_scitype, target_scitype,
    supports_weights, docstring, load_path

@testset "meta-pkg" begin
    @show methods(package_name)
    rr = RidgeRegressor()
    @test package_name(rr) == "MLJLinearModels"
    @test package_url(rr) ==
        "https://github.com/alan-turing-institute/MLJLinearModels.jl"
    @test package_license(rr) == "MIT"
    @test is_pure_julia(rr) == true
    @test is_wrapper(rr) == false
end

@testset "meta-reg" begin
    lr = LinearRegressor()
    @test input_scitype(lr) == MLJBase.Table(MLJBase.Continuous)
    @test target_scitype(lr) == AbstractVector{MLJBase.Continuous}
    @test supports_weights(lr) == false
    @test !isempty(docstring(lr))
    @test load_path(lr) == "MLJLinearModels.LinearRegressor"
end

@testset "meta-clf" begin
    lr = LogisticClassifier()
    @test input_scitype(lr) == MLJBase.Table(MLJBase.Continuous)
    @test target_scitype(lr) == AbstractVector{<:MLJBase.Finite}
    @test supports_weights(lr) == false
    @test !isempty(docstring(lr))
    @test load_path(lr) == "MLJLinearModels.LogisticClassifier"
end


text:
#
# mappings.jl -
#
# Provide basic mappings.
#
#-------------------------------------------------------------------------------
#
# This file is part of LazyAlgebra (https://github.com/emmt/LazyAlgebra.jl)
# released under the MIT "Expat" license.
#
# Copyright (c) 2017-2020 Éric Thiébaut.
#

#------------------------------------------------------------------------------
# IDENTITY AND UNIFORM SCALING

identical(::Identity, ::Identity) = true

@callable Identity

# Traits:
SelfAdjointType(::Identity) = SelfAdjoint()
MorphismType(::Identity) = Endomorphism()
DiagonalType(::Identity) = DiagonalMapping()

apply(::Type{<:Operations}, ::Identity, x, scratch::Bool=false) = x

# vcreate for identity always return x (see doc. of vcreate).
vcreate(::Type{<:Operations}, ::Identity, x, scratch::Bool) = x

apply!(α::Number, ::Type{<:Operations}, ::Identity, x, ::Bool, β::Number, y) =
    vcombine!(y, α, x, β, y)

# Rules to automatically convert UniformScaling from standard library module
# LinearAlgebra into λ*Id.  For other operators, there is no needs to extend ⋅
# (\cdot) and ∘ (\circ) as they are already converted in calls to *.  But in
# the case of UniformScaling, we must explicitly do that for * and for ∘ (not
# for ⋅ which is replaced by a * by existing rules).
for op in (:(+), :(-), :(*), :(∘), :(/), Symbol("\\"))
    @eval begin
        Base.$op(A::UniformScaling, B::Mapping) = $op(Mapping(A), B)
        Base.$op(A::Mapping, B::UniformScaling) = $op(A, Mapping(B))
    end
end

#------------------------------------------------------------------------------
# SYMBOLIC MAPPINGS (FOR TESTS)

struct SymbolicMapping{T} <: Mapping end
struct SymbolicLinearMapping{T} <: LinearMapping end
SymbolicMapping(id::AbstractString) = SymbolicMapping(Symbol(id))
SymbolicMapping(id::Symbol) = SymbolicMapping{Val{id}}()
SymbolicLinearMapping(id::AbstractString) = SymbolicLinearMapping(Symbol(id))
SymbolicLinearMapping(x::Symbol) = SymbolicLinearMapping{Val{x}}()

show(io::IO, A::SymbolicMapping{Val{T}}) where {T} = print(io, T)
show(io::IO, A::SymbolicLinearMapping{Val{T}}) where {T} = print(io, T)

identical(::T, ::T) where {T<:SymbolicMapping} = true
identical(::T, ::T) where {T<:SymbolicLinearMapping} = true

#------------------------------------------------------------------------------
# NON-UNIFORM SCALING

"""
    Diag(A) -> NonuniformScaling(A)

yields a non-uniform scaling linear mapping (of type `NonuniformScaling`) whose
effect is to apply elementwise multiplication of its argument by the scaling
factors `A`.  This mapping can be thought as a *diagonal* operator.

The `diag` method in `LinearAlgebra` can be called to retrieve the scaling
factors:

    using LinearAlgebra
    W = Diag(A)
    diag(W) === A  # this is true

!!! note
    Beware of the differences between the [`Diag`](@ref) (with an uppercase
    'D') and [`diag`](@ref) (with an lowercase 'd') methods.

"""
struct NonuniformScaling{T} <: LinearMapping
    diag::T
end

const Diag{T} = NonuniformScaling{T}

@callable NonuniformScaling

# Traits:
MorphismType(::NonuniformScaling) = Endomorphism()
DiagonalType(::NonuniformScaling) = DiagonalMapping()
SelfAdjointType(A::NonuniformScaling) =
    _selfadjointtype(eltype(coefficients(A)), A)
_selfadjointtype(::Type{<:Real}, ::NonuniformScaling) =
    SelfAdjoint()
_selfadjointtype(::Type{<:Complex}, ::NonuniformScaling) =
    NonSelfAdjoint()

coefficients(A::NonuniformScaling) = A.diag
LinearAlgebra.diag(A::NonuniformScaling) = coefficients(A)

identical(A::T, B::T) where {T<:NonuniformScaling} =
    coefficients(A) === coefficients(B)

function inv(A::NonuniformScaling{<:AbstractArray{T,N}}
             ) where {T<:AbstractFloat, N}
    q = coefficients(A)
    r = similar(q)
    @inbounds @simd for i in eachindex(q, r)
        r[i] = one(T)/q[i]
    end
    return NonuniformScaling(r)
end

eltype(::Type{<:NonuniformScaling{<:AbstractArray{T,N}}}) where {T, N} = T

input_ndims(::NonuniformScaling{<:AbstractArray{T,N}}) where {T, N} = N
input_size(A::NonuniformScaling{<:AbstractArray}) = size(coefficients(A))
input_size(A::NonuniformScaling{<:AbstractArray}, i) =
    size(coefficients(A), i)

output_ndims(::NonuniformScaling{<:AbstractArray{T,N}}) where {T, N} = N
output_size(A::NonuniformScaling{<:AbstractArray}) = size(coefficients(A))
output_size(A::NonuniformScaling{<:AbstractArray}, i) =
    size(coefficients(A), i)

# Simplify left multiplication (and division) by a scalar.
# FIXME: α = 0 should be treated specifically
*(α::Number, A::NonuniformScaling)::NonuniformScaling =
    (α == 1 ? A : NonuniformScaling(vscale(α, coefficients(A))))

# Extend composition of diagonal operators.
*(A::NonuniformScaling, B::NonuniformScaling) =
    NonuniformScaling(vproduct(coefficients(A), coefficients(B)))

function apply!(α::Number,
                ::Type{P},
                W::NonuniformScaling{<:AbstractArray{Tw,N}},
                x::AbstractArray{Tx,N},
                scratch::Bool,
                β::Number,
                y::AbstractArray{Ty,N}) where {P<:Operations,
                                               Tw<:Floats,
                                               Tx<:Floats,
                                               Ty<:Floats,N}
    w = coefficients(W)
    I = all_indices(w, x, y)
    if α == 0
        vscale!(y, β)
    elseif β == 0
        if α == 1
            _apply_diagonal!(P, axpby_yields_x, I, 1, w, x, 0, y)
        else
            a = promote_multiplier(α, Tw, Tx)
            _apply_diagonal!(P, axpby_yields_ax, I, a, w, x, 0, y)
        end
    elseif β == 1
        if α == 1
            _apply_diagonal!(P, axpby_yields_xpy, I, 1, w, x, 1, y)
        else
            a = promote_multiplier(α, Tw, Tx)
            _apply_diagonal!(P, axpby_yields_axpy, I, a, w, x, 1, y)
        end
    else
        b = promote_multiplier(β, Ty)
        if α == 1
            _apply_diagonal!(P, axpby_yields_xpby, I, 1, w, x, b, y)
        else
            a = promote_multiplier(α, Tw, Tx)
            _apply_diagonal!(P, axpby_yields_axpby, I, a, w, x, b, y)
        end
    end
    return y
end

function _apply_diagonal!(::Type{Direct}, axpby::Function, I,
                          α, w, x, β, y)
    @inbounds @simd for i in I
        y[i] = axpby(α, w[i]*x[i], β, y[i])
    end
end

function _apply_diagonal!(::Type{Adjoint}, axpby::Function, I,
                          α, w, x, β, y)
    @inbounds @simd for i in I
        y[i] = axpby(α, conj(w[i])*x[i], β, y[i])
    end
end

function _apply_diagonal!(::Type{Inverse}, axpby::Function, I,
                          α, w, x, β, y)
    @inbounds @simd for i in I
        y[i] = axpby(α, x[i]/w[i], β, y[i])
    end
end

function _apply_diagonal!(::Type{InverseAdjoint}, axpby::Function, I,
                          α, w, x, β, y)
    @inbounds @simd for i in I
        y[i] = axpby(α, x[i]/conj(w[i]), β, y[i])
    end
end

function vcreate(::Type{<:Operations},
                 W::NonuniformScaling{<:AbstractArray{Tw,N}},
                 x::AbstractArray{Tx,N},
                 scratch::Bool) where {Tw,Tx,N}
    inds = same_axes(coefficients(W), x)
    T = promote_type(Tw, Tx)
    return (scratch && Tx == T ? x : similar(Array{T}, inds))
end

#------------------------------------------------------------------------------
# RANK-1 OPERATORS

"""
    RankOneOperator(u, v) -> A

yields the rank one linear operator `A = u⋅v'` defined by the two *vectors* `u`
and `v` and behaving as:

    A*x  -> vscale(vdot(v, x)), u)
    A'*x -> vscale(vdot(u, x)), v)

See also: [`SymmetricRankOneOperator`](@ref), [`LinearMapping`](@ref),
          [`apply!`](@ref), [`vcreate`](@ref).

"""
struct RankOneOperator{U,V} <: LinearMapping
    u::U
    v::V
end

@callable RankOneOperator

function apply!(α::Number, ::Type{Direct}, A::RankOneOperator,
                x, scratch::Bool, β::Number, y)
    return _apply_rank_one!(α, A.u, A.v, x, β, y)
end

function apply!(α::Number, ::Type{Adjoint}, A::RankOneOperator,
                x, scratch::Bool, β::Number, y)
    return _apply_rank_one!(α, A.v, A.u, x, β, y)
end

function _apply_rank_one!(α::Number, u, v, x, β::Number, y)
    if α == 0
        # Lazily assume that y has correct type, dimensions, etc.
        vscale!(y, β)
    else
        vcombine!(y, α*vdot(v, x), u, β, y)
    end
    return y
end

# Lazily assume that x has correct type, dimensions, etc.
# FIXME: optimize when scratch=true
vcreate(::Type{Direct}, A::RankOneOperator, x, scratch::Bool) = vcreate(A.v)
vcreate(::Type{Adjoint}, A::RankOneOperator, x, scratch::Bool) = vcreate(A.u)

input_type(A::RankOneOperator{U,V}) where {U,V} = V
input_ndims(A::RankOneOperator) = ndims(A.v)
input_size(A::RankOneOperator) = size(A.v)
input_size(A::RankOneOperator, d...) = size(A.v, d...)
input_eltype(A::RankOneOperator) = eltype(A.v)

output_type(A::RankOneOperator{U,V}) where {U,V} = U
output_ndims(A::RankOneOperator) = ndims(A.u)
output_size(A::RankOneOperator) = size(A.u)
output_size(A::RankOneOperator, d...) = size(A.u, d...)
output_eltype(A::RankOneOperator) = eltype(A.u)

identical(A::T, B::T) where {T<:RankOneOperator} =
    ((A.u === B.u)&(A.v === B.v))

"""
    SymmetricRankOneOperator(u) -> A

yields the symmetric rank one operator `A = u⋅u'` defined by the *vector* `u`
and behaving as follows:

    A'*x -> A*x
    A*x  -> vscale(vdot(u, x)), u)

See also: [`RankOneOperator`](@ref), [`LinearMapping`](@ref),
          [`Trait`](@ref) [`apply!`](@ref), [`vcreate`](@ref).

"""
struct SymmetricRankOneOperator{U} <: LinearMapping
    u::U
end

@callable SymmetricRankOneOperator

# Traits:
MorphismType(::SymmetricRankOneOperator) = Endomorphism()
SelfAdjointType(::SymmetricRankOneOperator) = SelfAdjoint()

function apply!(α::Number, ::Type{<:Union{Direct,Adjoint}},
                A::SymmetricRankOneOperator, x, scratch::Bool, β::Number, y)
    return _apply_rank_one!(α, A.u, A.u, x, β, y)
end

function vcreate(::Type{<:Union{Direct,Adjoint}},
                 A::SymmetricRankOneOperator, x, scratch::Bool)
    # Lazily assume that x has correct type, dimensions, etc.
    return (scratch ? x : vcreate(x))
end

input_type(A::SymmetricRankOneOperator{U}) where {U} = U
input_ndims(A::SymmetricRankOneOperator) = ndims(A.u)
input_size(A::SymmetricRankOneOperator) = size(A.u)
input_size(A::SymmetricRankOneOperator, d...) = size(A.u, d...)
input_eltype(A::SymmetricRankOneOperator) = eltype(A.u)

output_type(A::SymmetricRankOneOperator{U}) where {U} = U
output_ndims(A::SymmetricRankOneOperator) = ndims(A.u)
output_size(A::SymmetricRankOneOperator) = size(A.u)
output_size(A::SymmetricRankOneOperator, d...) = size(A.u, d...)
output_eltype(A::SymmetricRankOneOperator) = eltype(A.u)

identical(A::T, B::T) where {T<:SymmetricRankOneOperator} =
    (A.u === B.u)

#------------------------------------------------------------------------------
# GENERALIZED MATRIX AND MATRIX-VECTOR PRODUCT

"""
    GeneralMatrix(A)

creates a linear mapping whose coefficients are given by a multi-dimensional
array `A` and which generalizes the definition of the matrix-vector product
without calling `reshape` to change the dimensions.

For instance, assuming that `G = GeneralMatrix(A)` with `A` a regular array,
then `y = G*x` requires that the dimensions of `x` match the trailing
dimensions of `A` and yields a result `y` whose dimensions are the remaining
leading dimensions of `A`, such that `axes(A) = (axes(y)..., axes(x)...)`.
Applying the adjoint of `G` as in `y = G'*x` requires that the dimensions of
`x` match the leading dimension of `A` and yields a result `y` whose dimensions
are the remaining trailing dimensions of `A`, such that `axes(A) = (axes(x)...,
axes(y)...)`.

See also: [`reshape`](@ref).

"""
struct GeneralMatrix{T<:AbstractArray} <: LinearMapping
    arr::T
end

@callable GeneralMatrix

coefficients(A) = A.arr

# Make a GeneralMatrix behaves like an ordinary array.
eltype(A::GeneralMatrix) = eltype(coefficients(A))
length(A::GeneralMatrix) = length(coefficients(A))
ndims(A::GeneralMatrix) = ndims(coefficients(A))
axes(A::GeneralMatrix) = axes(coefficients(A))
size(A::GeneralMatrix) = size(coefficients(A))
size(A::GeneralMatrix, i...) = size(coefficients(A), i...)
getindex(A::GeneralMatrix, i...) = getindex(coefficients(A), i...)
setindex!(A::GeneralMatrix, x, i...) = setindex!(coefficients(A), x, i...)
stride(A::GeneralMatrix, k) = stride(coefficients(A), k)
strides(A::GeneralMatrix) = strides(coefficients(A))
eachindex(A::GeneralMatrix) = eachindex(coefficients(A))

identical(A::T, B::T) where {T<:GeneralMatrix} =
    (coefficients(A) === coefficients(B))

function apply!(α::Number,
                P::Type{<:Operations},
                A::GeneralMatrix{<:AbstractArray{<:GenMult.Floats}},
                x::AbstractArray{<:GenMult.Floats},
                scratch::Bool,
                β::Number,
                y::AbstractArray{<:GenMult.Floats})
    return apply!(α, P, coefficients(A), x, scratch, β, y)
end

function vcreate(P::Type{<:Operations},
                 A::GeneralMatrix{<:AbstractArray{<:GenMult.Floats}},
                 x::AbstractArray{<:GenMult.Floats},
                 scratch::Bool)
    return vcreate(P, coefficients(A), x, scratch)
end

for (T, L) in ((:Direct, 'N'), (:Adjoint, 'C'))
    @eval begin
        function apply!(α::Number,
                        ::Type{$T},
                        A::AbstractArray{<:GenMult.Floats},
                        x::AbstractArray{<:GenMult.Floats},
                        scratch::Bool,
                        β::Number,
                        y::AbstractArray{<:GenMult.Floats})
            return lgemv!(α, $L, A, x, β, y)
        end
    end
end

# To have apply and apply! methods callable with an array (instead of a
# mapping), we have to provide the different possibilities.

apply(A::AbstractArray, x::AbstractArray, scratch::Bool) =
    apply(Direct, A, x, scratch)

apply(P::Type{<:Operations}, A::AbstractArray, x::AbstractArray, scratch::Bool) =
    apply!(1, P, A, x, scratch, 0, vcreate(P, A, x, scratch))

apply!(y::AbstractArray, A::AbstractArray, x::AbstractArray) =
    apply!(1, Direct, A, x, false, 0, y)

apply!(y::AbstractArray, P::Type{<:Operations}, A::AbstractArray, x::AbstractArray) =
    apply!(1, P, A, x, false, 0, y)

function vcreate(P::Type{<:Union{Direct,InverseAdjoint}},
                 A::AbstractArray{Ta,Na},
                 x::AbstractArray{Tx,Nx},
                 scratch::Bool) where {Ta,Na,Tx,Nx}
    # Non-transposed matrix.  Trailing dimensions of X must match those of A,
    # leading dimensions of A are those of the result.  Whatever the scratch
    # parameter, a new array is returned as the operation cannot be done
    # in-place.
    @noinline incompatible_dimensions() =
        bad_size("the indices of `x` do not match the trailing indices of `A`")
    1 ≤ Nx < Na || incompatible_dimensions()
    Ny = Na - Nx
    xinds = axes(x)
    Ainds = axes(A)
    @inbounds for d in 1:Nx
        xinds[d] == Ainds[Ny + d] || incompatible_dimensions()
    end
    shape = ntuple(d -> Ainds[d], Val(Ny)) # faster than Ainds[1:Ny]
    return similar(A, promote_type(Ta, Tx), shape)
end

function vcreate(P::Type{<:Union{Adjoint,Inverse}},
                 A::AbstractArray{Ta,Na},
                 x::AbstractArray{Tx,Nx},
                 scratch::Bool) where {Ta,Na,Tx,Nx}
    # Transposed matrix.  Leading dimensions of X must match those of A,
    # trailing dimensions of A are those of the result.  Whatever the scratch
    # parameter, a new array is returned as the operation cannot be done
    # in-place.
    @noinline incompatible_dimensions() =
        bad_size("the indices of `x` do not match the leading indices of `A`")
    1 ≤ Nx < Na || incompatible_dimensions()
    Ny = Na - Nx
    xinds = axes(x)
    Ainds = axes(A)
    @inbounds for d in 1:Nx
        xinds[d] == Ainds[d] || incompatible_dimensions()
    end
    shape = ntuple(d -> Ainds[Nx + d], Val(Ny)) # faster than Ainds[Nx+1:end]
    return similar(A, promote_type(Ta, Tx), shape)
end


text:
# 1D Shallow water equations with flat bottom
using ConservationLawsDiffEq
using OrdinaryDiffEq

const CFL = 0.1
const Tend = 0.2
const gr = 9.8

function Jf(u::AbstractVector)
  h = u[1]
  q = u[2]
  F =[0.0 1.0;-q^2/h^2+gr*h 2*q/h]
  F
end
f(u::AbstractVector) = [u[2];u[2]^2/u[1]+0.5*gr*u[1]^2]
f0(x) = x < 0.0 ? [2.0,0.0] : [1.0,0.0]

function Nflux(ϕl::AbstractVector, ϕr::AbstractVector)
  hl = ϕl[1]; hr = ϕr[1]
  ul = ϕl[2]/ϕl[1]; ur = ϕr[2]/ϕr[1];
  hm = 0.5*(hl+hr)
  um = 0.5*(ul+ur)
  return([hm*um;hm*um^2+0.5*gr*(0.5*(hl^2+hr^2))])
end
ve(u::AbstractVector) = [gr*u[1]-0.5*(u[2]/u[1])^2;u[2]/u[1]]

# Now discretizate the domain
mesh = Uniform1DFVMesh(100, [-5.0, 5.0])

function get_problem(alg, mesh)
  #Compute discrete initial data
  u0 = getInitialState(mesh,f0,use_threads = true)

  # Now get a explicit semidiscretization (discrete in space) du_h(t)/dt = f_h(u_h(t))
  f_h = getSemiDiscretization(f,alg,mesh,[Periodic()]; Df = Jf, use_threads = false,numvars = 2)

  #Setup ODE problem for a time interval = [0.0,1.0]
  ode_prob = ODEProblem(f_h,u0,(0.0,Tend))

  #Setup callback in order to fix CFL constant value
  cb = getCFLCallback(f_h, CFL)

  #Estimate an initial dt
  dt = update_dt!(u0, f_h, CFL)
  return ode_prob, cb, dt
end

ode_prob, cb, dt = get_problem(FVSKTScheme(), mesh)
sol = solve(ode_prob,SSPRK22(); dt = dt, callback = cb)

#Plot
using Plots
u_h = fv_solution(sol, mesh; vars = 2)
plot(u_h, tidx=1, vars=1, lab="ho",line=(:dot,2))
plot!(u_h, vars=1,lab="KT h")