text:
include("GrpcerrorsClients/GrpcerrorsClients.jl")
using .GrpcerrorsClients

# GrpcerrorsClients.Data mode values:
# 1: throw an error after seconds provided in `param`
# 2: no error, just wait until seconds provided in `param`, respond with SimulationParams
Base.show(io::IO, data::GrpcerrorsClients.Data) = print(io, string("[", data.mode, ", ", data.param, "]"))

# single request, single response
function test_simplerpc(client::GRPCErrorsBlockingClient)
    data = GrpcerrorsClients.Data(; mode=1, param=0)
    try
        _, status_future = GrpcerrorsClients.SimpleRPC(client, data)
        gRPCCheck(status_future)
        error("error not caught")
    catch ex
        @test isa(ex, gRPCServiceCallException)
        @test ex.message == "simulated error mode 1"
        @test ex.grpc_status == StatusCode.UNKNOWN.code
    end

    data = GrpcerrorsClients.Data(; mode=2, param=5)
    try
        _, status_future = GrpcerrorsClients.SimpleRPC(client, data)
        gRPCCheck(status_future)
        error("error not caught")
    catch ex
        @test isa(ex, gRPCServiceCallException)
        @test ex.message == StatusCode.DEADLINE_EXCEEDED.message
        @test ex.grpc_status == StatusCode.DEADLINE_EXCEEDED.code
    end
end

# single request, stream response
function test_stream_response(client::GRPCErrorsBlockingClient)
    data = GrpcerrorsClients.Data(; mode=1, param=0)
    try
        _, status_future = GrpcerrorsClients.StreamResponse(client, data)
        gRPCCheck(status_future)
        error("error not caught")
    catch ex
        @test isa(ex, gRPCServiceCallException)
        @test ex.message == "simulated error mode 1"
        @test ex.grpc_status == StatusCode.UNKNOWN.code
    end

    data = GrpcerrorsClients.Data(; mode=2, param=25)
    try
        _, status_future = GrpcerrorsClients.StreamResponse(client, data)
        gRPCCheck(status_future)
        error("error not caught")
    catch ex
        @test isa(ex, gRPCServiceCallException)
        @test ex.message == StatusCode.DEADLINE_EXCEEDED.message
        @test ex.grpc_status == StatusCode.DEADLINE_EXCEEDED.code
    end
end

# stream request, single response
function test_stream_request(client::GRPCErrorsBlockingClient)
    try
        data = GrpcerrorsClients.Data(; mode=1, param=0)
        in = Channel{GrpcerrorsClients.Data}(1)
        @async begin
            put!(in, data)
            close(in)
        end
        _, status_future = GrpcerrorsClients.StreamRequest(client, in)
        gRPCCheck(status_future)
        error("error not caught")
    catch ex
        @test isa(ex, gRPCServiceCallException)
        @test ex.message == "simulated error mode 1"
        @test ex.grpc_status == StatusCode.UNKNOWN.code
    end

    try
        data = GrpcerrorsClients.Data(; mode=2, param=5)
        in = Channel{GrpcerrorsClients.Data}(1)
        @async begin
            put!(in, data)
            close(in)
        end
        _, status_future = GrpcerrorsClients.StreamRequest(client, in)
        gRPCCheck(status_future)
        error("error not caught")
    catch ex
        @test isa(ex, gRPCServiceCallException)
        @test ex.message == StatusCode.DEADLINE_EXCEEDED.message
        @test ex.grpc_status == StatusCode.DEADLINE_EXCEEDED.code
    end
end

# stream request, stream response
function test_stream_request_response(client::GRPCErrorsBlockingClient)
    try
        data = GrpcerrorsClients.Data(; mode=1, param=0)
        in = Channel{GrpcerrorsClients.Data}(1)
        @async begin
            put!(in, data)
            close(in)
        end
        _, status_future = GrpcerrorsClients.StreamRequestResponse(client, in)
        gRPCCheck(status_future)
        error("error not caught")
    catch ex
        @test isa(ex, gRPCServiceCallException)
        @test ex.message == "simulated error mode 1"
        @test ex.grpc_status == StatusCode.UNKNOWN.code
    end

    try
        data = GrpcerrorsClients.Data(; mode=2, param=5)
        in = Channel{GrpcerrorsClients.Data}(1)
        @async begin
            put!(in, data)
            close(in)
        end
        _, status_future = GrpcerrorsClients.StreamRequestResponse(client, in)
        gRPCCheck(status_future)
        error("error not caught")
    catch ex
        @test isa(ex, gRPCServiceCallException)
        @test ex.message == StatusCode.DEADLINE_EXCEEDED.message
        @test ex.grpc_status == StatusCode.DEADLINE_EXCEEDED.code
    end
end

function test_blocking_client(server_endpoint::String)
    client = GRPCErrorsBlockingClient(server_endpoint; verbose=false, request_timeout=3)
    @testset "request response" begin
        test_simplerpc(client)
    end
    @testset "streaming recv" begin
        test_stream_response(client)
    end
    @testset "streaming send" begin
        test_stream_request(client)
    end
    @testset "streaming send recv" begin
        test_stream_request_response(client)
    end
end

function test_clients(server_endpoint::String)
    @info("testing blocking client")
    test_blocking_client(server_endpoint)
end


text:
#  Contains code that is based in part on Chebfun v5's chebfun/standardChop.m,
# which is distributed with the following license:

# Copyright (c) 2017, The Chancellor, Masters and Scholars of the University
# of Oxford, and the Chebfun Developers. All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
#     * Redistributions of source code must retain the above copyright
#       notice, this list of conditions and the following disclaimer.
#     * Redistributions in binary form must reproduce the above copyright
#       notice, this list of conditions and the following disclaimer in the
#       documentation and/or other materials provided with the distribution.
#     * Neither the name of the University of Oxford nor the names of its
#       contributors may be used to endorse or promote products derived from
#       this software without specific prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR
# ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
# (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
# LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
# ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
# (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

# Jared Aurentz and Nick Trefethen, July 2015.
#
# Copyright 2017 by The University of Oxford and The Chebfun Developers.
# See http://www.chebfun.org/ for Chebfun information.

function standardchoplength(coeffs, tol)
    # Check magnitude of TOL:
    if tol ≥ 1
        throw(ArgumentError("tolerance must be less than 1"))
    end

    # Make sure COEFFS has length at least 17:
    n = length(coeffs)

    if  n < 17
        # resort to naive chop
        mx = maximum(abs,coeffs)
        if mx == 0
            return 0
        end
        for k=n:-1:1
            if abs(coeffs[k]) > tol*mx
                return k
            end
        end
        return 0
    end

    # Step 1: Convert COEFFS to a new monotonically nonincreasing
    #         vector ENVELOPE normalized to begin with the value 1.

    b = abs.(coeffs)
    m = b[end]*fill(1.0,n)
    for j = n-1:-1:1
        m[j] = max(b[j], m[j+1]);
    end
    if  m[1] == 0
        return 1
    end
    envelope = m/m[1];

    plateauPoint = 0

    for j = 2:n
        j2 = round(Int,1.25*j + 5);
        if  j2 > n
            # there is no plateau: exit
            return n
        end
        e1 = envelope[j]
        e2 = envelope[j2]
        r = 3*(1 - log(e1)/log(tol))
        plateau = (e1 == 0) || (e2/e1 > r)
        if plateau
            # a plateau has been found: go to Step 3
            plateauPoint = j - 1
            break
        end
    end

    # Step 3: fix CUTOFF at a point where ENVELOPE, plus a linear function
    # included to bias the result towards the left end, is minimal.
    #

    if  plateauPoint ≠ 0 && envelope[plateauPoint] == 0
        return plateauPoint
    end

    j3 = sum(envelope .≥ tol^(7/6))
    if j3 < j2
        j2 = j3 + 1
        envelope[j2] = tol^(7/6)
    end
    cc = log10.(envelope[1:j2])
    cc .+= range(0, stop=(-1/3)*log10(tol), length=j2)
    d = argmin(cc)
    return max(d - 1, 1)
end

ext: jl


text:
using DataFrames
using LinearAlgebra

export Mack, claims2cum, cum2claims, cum2futureclaims, replacena!
"""
A reserve triangle with chain ladder information
"""
mutable struct Mack
  "Length of history"
  I::Int
  "Cummulative claims (upper triangle)"
  cum::Matrix{Float64}
  "Time vector of development factors"
  f::Vector{Float64}
  "Time vector of future claims payments"
  futureclaims::Vector{Float64}
  "Vector of reseves"
  res::Vector{Float64}
  """
  Total reserves: ``\\mathit{tot\\_res} =\\sum_{i=1}^I \\mathit{res}_i``
  """
  tot_res::Float64
  "Vector of mean square errors"
  mse::Vector{Float64}
  "Total mean square error"
  tot_mse::Float64
end

# function Mack(df_triang::DataFrame; cum::Bool = false)
# Mack(convert(DataArray{Real, 2}, df_triang); cum = cum)
# end
#
# function Mack(da_triang::DataArray{Real, 2}; cum::Bool = false)
#   Mack(convert(Array{Real,2}, da_triang, 0); cum = cum)
# end

function Mack(triang::Matrix{Float64}; cum::Bool = false)
  if cum
    C = deepcopy(triang)
  else
    C = claims2cum(triang)
  end
  # Number of accident / development years
  I = size(C)[1]
  # Development factors
  f = [sum(C[1:(I-𝑘),𝑘+1])/sum(C[1:(I-𝑘),𝑘]) for 𝑘 ∈ 1:(I-1)]
  # Ultimate claim amounts
  for 𝑘  ∈ 2:I
    for 𝑙 ∈ (I+1-𝑘+1):I
      C[𝑘,𝑙] = C[𝑘,I+1-𝑘] * prod(f[(I+1-𝑘):( 𝑙 -1)])
    end
  end
  futclaims = cum2futureclaims(C)
  # Reserves
  R = [C[𝑘,I]-C[𝑘,I+1-𝑘] for 𝑘 ∈ 1:I]
  σ² =
    [1/(I-𝑘-1) *
      C[1:(I-𝑘),𝑘] ⋅ (C[1:(I-𝑘),𝑘+1] ./ C[1:(I-𝑘),𝑘] .- f[𝑘]).^2
    for 𝑘 ∈ 1:(I-1)]
  σ²[I-1] = min(σ²[I-2]^2 / σ²[I-3], min(σ²[I-3],σ²[I-2]))
  # mean square error
  mse = zeros(Float64, I)
  for 𝑖 ∈ 2:I
    mse[𝑖] = 0
    for 𝑘 ∈ (I+1-𝑖):(I-1)
      mse[𝑖] += σ²[𝑘]/f[𝑘]^2 * (1/C[𝑖,𝑘] + 1/ sum(C[1:I-𝑘,𝑘]))
    end
    mse[𝑖] *= C[𝑖,I]^2
  end
  # total mean square error
  mse_total = 0.0
  for 𝑖 ∈ 2:I
    tmp = 0.0
    for 𝑘 ∈ collect(I+1-𝑖 : I-1)
      tmp += (2σ²[𝑘] / f[𝑘]^2) / sum(C[1 : I-𝑘, 𝑘])
    end
    tmp *= C[𝑖,I] * sum(C[𝑖+1 : I, I])
    tmp += mse[𝑖]
    mse_total += tmp
  end
  return Mack(I, C, f, futclaims, R, sum(R), mse, mse_total)
end

"""
replacena!(df::DataFrame, replacement::Any)

Replaces all 'NA'-values with 'replacement'
"""
function replacena!(df::DataFrame, replacement::Any)
  nrows, ncols = size(df)
  for 𝑗 ∈ 1:ncols; for 𝑖 ∈ 1:nrows
    if ismissing(df[𝑖,𝑗]); df[𝑖,𝑗] = replacement; end
    end
  end
end

"""
upperleft(quadrat_mat::Array{Real,2})

Sets all values in the strict lower right triangle to 0.
"""
function upperleft(quadrat_mat::Matrix{T}) where T <: Real
  c = deepcopy(quadrat_mat)
  n = size(c,1)
  if n ≠ size(c,2)
    error("known: Matrix is not quadratic")
  end
  for 𝑗 ∈ 1:n; for 𝑖 ∈ 1:n
    if 𝑖 > n-𝑗+1; c[𝑖,𝑗] = 0 end
  end; end
  c
end

"""
lowerright(quadrat_mat::Array{Real,2})

Sets all values in the upperleft triangle (incl. diagonal) to 0.
"""
function lowerright(quadrat_mat::Matrix{T}) where T <: Real
  quadrat_mat - upperleft(quadrat_mat)
end

"""
    claims2cum(c::Matrix)

    Convert a quadratic upper left triangular reserve matrix to
    the correponding cumulative triangular reserve matrix.
    The (strictly) lower triangular part is set to 0
"""
function claims2cum(c::Matrix{T}) where T <: Real
  if size(c,1) ≠ size(c,2)
    error("claims2cum: Matrix is not quadratic")
  end
 # cum = zeros(Real,size(c))
  #for i in 1:size(c, 1)
    #cum[i,1:size(c, 1)-i+1] = cumsum(c[i,1:size(c, 1)-i+1], 2)
  #end
  upperleft(cumsum(c, dims=2))
end

"""
    cum2claims(c::Matrix)

    Convert a quadratic cumulative claims matrix into a
    claims matrix.
"""
function cum2claims(c::Matrix{T}) where T <: Real
  I = size(c,2)
  if size(c,1) ≠ I
    error("cum2claims: Matrix is not quadratic")
  end
  claims = zero(c)
  for 𝑖 ∈ 1:I
    claims[𝑖,1] = c[𝑖,1]
    for 𝑘 ∈ 2:I
     claims[𝑖,𝑘] = c[𝑖,𝑘]-c[𝑖,𝑘-1]
    end
  end
  claims
end

"""
    cum2futclaims(c::Array{Real,2})

    Calculate the future claims vector from the quadratic,
    (strictly) lower triangle of a projected cumulative claims
    matrix
"""
function cum2futureclaims(c::Matrix{T}) where T <: Real
  I = size(c,2)
  if size(c,1) ≠ I
    error("cumclaims: Matrix is not quadratic")
  end
  claims = cum2claims(c)
  futureclaims = zero(claims[1:(I-1), 1])
   for 𝑘 ∈ 2:I
    futureclaims[𝑘-1] = sum(claims[I-𝑘+2:I,𝑘])
  end
  futureclaims
end


text:
using JSON

demandsDict = Dict{Int64, Dict{(Int64,Int64), Float64}}()
# get ground trueth demands, indexed by 0
demandsDict[0] = iniDemand("../data_original/SiouxFalls_trips_simp.txt")
# get initial demands, indexed by 1
demandsDict[1] = iniDemand("../data_original/SiouxFalls_trips_simp.txt", 1)

numNodes, numLinks, numODpairs, capacity, free_flow_time, ta_data = paraNetwork("Sioux_simp")
numRoutes, odPairRoute, linkRoute, link_label_dict = furInfo()

start_node = ta_data.start_node
end_node = ta_data.end_node

#load OD pair labels
odPairLabel = readall("../temp_files/od_pair_label_dict_Sioux_simp_refined.json")
odPairLabel = JSON.parse(odPairLabel)

odPairLabel_ = readall("../temp_files/od_pair_label_dict__Sioux_simp_refined.json")
odPairLabel_ = JSON.parse(odPairLabel_)

#load node-link incidence
nodeLink = readall("../temp_files/node_link_incidence_Sioux_simp.json");
nodeLink = JSON.parse(nodeLink);

max_stars_repo_name: jingzbu/InverseVITraffic

max_issues_repo_name: jingzbu/InverseVITraffic

max_forks_repo_name: jingzbu/InverseVITraffic


text:
module Test_psanalysis

using PeriodicSystems
using DescriptorSystems
using Symbolics
using Test
using LinearAlgebra
using ApproxFun

println("Test_analysis")

# # selected sequence to illustrate failure of BlockMatrices.jl v0.16.16
# using LinearAlgebra
# using ApproxFun
# Af = Fun(t -> [0 1; -10*cos(t) -24-10*sin(t)],Fourier(0..2π));
# D = Derivative(domain(Af));
# ND = [D 0I; 0I D];
# Aop = Af - ND;
# NA = 322
# RW = Aop[1:NA,1:NA]

@testset "test_analysis" begin

# using Floquet based approach
At = PeriodicFunctionMatrix(t -> [0 1; -10*cos(t) -24-10*sin(t)],2pi);
psys = ps(At,HarmonicArray(rand(2,1),2*pi),rand(1,2))
@time ev = pspole(psys,500)
@test sort(ev) ≈ sort([0;-24])

# using Fourier series
Afun = FourierFunctionMatrix(Fun(t -> [0 1; -10*cos(t) -24-10*sin(t)],Fourier(0..2π)))
psys = ps(Afun,FourierFunctionMatrix(rand(2,1),2*pi),rand(1,2))
@time ev1 = pspole(psys,50)
@test isapprox(sort(real(ev)),sort(real(ev1)),rtol=1.e-6) && norm(imag(ev1)) < 1.e-10

# using Toeplitz operator truncation
Ahr = convert(HarmonicArray,PeriodicFunctionMatrix(t -> [0 1; -10*cos(t) -24-10*sin(t)],2pi));
psys = ps(Ahr,rand(2,1),rand(1,2))
@time ev2 = pspole(psys,50)
@test sort(ev) ≈ sort(real(ev2[sortperm(imag(ev2),by=abs)][1:2]))

# example Zhou, Hagiwara SCL 2002 period pi/2 and pi
# using Floquet based approach
At = PeriodicFunctionMatrix(t -> [-1-sin(2*t)^2 2-0.5*sin(4*t); -2-0.5*sin(4*t) -1-cos(2*t)^2],pi/2);
At2 = PeriodicFunctionMatrix(t -> [-1-sin(2*t)^2 2-0.5*sin(4*t); -2-0.5*sin(4*t) -1-cos(2*t)^2],pi);
psys = ps(At,rand(2,1),rand(1,2))
psys2 = ps(At2,rand(2,1),rand(1,2))
@time ev = pspole(psys; reltol = 1.e-10)
@test sort(real(ev)) ≈ sort([-1;-2]) && sort(imag(ev)) ≈ [2;2]
@time ev = pspole(psys, 500; reltol = 1.e-10)
@test sort(real(ev)) ≈ sort([-1;-2]) && sort(imag(ev)) ≈ [2;2]
@time ev2 = pspole(psys2; reltol = 1.e-10)
@test sort(ev2) ≈ sort([-1;-2])
@time ev2 = pspole(psys2,500; reltol = 1.e-10)
@test sort(ev2) ≈ sort([-1;-2])

# using Toeplitz operator truncation period pi/2
Ahr = convert(HarmonicArray,PeriodicFunctionMatrix(t -> [-1-sin(2*t)^2 2-0.5*sin(4*t); -2-0.5*sin(4*t) -1-cos(2*t)^2],pi/2));
psys = ps(Ahr,rand(2,0),rand(0,2))
ev3 = pspole(psys,5)
@test sort(real(ev3)) ≈ sort([-1;-2]) && sort(imag(ev3)) ≈ [2;2]

# simple period
ev4 = pspole(psys,50)
@test sort(real(ev4)) ≈ sort([-1;-2]) && sort(imag(ev4)) ≈ [2;2]

# double period
ev5 = pspole(psys,50; P = 2)
@test sort(real(ev5)) ≈ sort([-1;-2]) && norm(imag(ev5)) < 1.e-10

# using Fourier series
Afun = FourierFunctionMatrix(Fun(t -> [-1-sin(2*t)^2 2-0.5*sin(4*t); -2-0.5*sin(4*t) -1-cos(2*t)^2],Fourier(0..π/2)))
psys = ps(Afun,rand(2,1),rand(1,2))
ev3 = pspole(psys)
@test sort(real(ev3)) ≈ sort([-1;-2]) && sort(imag(ev3)) ≈ [2;2]

# using Fourier series truncation period pi
ev4 = pspole(psys, P = 2)
@test sort(real(ev4)) ≈ sort([-1;-2]) && norm(imag(ev4)) < 1.e-10

# using Toeplitz operator truncation period pi
Ahr2 = convert(HarmonicArray,PeriodicFunctionMatrix(t -> [-1-sin(2*t)^2 2-0.5*sin(4*t); -2-0.5*sin(4*t) -1-cos(2*t)^2],pi));
psys2 = ps(Ahr2,rand(2,0),rand(0,2))
ev3 = pspole(psys2,5)
@test sort(real(ev3)) ≈ sort([-1;-2]) && norm(imag(ev3)) < 1.e-10

ev4 = pspole(psys2,50; P = 2)
@test sort(real(ev4)) ≈ sort([-1;-2]) && norm(imag(ev4)) < 1.e-10

# constant system case
Q = convert(HarmonicArray,PeriodicFunctionMatrix([-1 0; 0 -2],pi/2));
psys = ps(Q,rand(2,0),rand(0,2))
ev1 = pspole(psys,5)
@test sort(real(ev1)) ≈ sort([-1;-2]) && norm(imag(ev1)) < 1.e-10

QF = convert(FourierFunctionMatrix,Q)
#psys = ps(QF,rand(2,0),rand(0,2)) # null dimensions not supported in ApproxFun
psys = ps(QF,rand(2,1),rand(1,2))
ev2 = pspole(psys,5)
@test sort(real(ev2)) ≈ sort([-1;-2]) && norm(imag(ev2)) < 1.e-10

# Floquet approach
A = PeriodicFunctionMatrix(t -> [-1-sin(2*t)^2 2-0.5*sin(4*t); -2-0.5*sin(4*t) -1-cos(2*t)^2],pi);
psys = ps(A,rand(2,1),rand(1,2))
@time ev = pspole(psys; reltol = 1.e-10)
@time ev = pspole(psys,50; reltol = 1.e-10)

# Fourier approach
Afun = FourierFunctionMatrix(Fun(t -> [-1-sin(2*t)^2 2-0.5*sin(4*t); -2-0.5*sin(4*t) -1-cos(2*t)^2],Fourier(0..π)))
psys1 = ps(Afun,rand(2,1),rand(1,2))
ev1 = pspole(psys1)
@test sort(real(ev)) ≈ sort(real(ev1)) && norm(imag(ev1)) < 1.e-10

ev2 = pspole(psys1,5)
@test sort(real(ev)) ≈ sort(real(ev2)) && norm(imag(ev2)) < 1.e-10

ev3 = pspole(psys1,20)
@test sort(real(ev)) ≈ sort(real(ev3)) && norm(imag(ev3)) < 1.e-10

# # lossy Mathieu differential equation
# k = -1.; ξ = 0.05; β =0.2; ωh = 2; # unstable feedback
# k = -.5; ξ = 0.05; β =0.2; ωh = 2; # stable feedback
# Ahr = convert(HarmonicArray,PeriodicFunctionMatrix(t -> [0 1; k*(1-2β*cos(ωh*t)) -2*ξ],2pi));
# Asym = convert(PeriodicSymbolicMatrix,Ahr); Asym.F
# N = 5; ev = eigvals(hr2bt(Ahr,N)-phasemat(Ahr,N))
# ev[sortperm(imag(ev),by=abs)][1:2]

# zeros example: infinite zeros present but not for the averaged system
a = PeriodicFunctionMatrix(t -> [sin(t)],2*pi);
b = PeriodicFunctionMatrix(t ->  [cos(t)], 2*pi);
c = PeriodicFunctionMatrix(t ->  [cos(t)+sin(t)], 2*pi);
d = [0];
psys = ps(HarmonicArray,a,b,c,d);
z = pszero(psys,atol=1.e-7)
@test z[1] == Inf
psys1 = ps(FourierFunctionMatrix,a,b,c,d)
z1 = pszero(psys1,30,atol=1.e-7)
@test any(isinf.(z1))

psys2 = ps(a,b,c,d);
z2 = pszero(psys2,30,atol=1.e-7)
@test any(isinf.(z2))

# Bitanti-Colaneri's book p.26
a = PeriodicFunctionMatrix(t -> [-1+sin(t) 0; 1-cos(t) -3],2*pi);
b = PeriodicFunctionMatrix(t ->  [-1-cos(t);2-sin(t)], 2*pi);
c = [0 1]; d = [0];

psys = ps(a,b,c,d)
evref = pspole(psys,100)
@test sort(real(evref)) ≈ sort([-1,-3])  && norm(imag(evref)) < 1.e-10

# using Harmonic Array based computation
psyschr = ps(HarmonicArray,a,b,c,d)
ev = pspole(psyschr,20)
@test sort(real(evref)) ≈ sort(real(ev))  && norm(imag(ev)) < 1.e-10

Z = convert(HarmonicArray,PeriodicFunctionMatrix(t -> [(-2+3*sin(t))/(2-sin(t))],2*pi))
ρ = real(Z.values[:,:,1])[1,1]

z = pszero(psyschr,atol=1.e-7)
@test z ≈ [ρ;Inf]

# using Fourier Function Matrix based computation
psysc = ps(FourierFunctionMatrix,a,b,c,d);
ev1 = pspole(psysc,30)
@test sort(real(evref)) ≈ sort(real(ev1))  && norm(imag(ev1)) < 1.e-10

z1 = pszero(psysc,atol=1.e-7)
@test minimum(abs.(z1.-ρ)) < 1.e-10  && any(isinf.(z1))

psys2 = ps(a,b,c,d);
z2 = pszero(psys2,30,atol=1.e-7)
@test z2 ≈ [ρ;Inf]

# Zhou-Hagiwara Automatica 2002
β = 0.5
a1 = PeriodicFunctionMatrix(t -> [-1-sin(2*t)^2 2-0.5*sin(4*t); -2-0.5*sin(4*t) -1-cos(2*t)^2],pi);
#γ(t) = mod(t,pi) < pi/2 ? sin(2*t) : 0
γ = chop(Fun(t -> mod(t,pi) < pi/2 ? sin(2*t) : 0, Fourier(0..pi)),1.e-7);
#b1 = PeriodicFunctionMatrix(t ->  [0; 1-2*β*(mod(t,float(pi)) < pi/2 ? sin(2*t) : 0 )], pi);
b1 = PeriodicFunctionMatrix(t ->  [0; 1-2*β*γ(t)], pi);
c = [1 1]; d = [0];

# using Harmonic Array based lifting
psyschr = ps(HarmonicArray,a1,b1,c,d);
ev = pspole(psyschr)
@test sort([-1,-2]) ≈ sort(real(ev))  && norm(imag(ev)) < 1.e-10

z = pszero(psyschr,atol=1.e-7)
@test z ≈ [-3.5;Inf]

# using Fourier Function Matrix based lifting
psysc = ps(FourierFunctionMatrix,a1,b1,c,d);
ev = pspole(psysc)
@test sort([-1,-2]) ≈ sort(real(ev))  && norm(imag(ev)) < 1.e-10

z = pszero(psysc,30,atol=1.e-7)
@test z ≈ [-3.5;Inf]

psys2 = ps(a1,b1,c,d);
z2 = pszero(psys2,30,atol=1.e-7)
@test z2 ≈  [-3.5;Inf]

# Zhou IEEE TAC 2009
β = 0.5
a1 = PeriodicFunctionMatrix(t -> [-1-sin(2*t)^2 2-0.5*sin(4*t); -2-0.5*sin(4*t) -1-cos(2*t)^2],pi);
b1 = PeriodicFunctionMatrix(t ->  [0; 1-2*β*cos(2*t)], pi);
c = [0 1]; d = [0];

# using Harmonic Array based lifting
psyschr = ps(HarmonicArray,a1,b1,c,d);
ev = pspole(psyschr)
@test sort([-1,-2]) ≈ sort(real(ev))  && norm(imag(ev)) < 1.e-10

z = pszero(psyschr,atol=1.e-7)
@test z ≈ [-1.5;Inf]

psysc = ps(FourierFunctionMatrix,a1,b1,c,d);
ev = pspole(psysc)
@test sort([-1,-2]) ≈ sort(real(ev))  && norm(imag(ev)) < 1.e-10

z = pszero(psysc,30,atol=1.e-7)
@test z ≈ [-1.5;Inf]

psys2 = ps(a1,b1,c,d);
z2 = pszero(psys2,30,atol=1.e-7)
@test z2 ≈  [-1.5;Inf]

# # Ziegler's column

# β = 2pi; η =0.5; λ = 0.5;
# M = [1 3/8; 3/2 1];
# K0 = [3/8 3/16; -3/4 3/4]; K1 = λ*[-1 η; 0 4η-4];
# A = PeriodicFunctionMatrix(τ -> [eye(2) zeros(2,2); zeros(2,2) -M\(K0+cos(β*τ).*K1)],2*pi/β)

# # Mathieu equation: determination of minimum eigenvalue
# q = 2; a = 2; a = -1.513956885056448; a = -1
# Ahr = convert(HarmonicArray,PeriodicFunctionMatrix(t -> [0 1; a-2*q*cos(2*t) 0],pi));
# Asym = convert(PeriodicSymbolicMatrix,Ahr); Asym.F
# ev = psceig(Ahr)
# N = 20; ev = eigvals(hr2bt(Ahr,N)-phasemat(Ahr,N))
# ev[sortperm(imag(ev),by=abs)][1:2]

# N = 5; ev1 = eigvals(hr2bt1(Ahr,N,2)-phasemat1(Ahr,N,2))
# ev1[sortperm(imag(ev1),by=abs)][1:2]

# discrete systems
Ad = PeriodicMatrix([[8. 0; 0 0], [1 1;1 1], [0 1; 1 0]], 6, nperiod = 2);
Bd = PeriodicMatrix( [[ 1; 0 ], [ 1; 1]] ,2);
Cd = PeriodicMatrix( [[ 1 1 ], [ 1 0]] ,2);
Dd = PeriodicMatrix( [[ 1 ]], 1);
psys = PeriodicStateSpace(Ad,Bd,Cd,Dd);
ev = pspole(psys)
ev1 = pspole(psys,7)
@test sort(ev) ≈ sort([2,0]) ≈ sort(ev1) ≈ sort(gpole(ps2ls(psys))).^(1/6)

z = pszero(psys)
z1 = gzero(ps2ls(psys))
@test all(abs.(z[isfinite.(z)].^6-z1[isfinite.(z1)]) .< 1.e-8) && length(z[isinf.(z)]) == length(z1[isinf.(z1)])

Ad = PeriodicMatrix([[4. 0], [1;1], [2 0;0 1]],3);
Bd = PeriodicMatrix( [[ 1 ], [ 1; 1], [1;2]] ,3);
Cd = PeriodicMatrix( [[ 1 1 ], [ 1 ], [1 0]] ,3);
Dd = PeriodicMatrix( [[ 1 ]], 1);
psys = PeriodicStateSpace(Ad,Bd,Cd,Dd);
ev = pspole(psys)
@test sort(ev) ≈ sort([2,0])
ev1 = pspole(psys,2)
@test ev1 ≈ [2] ≈ gpole(ps2ls(psys,2)).^(1/3)
@test !isstable(psys)

z = pszero(psys)
z1 = gzero(ps2ls(psys))
@test all(abs.(z[isfinite.(z)].^3-z1[isfinite.(z1)]) .< 1.e-8) && length(z[isinf.(z)]) == length(z1[isinf.(z1)])

Ad = PeriodicMatrix([[1. 2], [7;8]],2);
Bd = PeriodicMatrix( [[ 3 ], [ 9; 10]] ,2);
Cd = PeriodicMatrix( [[ 4 5 ], [ 11 ]] ,2);
Dd = PeriodicMatrix( [[ 0 ]], 1);
psys = PeriodicStateSpace(Ad,Bd,Cd,Dd);

z = pszero(psys)
z1 = gzero(ps2ls(psys))
@test all(abs.(z[isfinite.(z)].^2-z1[isfinite.(z1)]) .< 1.e-8) && length(z[isinf.(z)]) == length(z1[isinf.(z1)])

z = pszero(psys,2)
z1 = gzero(ps2ls(psys,2))
@test all(abs.(z[isfinite.(z)].^2-z1[isfinite.(z1)]) .< 1.e-8) && length(z[isinf.(z)]) == length(z1[isinf.(z1)])

psysa = convert(PeriodicStateSpace{PeriodicArray},psys)
za = pszero(psysa)
za1 = gzero(ps2ls(psysa))
@test all(abs.(za[isfinite.(za)].^2-za1[isfinite.(za1)]) .< 1.e-8) && length(za[isinf.(za)]) == length(za1[isinf.(za1)])

za = pszero(psysa,2)
za1 = gzero(ps2ls(psysa,2))
@test all(abs.(za[isfinite.(za)].^2-za1[isfinite.(za1)]) .< 1.e-8) && length(za[isinf.(za)]) == length(za1[isinf.(za1)])

# periodic array representation
Ad = PeriodicArray(rand(Float32,2,2,10),10);
Bd = PeriodicArray(rand(2,1,2),2);
Cd = PeriodicArray(rand(1,2,3),3);
Dd = PeriodicArray(zeros(1,1,1), 1);
psys = PeriodicStateSpace(Ad,Bd,Cd,Dd);
ev = pspole(psys)
ev1 = pspole(psys; fast = true)
@test sort(real(ev)) ≈ sort(real(ev1)) && norm(sort(imag(ev))-sort(imag(ev1))) < 1.e-7

z = pszero(psys,atol=1.e-7)
z1 = gzero(ps2ls(psys),atol=1.e-7)
@test all(abs.(z[isfinite.(z)].^30-z1[isfinite.(z1)]) .< 1.e-8) && length(z[isinf.(z)]) == length(z1[isinf.(z1)])

# Pitelkau's example
ω = 0.00103448
period = 2*pi/ω
a = [  0            0     5.318064566757217e-02                         0
       0            0                         0      5.318064566757217e-02
       -1.352134256362805e-03   0             0    -7.099273035392090e-02
       0    -7.557182037479544e-04     3.781555288420663e-02             0
];
b = t->[0; 0; 0.1389735*10^-6*sin(ω*t); -0.3701336*10^-7*cos(ω*t)];
c = [1 0 0 0;0 1 0 0];
d = zeros(2,1);
psysc = ps(a,PeriodicFunctionMatrix(b,period),c,d);
@time evc = pspole(psysc)
@test all(abs.(real(evc)) .< 1.e-10)
@test !isstable(psysc)

zc = pszero(psysc)
@test all(isinf.(zc))

K = 120;
@time psys = psc2d(psysc,period/K,reltol = 1.e-10);
@time evd = pspole(psys)
@test all(abs.(evd) .≈ 1)
@test !isstable(psys)

zd = pszero(psys,atol=1.e-7)
@test all(isinf.(zd))

end # test

end # module


text:
@info("Basic test of PIPotential construction and evaluation")
maxdeg = 10
r0 = 1.0
rcut = 3.0
trans = PolyTransform(1, r0)
Pr = transformed_jacobi(maxdeg, trans, rcut; pcut = 2)
D = ACE.SparsePSHDegree()
P1 = ACE.RnYlm1pBasis(Pr; species = :X, D = D)
basis = ACE.PIBasis(P1, 2, D, maxdeg)
c = ACE.Random.randcoeffs(basis)
V = combine(basis, c)
Nat = 15
Rs, Zs, z0 = ACE.rand_nhd(Nat, Pr, :X)
val_basis = real(sum(c .* evaluate(basis, Rs, Zs, z0)))
val_V = evaluate(V, Rs, Zs, z0)
println(@test(val_basis ≈ val_V))
J = evaluate_d(basis, Rs, Zs, z0)
grad_basis = real(sum(c[i] * J[i,:] for i = 1:length(c)))[:]
grad_V = evaluate_d(V, Rs, Zs, z0)
println(@test(grad_basis ≈ grad_V))

#---
# using BenchmarkTools
# Vgr = ACE.GraphPIPot(V)
# tmp = ACE.alloc_temp(V, Nat)
# tmpgr = ACE.alloc_temp(Vgr, Nat)
# @btime ACE.evaluate!($tmp, $V, $Rs, $Zs, $z0)
# @btime ACE.evaluate!($tmpgr, $Vgr, $Rs, $Zs, $z0)
# @btime ACE.evaluate_old!($tmpgr, $Vgr, $Rs, $Zs, $z0)
#
# #---
#
# @code_warntype ACE.evaluate!(tmpgr, Vgr, Rs, Zs, z0)


text:
using Gnuplot, ColorSchemes, RDatasets
volcano = Matrix{Float64}(dataset("datasets", "volcano"))
@gp title = "Auckland s Maunga Whau Volcano" :-
@gp :- volcano "w image notit" "set auto fix" "set size ratio -1" palette(:inferno)
save(term="pngcairo size 900,600", output="./code/plt3d_ex6.png") # hide

max_stars_repo_name: lazarusA/gnuplot-examples

max_issues_repo_name: lazarusA/gnuplot-examples

max_forks_repo_name: lazarusA/gnuplot-examples


text:
using TauPy
using Test

@testset "Phase" begin
    @testset "Time" begin
        let p = travel_time(300, 60, "PcS", model="prem")
            @test p isa Vector{<:Phase}
            @test length(p) == 1
            @test p[1].model == "prem"
            @test p[1].name == "PcS"
            @test p[1].delta ≈ 60.0
            @test p[1].depth ≈ 300.0
            @test p[1].time ≈ 860.2281352411015
            @test p[1].dtdd ≈ 4.425149199692147
            @test p[1].inc ≈ 7.316383583264879
            @test p[1].takeoff ≈ 21.338284647101577
            @test p[1].pierce ≈ Float64[]
            @test p[1].distance ≈ Float64[]
            @test p[1].radius ≈ Float64[]
        end

        # Triplications around 20° should many arrivals
        let p = travel_time(0, 20, "P", model="iasp91")
            @test length(p) == 5
            @test all([pp.name for pp in p] .== "P")
            @test all([pp.model for pp in p] .== "iasp91")
            @test p[end].delta ≈ 20.0
            @test p[end].depth ≈ 0.0
            @test p[end].time ≈ 279.8555618981382
            @test p[end].dtdd ≈ 9.484314474162678
            @test p[end].inc ≈ 29.65049984974929
            @test p[end].takeoff ≈ 29.65049984974929
            @test p[end].pierce ≈ Float64[]
            @test p[end].distance ≈ Float64[]
            @test p[end].radius ≈ Float64[]
        end

        let h = 0, Δ = 60
            pp = travel_time(h, Δ, "P")
            ps = travel_time(h, Δ, "S")
            p = travel_time(h, Δ, ["P", "S"])
            @test length(p) == 2
            @test p[1].name == "P"
            @test p[2].name == "S"
            for f in fieldnames(Phase)
                @test getfield(pp[1], f) == getfield(p[1], f)
                @test getfield(ps[1], f) == getfield(p[2], f)
            end
        end
    end

    @testset "Path" begin
        let p = travel_time(100, 60, "ScS")[1]
            @test_throws ErrorException turning_depth(p)
        end

        let p = path(100, 120, "PKiKP")
            @test length(p) == 1
            @test p[1].model == "ak135"
            @test p[1].name == "PKiKP"
            @test p[1].delta ≈ 120.0
            @test p[1].depth ≈ 100.0
            @test p[1].time ≈ 1118.8440838099027
            @test p[1].dtdd ≈ 1.9568726794035047
            @test p[1].inc ≈ 5.858483989177137
            @test p[1].takeoff ≈ 8.272751909308088
            @test p[1].pierce ≈ Float64[]
            distance = [0.0, 0.000731725, 0.0266589, 0.029084, 0.0315151, 0.0363955, 0.0462292, 0.0651594, 0.0844533, 0.0861091, 0.0877676, 0.105556, 0.123657, 0.150756, 0.156101, 0.161479, 0.178193, 0.195233, 0.223203, 0.226581, 0.229972, 0.247623, 0.256577, 0.265619, 0.281917, 0.298502, 0.300332, 0.302166, 0.320744, 0.330169, 0.339685, 0.358042, 0.376744, 0.377477, 0.378211, 0.397844, 0.407804, 0.417862, 0.43775, 0.458023, 0.464182, 0.470373, 0.482851, 0.507072, 0.531778, 0.545519, 0.581614, 0.634116, 0.636169, 0.68739, 0.730078, 0.743049, 0.800721, 0.827352, 0.8605, 0.921914, 0.928103, 0.935293, 0.94251, 0.957024, 0.99736, 1.03852, 1.04006, 1.04159, 1.08617, 1.13174, 1.15241, 1.19942, 1.24722, 1.26839, 1.35612, 1.38744, 1.46901, 1.50969, 1.59274, 1.59538, 1.59802, 1.6033, 1.6139, 1.62454, 1.63522, 1.65104, 1.66695, 1.69901, 1.73139, 1.76411, 1.7967, 1.82961, 1.86286, 1.89645, 1.96372, 2.03233, 2.10141, 2.17186, 2.24279, 2.31512, 2.39346, 2.39615, 2.39884, 2.40423, 2.41503, 2.43673, 2.45856, 2.46223, 2.48066, 2.52532, 2.54784, 2.5705, 2.6133, 2.6969, 2.76179, 2.76844, 2.82813, 2.92778, 2.9632, 3.03275, 3.09144, 3.10345, 3.24376, 3.25956, 3.31648, 3.39062, 3.43229, 3.53735, 3.60977, 3.6125, 3.68889, 3.79216, 3.84144, 3.92027, 3.97963, 4.00037, 4.15911, 4.17234, 4.24052, 4.32324, 4.37049, 4.48792, 4.57427, 4.65828, 4.78389, 4.99127, 4.99956, 5.00715, 5.22151, 5.3673, 5.45001, 5.6685, 5.68529, 5.73935, 5.92767, 6.12544, 6.1774, 6.41211, 6.43481, 6.52613, 6.70021, 6.94192, 6.97396, 7.24102, 7.25646, 7.37187, 7.5481, 7.82, 7.84929, 8.16052, 8.1682, 8.28289, 8.4823, 8.76297, 8.80831, 9.15608, 9.20474, 9.26455, 9.51632, 9.60282, 9.68997, 9.86622, 9.88627, 10.0755, 10.2678, 10.5686, 10.6266, 10.682, 10.6852, 10.6885, 10.6951, 10.7082, 10.7345, 10.761, 10.7877, 10.8414, 10.8685, 10.8957, 11.0106, 11.0695, 11.1293, 11.2493, 11.3108, 11.3733, 11.4986, 11.5628, 11.6281, 11.7591, 11.8262, 11.8944, 12.0312, 12.1013, 12.1725, 12.3156, 12.3888, 12.4632, 12.6128, 12.7671, 13.0705, 13.078, 13.0849, 13.2829, 13.4174, 13.4947, 13.6989, 13.7146, 13.7652, 13.9436, 14.1293, 14.1787, 14.4023, 14.4239, 14.5107, 14.6777, 14.9103, 14.9414, 15.2022, 15.2173, 15.3291, 15.5013, 15.7686, 15.7977, 16.1074, 16.1151, 16.2299, 16.4313, 16.7145, 16.7606, 17.1137, 17.1632, 17.2241, 17.4822, 17.5509, 17.6202, 17.7604, 17.7765, 17.9285, 18.0832, 18.3253, 18.3726, 18.4177, 18.433, 18.4483, 18.4791, 18.51, 18.541, 18.6033, 18.666, 18.7291, 18.7926, 18.921, 19.0745, 19.2304, 19.55, 19.7851, 19.9976, 20.2146, 20.5605, 20.918, 21.4383, 21.6636, 22.4552, 23.2971, 23.391, 24.1944, 25.1528, 25.7264, 26.1789, 27.2809, 28.4675, 28.5628, 29.7498, 31.1414, 32.6584, 34.3156, 36.156, 36.5983, 38.1971, 40.4906, 42.7234, 43.1012, 46.1224, 49.7054, 52.2481, 54.1233, 59.9402, 65.7571, 67.6323, 70.175, 73.758, 76.7791, 77.157, 79.3897, 81.6832, 83.282, 83.7243, 85.5647, 87.222, 88.739, 90.1305, 91.3175, 91.4129, 92.5995, 93.7014, 94.1539, 94.7276, 95.686, 96.4894, 96.5833, 97.4252, 98.2168, 98.4421, 98.9623, 99.3199, 99.6657, 99.8828, 100.095, 100.33, 100.65, 100.806, 100.959, 101.088, 101.151, 101.214, 101.277, 101.339, 101.37, 101.401, 101.432, 101.447, 101.463, 101.508, 101.555, 101.797, 101.952, 102.104, 102.12, 102.26, 102.329, 102.398, 102.656, 102.717, 102.767, 103.12, 103.166, 103.449, 103.65, 103.765, 103.773, 104.083, 104.112, 104.379, 104.551, 104.663, 104.678, 104.939, 104.97, 105.203, 105.37, 105.456, 105.478, 105.702, 105.751, 105.937, 106.115, 106.166, 106.181, 106.386, 106.463, 106.597, 106.795, 106.802, 106.81, 107.113, 107.268, 107.417, 107.492, 107.565, 107.708, 107.779, 107.849, 107.986, 108.054, 108.121, 108.252, 108.318, 108.382, 108.507, 108.57, 108.631, 108.751, 108.811, 108.87, 108.985, 109.012, 109.039, 109.093, 109.119, 109.146, 109.172, 109.185, 109.192, 109.195, 109.198, 109.254, 109.312, 109.613, 109.805, 109.994, 110.014, 110.19, 110.278, 110.364, 110.616, 110.676, 110.724, 111.072, 111.117, 111.398, 111.597, 111.712, 111.72, 112.031, 112.06, 112.332, 112.508, 112.624, 112.639, 112.906, 112.938, 113.18, 113.354, 113.446, 113.468, 113.703, 113.755, 113.953, 114.141, 114.195, 114.212, 114.43, 114.513, 114.659, 114.873, 114.881, 114.889, 115.096, 115.222, 115.306, 115.392, 115.51, 115.557, 115.64, 115.708, 115.721, 115.88, 115.901, 115.96, 116.039, 116.088, 116.191, 116.268, 116.271, 116.343, 116.448, 116.49, 116.564, 116.621, 116.637, 116.777, 116.789, 116.848, 116.917, 116.953, 117.052, 117.112, 117.119, 117.183, 117.267, 117.31, 117.333, 117.355, 117.4, 117.418, 117.422, 117.444, 117.465, 117.476, 117.482, 117.484, 117.487, 117.565, 117.638, 117.708, 117.779, 117.848, 117.917, 117.984, 118.017, 118.051, 118.084, 118.116, 118.149, 118.181, 118.213, 118.229, 118.245, 118.256, 118.266, 118.277, 118.282, 118.285, 118.288, 118.371, 118.411, 118.493, 118.524, 118.612, 118.633, 118.681, 118.728, 118.749, 118.794, 118.839, 118.84, 118.842, 118.883, 118.923, 118.938, 118.945, 118.952, 118.958, 119.02, 119.053, 119.08, 119.137, 119.15, 119.193, 119.244, 119.246, 119.299, 119.335, 119.349, 119.373, 119.398, 119.41, 119.416, 119.422, 119.443, 119.462, 119.473, 119.483, 119.502, 119.503, 119.504, 119.522, 119.541, 119.55, 119.56, 119.578, 119.58, 119.582, 119.598, 119.615, 119.624, 119.633, 119.65, 119.654, 119.657, 119.685, 119.702, 119.719, 119.724, 119.73, 119.757, 119.775, 119.793, 119.794, 119.796, 119.815, 119.834, 119.844, 119.849, 119.851, 119.854, 119.88, 119.88, 119.893, 119.905, 119.91, 119.924, 119.938, 119.945, 119.952, 119.959, 119.962, 119.964, 119.966, 119.974, 119.978, 119.979, 119.98, 119.981, 119.989, 119.994, 119.997, 119.999, 120.0]
            @test p[1].distance ≈ distance atol=0.01
            radius = [6271.0, 6270.45, 6251.0, 6249.19, 6247.37, 6243.74, 6236.44, 6222.49, 6208.41, 6207.2, 6206.0, 6193.15, 6180.19, 6161.0, 6157.24, 6153.48, 6141.84, 6130.07, 6111.0, 6108.72, 6106.43, 6094.57, 6088.6, 6082.6, 6071.85, 6061.0, 6059.81, 6058.62, 6046.59, 6040.53, 6034.44, 6022.78, 6011.0, 6010.54, 6010.08, 5997.84, 5991.67, 5985.47, 5973.29, 5961.0, 5957.42, 5953.83, 5946.63, 5932.75, 5918.74, 5911.0, 5890.88, 5862.11, 5861.0, 5833.51, 5811.0, 5804.23, 5774.52, 5761.0, 5744.36, 5714.02, 5711.0, 5707.7, 5704.39, 5697.75, 5679.46, 5661.0, 5660.32, 5659.63, 5639.89, 5619.96, 5611.0, 5590.78, 5570.44, 5561.5, 5524.89, 5512.0, 5478.82, 5462.5, 5429.6, 5428.57, 5427.53, 5425.46, 5421.31, 5417.16, 5413.0, 5406.85, 5400.69, 5388.34, 5375.94, 5363.5, 5351.19, 5338.83, 5326.44, 5314.0, 5289.33, 5264.5, 5239.83, 5215.0, 5190.33, 5165.5, 5138.98, 5138.08, 5137.17, 5135.36, 5131.74, 5124.49, 5117.22, 5116.0, 5109.89, 5095.17, 5087.79, 5080.39, 5066.5, 5039.66, 5019.1, 5017.0, 4998.3, 4967.5, 4956.67, 4935.6, 4918.0, 4914.42, 4873.09, 4868.5, 4852.04, 4830.82, 4819.0, 4789.53, 4769.5, 4768.75, 4747.86, 4720.0, 4706.85, 4686.03, 4670.5, 4665.11, 4624.35, 4621.0, 4603.81, 4583.18, 4571.5, 4542.81, 4522.0, 4501.99, 4472.5, 4424.88, 4423.0, 4421.28, 4373.5, 4341.75, 4324.0, 4277.99, 4274.5, 4263.33, 4225.0, 4185.68, 4175.5, 4130.3, 4126.0, 4108.8, 4076.5, 4032.71, 4027.0, 3980.17, 3977.5, 3957.72, 3928.0, 3883.24, 3878.5, 3829.0, 3827.8, 3809.98, 3779.5, 3737.65, 3731.0, 3681.0, 3674.14, 3665.76, 3631.0, 3619.25, 3607.49, 3583.98, 3581.33, 3556.5, 3531.67, 3493.61, 3486.37, 3479.5, 3478.76, 3478.01, 3476.52, 3473.55, 3467.59, 3461.62, 3455.64, 3443.67, 3437.67, 3431.67, 3406.65, 3394.03, 3381.34, 3356.32, 3343.7, 3331.01, 3305.98, 3293.37, 3280.68, 3255.64, 3243.02, 3230.34, 3205.3, 3192.69, 3180.01, 3154.97, 3142.35, 3129.68, 3104.63, 3079.35, 3031.26, 3030.09, 3029.02, 2998.77, 2978.69, 2967.31, 2937.78, 2935.54, 2928.36, 2903.42, 2878.03, 2871.39, 2841.76, 2838.94, 2827.7, 2806.37, 2777.36, 2773.54, 2742.03, 2740.24, 2727.03, 2707.03, 2676.7, 2673.45, 2639.5, 2638.67, 2626.37, 2605.15, 2576.04, 2571.39, 2536.4, 2531.59, 2525.71, 2501.16, 2494.73, 2488.29, 2475.38, 2473.91, 2460.12, 2446.3, 2425.05, 2420.96, 2417.07, 2415.75, 2414.43, 2411.8, 2409.16, 2406.52, 2401.23, 2395.94, 2390.64, 2385.34, 2374.72, 2362.18, 2349.61, 2324.38, 2306.26, 2290.18, 2274.05, 2248.94, 2223.72, 2188.29, 2173.39, 2123.06, 2072.73, 2067.31, 2022.4, 1972.07, 1943.5, 1921.74, 1871.4, 1821.07, 1817.19, 1770.74, 1720.41, 1670.08, 1619.91, 1569.42, 1558.03, 1519.09, 1468.76, 1425.29, 1418.42, 1368.09, 1317.76, 1287.35, 1267.43, 1217.5, 1267.43, 1287.35, 1317.76, 1368.09, 1418.42, 1425.29, 1468.76, 1519.09, 1558.03, 1569.42, 1619.91, 1670.08, 1720.41, 1770.74, 1817.19, 1821.07, 1871.4, 1921.74, 1943.5, 1972.07, 2022.4, 2067.31, 2072.73, 2123.06, 2173.39, 2188.29, 2223.72, 2248.94, 2274.05, 2290.18, 2306.26, 2324.38, 2349.61, 2362.18, 2374.72, 2385.34, 2390.64, 2395.94, 2401.23, 2406.52, 2409.16, 2411.8, 2414.43, 2415.75, 2417.07, 2420.96, 2425.05, 2446.3, 2460.12, 2473.91, 2475.38, 2488.29, 2494.73, 2501.16, 2525.71, 2531.59, 2536.4, 2571.39, 2576.04, 2605.15, 2626.37, 2638.67, 2639.5, 2673.45, 2676.7, 2707.03, 2727.03, 2740.24, 2742.03, 2773.54, 2777.36, 2806.37, 2827.7, 2838.94, 2841.76, 2871.39, 2878.03, 2903.42, 2928.36, 2935.54, 2937.78, 2967.31, 2978.69, 2998.77, 3029.02, 3030.09, 3031.26, 3079.35, 3104.63, 3129.68, 3142.35, 3154.97, 3180.01, 3192.69, 3205.3, 3230.34, 3243.02, 3255.64, 3280.68, 3293.37, 3305.98, 3331.01, 3343.7, 3356.32, 3381.34, 3394.03, 3406.65, 3431.67, 3437.67, 3443.67, 3455.64, 3461.62, 3467.59, 3473.55, 3476.52, 3478.01, 3478.76, 3479.5, 3486.37, 3493.61, 3531.67, 3556.5, 3581.33, 3583.98, 3607.49, 3619.25, 3631.0, 3665.76, 3674.14, 3681.0, 3731.0, 3737.65, 3779.5, 3809.98, 3827.8, 3829.0, 3878.5, 3883.24, 3928.0, 3957.72, 3977.5, 3980.17, 4027.0, 4032.71, 4076.5, 4108.8, 4126.0, 4130.3, 4175.5, 4185.68, 4225.0, 4263.33, 4274.5, 4277.99, 4324.0, 4341.75, 4373.5, 4421.28, 4423.0, 4424.88, 4472.5, 4501.99, 4522.0, 4542.81, 4571.5, 4583.18, 4603.81, 4621.0, 4624.35, 4665.11, 4670.5, 4686.03, 4706.85, 4720.0, 4747.86, 4768.75, 4769.5, 4789.53, 4819.0, 4830.82, 4852.04, 4868.5, 4873.09, 4914.42, 4918.0, 4935.6, 4956.67, 4967.5, 4998.3, 5017.0, 5019.1, 5039.66, 5066.5, 5080.39, 5087.79, 5095.17, 5109.89, 5116.0, 5117.22, 5124.49, 5131.74, 5135.36, 5137.17, 5138.08, 5138.98, 5165.5, 5190.33, 5215.0, 5239.83, 5264.5, 5289.33, 5314.0, 5326.44, 5338.83, 5351.19, 5363.5, 5375.94, 5388.34, 5400.69, 5406.85, 5413.0, 5417.16, 5421.31, 5425.46, 5427.53, 5428.57, 5429.6, 5462.5, 5478.82, 5512.0, 5524.89, 5561.5, 5570.44, 5590.78, 5611.0, 5619.96, 5639.89, 5659.63, 5660.32, 5661.0, 5679.46, 5697.75, 5704.39, 5707.7, 5711.0, 5714.02, 5744.36, 5761.0, 5774.52, 5804.23, 5811.0, 5833.51, 5861.0, 5862.11, 5890.88, 5911.0, 5918.74, 5932.75, 5946.63, 5953.83, 5957.42, 5961.0, 5973.29, 5985.47, 5991.67, 5997.84, 6010.08, 6010.54, 6011.0, 6022.78, 6034.44, 6040.53, 6046.59, 6058.62, 6059.81, 6061.0, 6071.85, 6082.6, 6088.6, 6094.57, 6106.43, 6108.72, 6111.0, 6130.07, 6141.84, 6153.48, 6157.24, 6161.0, 6180.19, 6193.15, 6206.0, 6207.2, 6208.41, 6222.49, 6236.44, 6243.74, 6247.37, 6249.19, 6251.0, 6270.45, 6271.0, 6280.17, 6289.89, 6293.5, 6304.13, 6314.76, 6320.07, 6325.38, 6330.69, 6333.35, 6334.67, 6336.0, 6343.5, 6347.25, 6349.13, 6350.06, 6351.0, 6359.05, 6365.02, 6368.01, 6369.51, 6371.0]
            @test p[1].radius ≈ radius atol=0.1
            @test turning_depth(p[1]) == TauPy.RADIUS[p[1].model] - minimum(p[1].radius)
        end
    end

    @testset "Operators" begin
        let p = travel_time(100, 60, "ScS"), p′ = travel_time(100, 60, "ScS"),
                q = travel_time(100, 60, "PcP")
            # Vectors of Phases
            @test p == p
            @test p == p′
            @test p == deepcopy(p)
            @test (p == q) == false
            # Individual Phases
            @test p[1] == p[1]
            @test p[1] == p′[1]
            @test p[1] == deepcopy(p[1])
            @test (p[1] == q[1]) == false
        end
    end
end

@testset "PhaseGeog" begin
    @testset "Time" begin
        # Test against taup_time, and assume we're on a spherical Earth
        let p = travel_time(0, 0, 300, 30, 30, "PcS", model="prem")
            @test p isa Vector{<:PhaseGeog}
            @test length(p) == 1
            @test p[1].model == "prem"
            @test p[1].name == "PcS"
            @test p[1].evlon ≈ 0.0
            @test p[1].evlat ≈ 0.0
            @test p[1].stlon ≈ 30.0
            @test p[1].stlat ≈ 30.0
            @test p[1].depth ≈ 300.0
            @test p[1].delta ≈ 41.4096 atol=0.001
            @test p[1].time ≈ 780.67 atol=0.01
            @test p[1].dtdd ≈ 3.999 atol=0.001
            @test p[1].inc ≈ 6.6088 atol=0.001
            @test p[1].takeoff ≈ 19.1989 atol=0.01
            @test p[1].pierce ≈ Float64[]
            @test p[1].lon ≈ Float64[]
            @test p[1].lat ≈ Float64[]
            @test p[1].radius ≈ Float64[]
        end
    end

    @testset "Operators" begin
        let p = travel_time(0, 0, 300, 30, 30, "ScS"),
                p′ = travel_time(0, 0, 300, 30, 30, "ScS"),
                q = travel_time(0, 0, 300, 30, 30, "PcP")
            # Vectors of Phases
            @test p == p
            @test p == p′
            @test p == deepcopy(p)
            @test (p == q) == false
            # Individual Phases
            @test p[1] == p[1]
            @test p[1] == p′[1]
            @test p[1] == deepcopy(p[1])
            @test (p[1] == q[1]) == false
        end
    end

    @testset "Path" begin
        let elon = 0, elat = 0, depth = 100, slon = 45, slat = 45, phase = "P", model = "ak135"
            p = path(elon, elat, depth, slon, slat, "P", model=model)
            @test p isa Vector{<:PhaseGeog}
            @test length(p) == 1
            @test p[1].delta ≈ 60.0
            @test p[1].dtdd ≈ 6.835 atol=0.001
            @test p[1].evlon ≈ elon
            @test p[1].evlat ≈ elat
            @test p[1].stlon ≈ slon
            @test p[1].stlat ≈ slat
            @test p[1].inc ≈ 20.888 atol=0.001
            @test p[1].model == model
            @test p[1].takeoff ≈ 30.173 atol=0.001
            @test p[1].time ≈ 595.99 atol=0.01
            @test p[1].radius[1:5] ≈ [6271.00, 6270.44, 6251.00, 6249.18, 6247.37] atol=0.1
            @test length(p[1].lon) == length(p[1].lat) == length(p[1].radius)
        end
    end

    @testset "Cache" begin
        let Δ = rand(1:90), h = rand(0:700), phase = "P", model = "sp6"
            TauPy.clear_cache!()
            @test length(TauPy.RAY_CACHE) == length(TauPy.RAY_NCALLS) == 0
            p = travel_time(h, Δ, phase, model=model)
            @test first(values(TauPy.RAY_NCALLS)) == 1
            p′ = travel_time(h, Δ, phase, model=model)
            @test first(values(TauPy.RAY_NCALLS)) == 2
            @test p == p′
            @test travel_time(h, Δ, phase, model=model, cache=false) == p
            @test first(values(TauPy.RAY_NCALLS)) == 2
        end
    end
end

max_stars_repo_name: anowacki/TauPy


text:
using Slomo

model = JeansModel(Halos.NFWModel(), SersicModel(), ConstantBetaModel(0.5))
R = 10 .^ collect(-1:0.1:1)
betas = collect(-0.5:0.1:0.5)
params = [Dict(:beta => b) for b in betas]

s = sigma_los_parallel(model, R, params)


text:
using FFTW
using StatsBase

struct PoissonSolverFFT{DT <: Real} <: PoissonSolver{DT}
    nx::Int
    Δx::DT
    xgrid::Vector{DT}
    cells::Vector{DT}
    L::DT
    ρ::Vector{DT}
    ϕ::Vector{DT}

    function PoissonSolverFFT{DT}(nx::Int, L::DT) where {DT}
        Δx = L/(nx+1)
        xgrid = LinRange(Δx/2, L-Δx/2, nx)
        cells = LinRange(0, L, nx+1)
        new(nx, Δx, xgrid, cells, L, zeros(DT, nx), zeros(DT, nx))
    end
end

PoissonSolverFFT(nx::Int, L::DT) where {DT} = PoissonSolverFFT{DT}(nx, L)

Base.length(p::PoissonSolverFFT) = p.nx

function solve!(p::PoissonSolverFFT{DT}, x::AbstractVector{DT}, w::AbstractVector{DT} = one.(x) ./ length(x)) where {DT}
    h = fit(Histogram, mod.(x, p.L), p.cells)
    p.ρ .= h.weights ./ length(x)
    ρ̂ = rfft(p.ρ)
    k² = [(i-1)^2 for i in eachindex(ρ̂)]
    ϕ̂ = - ρ̂ ./ k²
    ϕ̂[1] = 0
    p.ϕ .= irfft(ϕ̂, length(p.ρ))
    return p
end

function get_indices(p::PoissonSolverFFT, x)
    y = mod(x, p.L)

    i1 = floor(Int, y / p.Δx)
    i2 = mod( ceil(Int, y / p.Δx) - 1, p.nx) + 1

    i1 == 0 && (i1 = p.nx)
    i2 == 0 && (i2 = p.nx)

    i1 == p.nx+1 && (i1 = 1)
    i2 == p.nx+1 && (i2 = 1)

    return (i1, i2)
end

function get_index(p::PoissonSolverFFT, x)
    i1, i2 = get_indices(p, x)
    return (abs(p.xgrid[i1] - x) ≤ abs(p.xgrid[i2] - x) ? i1 : i2)
end

function eval_density(p::PoissonSolverFFT{DT}, x::DT) where {DT}
    p.ρ[get_index(p, x)]
end

function eval_potential(p::PoissonSolverFFT{DT}, x::DT) where {DT}
    p.ϕ[get_index(p, x)]
end

function eval_field(p::PoissonSolverFFT{DT}, x::DT) where {DT}
    i1, i2 = get_indices(p, x)
    return - (p.ϕ[i2] - p.ϕ[i1]) / p.Δx
end


text:
############################################################
# Softmax Loss
############################################################
@defstruct SoftmaxLossLayer Layer (
  name :: AbstractString = "softmax-loss",
  (weight :: AbstractFloat = 1.0, weight >= 0),
  weights :: Array = [],
  normalize:: Symbol = :local,
  (dim :: Int = -2, dim != 0),
  (bottoms :: Vector{Symbol} = Symbol[], length(bottoms) == 2),
)
@characterize_layer(SoftmaxLossLayer,
  has_loss  => true,
  can_do_bp => true,
  is_sink   => true
)

type SoftmaxLossLayerState{T} <: LayerState
  layer    :: SoftmaxLossLayer
  loss     :: T

  softmax  :: SoftmaxLayerState
  logistic :: MultinomialLogisticLossLayerState
end

function setup(backend::Backend, layer::SoftmaxLossLayer, inputs::Vector{Blob}, diffs::Vector{Blob})
  data_type = eltype(inputs[1])

  softmax_layer = SoftmaxLayer(tops=Array{Symbol}(length(inputs)), bottoms=Array{Symbol}(length(inputs)), dim=layer.dim)
  softmax = setup(backend, softmax_layer, Blob[inputs[1]], Blob[])

  logistic_layer = MultinomialLogisticLossLayer(bottoms=Array{Symbol}(2),
      weights=layer.weights, normalize=layer.normalize, dim=layer.dim)
  logistic = setup(backend, logistic_layer, inputs, Blob[])

  state = SoftmaxLossLayerState(layer, convert(data_type, 0), softmax, logistic)
  return state
end
function shutdown(backend::Backend, state::SoftmaxLossLayerState)
  shutdown(backend, state.softmax)
  shutdown(backend, state.logistic)
end

function forward(backend::Backend, state::SoftmaxLossLayerState, inputs::Vector{Blob})
  forward(backend, state.softmax, Blob[inputs[1]])
  forward(backend, state.logistic, Blob[state.softmax.blobs[1], inputs[2]])
  state.loss = state.logistic.loss * state.layer.weight
end

function backward(backend::CPUBackend, state::SoftmaxLossLayerState, inputs::Vector{Blob}, diffs::Vector{Blob})
  diff = diffs[1]
  if isa(diff, CPUBlob)
    dims = size(diff)
    label_dim = [i == state.logistic.op_dim ? 1 : dims[i] for i = 1:length(dims)]
    label = reshape(inputs[2].data, label_dim...)

    idx_all = map(1:length(dims)) do i
      if i == state.logistic.op_dim
        map(x -> round(Int, x), label) .+ 1
      else
        dim = dims[i]
        reshape(1:dim, [j == i? dim : 1 for j = 1:length(dims)]...)
      end
    end

    if isa(state.logistic.weights_blob, NullBlob)
      copy!(diff, state.softmax.blobs[1])
    else
      copy!(diff, state.softmax.blobs[1].data .*
          broadcast_getindex(state.logistic.weights_blob.data, idx_all...))
    end

    diff_data = reshape(diff.data, dims)
    if isa(state.logistic.weights_blob, NullBlob)
      broadcast_setindex!(diff_data, broadcast_getindex(diff_data, idx_all...)-1, idx_all...)
    else
      broadcast_setindex!(diff_data, broadcast_getindex(diff_data, idx_all...) .-
          broadcast_getindex(state.logistic.weights_blob.data, idx_all...), idx_all...)
    end
    Vec.mul_scal!(diff.data, state.layer.weight * dims[state.logistic.op_dim]/prod(dims))
  end
end


text:
#!/usr/bin/env python
# coding: utf-8

# In[1]:

using Pkg
Pkg.add("Latexify")
using Latexify

Pkg.add("Plots")
using Plots
Pkg.add("PyPlot")
pyplot()

# In[56]:

# Linear Stommel 1948 solution from Giraldo and Warburton, 2008 "High-order Triangular DG Oceanic Shallow Water Model", Int. J. Numer. Meth. Fluids, 56:899-925
# see section 5.1.3 "Linear Stommel Model"

# Problem setup
γ=1.e-6                 # Linear drag (s-1)
L=1e6                   # Basin size in X and Y (m)
L_x=L
L_y=L
β=1.e-11                # Rate of change of Coriolis with distance  ( (s-1)/m )
f₀=1.e-4                # Coriolis ( s-1 )
τ=0.2                   # Wind stress ( N m-2 )
ρ=1000                  # Density ( kg m-3 )
H=1000                  # Fluid depth ( m )
grav=10                 # grav ( m^2 s-1 )
Δ=L/100                 # Nominal spacing

# Working variables
β_p=β*L_y/f₀
ϵ=γ/(L*β)

# In[57]:

ETMP=:( sqrt(1+(2*π*ϵ)^2)    )
R₊exp=:( ( -1+$ETMP )/(2*ϵ) )
R₋exp=:( ( -1-$ETMP )/(2*ϵ) )
R₊def=:( R₊=$R₊exp )
R₋def=:( R₋=$R₋exp )
display(latexify(R₊def))
display(latexify(R₋def))

# In[58]:

ETMP01=:( R₊*(exp(R₋)-1)   )
ETMP02=:( R₋*(1-exp(R₊))   )
ETMP03=:( exp(R₊)-exp(R₋)  )
Dexpr=:( ( $ETMP01+$ETMP02 ) / $ETMP03 )
Ddef=:( D = $Dexpr )
display(latexify(Ddef))

# In[59]:

ETMP04=:($ETMP01*exp(R₊*x₁))
ETMP05=:($ETMP02*exp(R₋*x₁))
R₂expr=:( 1/D*( ($ETMP04+$ETMP05 )/$ETMP03 ) )
R₂def=:( R₂ = $R₂expr )
display(latexify(R₂def))

# In[60]:

ETMP06=:((exp(R₋)-1)*exp(R₊*x₁))
ETMP07=:((1-exp(R₊))*exp(R₋*x₁))
R₁expr=:( π/D*((1.0.+(($ETMP06.+$ETMP07))/($ETMP03))) )
R₁def=:( R₁ = $R₁expr )
display(latexify(R₁def,env=:inline))

# In[61]:

u_scale_expr=:(τ*D/(ρ*H*γ*π))
u_scale_def=:(u_scale = $u_scale_expr )
display(latexify(u_scale_def,env=:inline))

# In[62]:

h_scale_expr=:(f₀*L*u_scale/grav)
h_scale_def=:( h_scale = $h_scale_expr)
display(latexify(h_scale_def,env=:inline))

# In[63]:

U_curly_expr=:(R₁*sin.(π*y₁))
U_curly_def=:( U_curly = $U_curly_expr )
display(latexify(U_curly_def,env=:inline))

# In[64]:

V_curly_expr=:(R₂*cos.(π*y₁))
V_curly_def=:( V_curly = $V_curly_expr )
display(latexify(V_curly_def,env=:inline))

# In[65]:

Hexpr01=:(-R₂/(π*f₀)*γ*sin(π*y₁))
Hexpr02=:( R₁/(π)*( ( cos(π*y₁)*(1.0+β_p*y₁) ) - β_p/π*sin(π*y₁) )  )
Hexpr=:($Hexpr01 + $Hexpr02)
Hexpr_def=:( H_curly = $Hexpr)

# In[66]:

function showdef(m)
 display(latexify(Hexpr_def,env=m))
 display(latexify(U_curly_def,env=m))
 display(latexify(V_curly_def,env=m))
 display(latexify(u_scale_def,env=m))
 display(latexify(h_scale_def,env=m))
 display(latexify(R₁def,env=m))
 display(latexify(R₂def,env=m))
 display(latexify(Ddef,env=m))
 display(latexify(R₊def,env=m))
 display(latexify(R₋def,env=m))

 display(Hexpr_def)
 display(U_curly_def)
 display(V_curly_def)
 display(u_scale_def)
 display(h_scale_def)
 display(R₁def)
 display(R₂def)
 display(Ddef)
 display(R₊def)
 display(R₋def)
end
showdef(:inline)

# In[67]:

R₋ = (-1 - sqrt(1 + (2 * π * ϵ) ^ 2)) / (2ϵ)
R₊ = (-1 + sqrt(1 + (2 * π * ϵ) ^ 2)) / (2ϵ)
D = (R₊ * (exp(R₋) - 1) + R₋ * (1 - exp(R₊))) / (exp(R₊) - exp(R₋))
R₂(x₁) = (1 / D) * (((R₊ * (exp(R₋) - 1)) * exp(R₊ * x₁) + (R₋ * (1 - exp(R₊))) * exp(R₋ * x₁)) / (exp(R₊) - exp(R₋)))
R₁(x₁) = (π / D) * (1.0 .+ ((exp(R₋) - 1) * exp(R₊ * x₁) .+ (1 - exp(R₊)) * exp(R₋ * x₁)) / (exp(R₊) - exp(R₋)))
u_scale = (τ * D) / (ρ * H * γ * π)
h_scale = (f₀ * L * u_scale) / grav
V_curly(x₁,y₁) = R₂(x₁) * cos.(π * y₁)
U_curly(x₁,y₁) = R₁(x₁) * sin.(π * y₁)
H_curly(x₁,y₁) = (-R₂(x₁) / (π * f₀)) * γ * sin(π * y₁) + (R₁(x₁) / π) * (cos(π * y₁) * (1.0 + β_p * y₁) - (β_p / π) * sin(π * y₁))

# In[68]:

# Temporary grid locations for plotting (should come from DG mesh function)
x=collect(range(0.,L_x,step=Δ))
x̃=x./L_x
y=collect(range(0.,L_y,step=Δ))
ỹ=y./L_y
nx=length( x̃ )
ny=length( ỹ )
Hval=zeros(nx,ny)
Uval=zeros(nx,ny)
Vval=zeros(nx,ny)
for j=1:ny
 for i=1:nx
  x₁=x̃[i]
  y₁=ỹ[j]
  Hval[i,j]=H_curly(x₁,y₁)
  Uval[i,j]=U_curly(x₁,y₁)
  Vval[i,j]=V_curly(x₁,y₁)
 end
end

# In[69]:

p1 = contour(y, x, h_scale.*Hval')
p2 = contour(y, x, u_scale.*Uval')
p3 = contour(y, x, u_scale.*Vval')
plot(p1,p2,p3)

# In[ ]:

max_stars_repo_name: christophernhill/climahacks

max_issues_repo_name: christophernhill/climahacks

max_forks_repo_name: christophernhill/climahacks


text:
# HSA dispatch packet offsets
_packet_names = fieldnames(HSA.KernelDispatchPacket)
_packet_offsets = fieldoffset.(HSA.KernelDispatchPacket, 1:length(_packet_names))

include(joinpath("gcn", "math.jl"))
include(joinpath("gcn", "indexing.jl"))
include(joinpath("gcn", "assertion.jl"))
include(joinpath("gcn", "synchronization.jl"))
include(joinpath("gcn", "memory_static.jl"))
include(joinpath("gcn", "hostcall.jl"))
include(joinpath("gcn", "output.jl"))
include(joinpath("gcn", "memory_dynamic.jl"))
include(joinpath("gcn", "execution_control.jl"))
include(joinpath("gcn", "atomics.jl"))


text:
# Cartesian product. See https://en.wikipedia.org/wiki/Cartesian_product

# nextTuple is a helper function which extracts all combinations of arrays by
# using the length of each array Yn and counting up from 0 at each index. For
# example, if we had an array X of lengths [2, 3, 3] nextTuple will adjust indices
# as follows: [0, 0, 0] -> [1, 0, 0] -> [0, 1, 0] -> [1, 1, 0] -> [0, 2, 0] ... [1, 2, 2]

function nextTuple!(indices::Array{Int,1}, sizes::Array{Int,1})::Bool
    for i in length(indices):-1:1
        if indices[i] < sizes[i]
            indices[i] += 1
            return true
        else
            indices[i] = 1
        end
    end
    false
end

# N fold Cartesian product
function cartesianProduct(aa::Array{Array{T, 1}, 1})::Array{Array{T, 1}, 1} where {T}
    result::Array{Array{T, 1}, 1} = []
    indices = ones(Int, length(aa))
    sizes = map(a -> length(a), aa)
    hasMore = length(indices) > 0
    while hasMore
        elt = [aa[i][indices[i]] for i in 1:length(indices)]
        push!(result, elt)
        hasMore = nextTuple!(indices, sizes)
    end

  return result
end


text:
export idx, cpx, exp_ikx, propagator
export gaussian, normal
export ramp, box, disc

function ramp(::Type{T}, dim::Int, dim_size::Int;
    offset=CtrFT, scale=ScaUnit) where {T}
    size = single_dim_size(dim,dim_size)
    offset = get_offset(size, offset)
    scale_n = get_scale(size, scale)
    f = ((x) -> T.(scale_n[dim] .* (x[dim] .- offset[dim])))
    IndexFunArray(T, f, size)
end

function ramp(dim::Int, dim_size::Int; offset=CtrFT, scale=ScaUnit)
    ramp(DEFAULT_T, dim, dim_size; offset=offset, scale=scale)
end

# values in the complex plane
function cpx(::Type{T}, size::NTuple{N, Int}; offset=CtrFT, scale=ScaUnit, dims=ntuple(+, N)) where {N,T}
    ig = idx(T, size, offset=offset, scale=scale, dims=dims).generator
    f(x) = complex(ig(x)...)
    return IndexFunArray(Complex{T}, f, size)
end
# values in the complex plane
function cpx(size::NTuple{N, Int}; offset=CtrFT, scale=ScaUnit, dims=ntuple(+, N)) where {N,T}
    cpx(DEFAULT_T, size, offset=offset, scale=scale, dims=dims)
end
function cpx(arr::AbstractArray{T, N}; offset=CtrFT, scale=ScaUnit, dims=ntuple(+, N)) where {N,T}
    cpx(T, size(arr), offset=offset, scale=scale, dims=dims)
end

function exp_ikx(::Type{T}, size::NTuple{N, Int}; shift_by=size.÷2, offset=CtrFT, scale=ScaFT, dims=ntuple(+, N),
                accumulator=sum, weight=1) where {N,T}
    myscale = apply_tuple_list((x,y)-> T.(-2pi .* x .* y), get_scale(size, scale), optional_mat_to_iter(shift_by))
    return exp_is(T, size, scale = myscale, offset=offset, dims = dims, accumulator=accumulator, weight=weight)
end

function exp_ikx(size::NTuple{N, Int}; shift_by=size.÷2, offset=CtrFT, scale=ScaFT, dims=ntuple(+, N), accumulator=sum, weight=1) where {N,T}
    myscale = apply_tuple_list((x,y)-> DEFAULT_T.(-2pi .* x .* y), get_scale(size, scale), optional_mat_to_iter(shift_by))
    return exp_is(complex(DEFAULT_T), size, scale = myscale, offset=offset, dims = dims, accumulator=accumulator, weight=weight)
end

function exp_ikx(arr::AbstractArray{T, N}; shift_by=size(arr).÷2, offset=CtrFT, scale=ScaFT, dims=ntuple(+, N), accumulator=sum, weight=1) where {N,T}
    myscale = apply_tuple_list((x,y)-> T.(-2pi .* x .* y), get_scale(size(arr), scale), optional_mat_to_iter(shift_by))
    return exp_is(complex(eltype(arr)),size(arr), scale = myscale,  offset=offset, dims = dims, accumulator=accumulator)
end

function propagator(::Type{T}, size::NTuple{N, Int}; Δz=1.0, shift_by=0, k_max=0.5, offset=CtrFT, scale=ScaFT, dims=ntuple(+, N),
    accumulator=sum, weight=1) where {N,T}
myscale = apply_tuple_list((x,y)-> T.(x ./ y), get_scale(size, scale), optional_mat_to_iter(k_max))
myscale2 = apply_tuple_list((x,y)-> T.(-2pi .* x .* y), get_scale(size, scale), optional_mat_to_iter(shift_by))
return cispi.((2 .*Δz) .*
    phase_kz(T, size, scale = myscale, offset=offset, dims = dims, accumulator=accumulator, weight=weight) .+
    2 .* phase_kxy(T, size, scale = myscale2, offset=offset, dims = dims, accumulator=accumulator, weight=weight))
end

function propagator(size::NTuple{N, Int}; Δz=1.0,  k_max=0.5, shift_by=0, offset=CtrFT, scale=ScaFT, dims=ntuple(+, N), accumulator=sum, weight=1) where {N,T}
myscale = apply_tuple_list((x,y)-> DEFAULT_T.(x ./ y), get_scale(size, scale), optional_mat_to_iter(k_max))
myscale2 = apply_tuple_list((x,y)-> DEFAULT_T.(-2pi .* x .* y), get_scale(size, scale), optional_mat_to_iter(shift_by))
return cispi.((2 .*Δz) .*
        phase_kz(DEFAULT_T, size, scale = myscale, offset=offset, dims = dims, accumulator=accumulator, weight=weight) .+
        2 .* phase_kxy(DEFAULT_T, size, scale = myscale2, offset=offset, dims = dims, accumulator=accumulator, weight=weight))
    end

function propagator(arr::AbstractArray{T, N}; Δz=1.0, k_max=0.5, shift_by=0, offset=CtrFT, scale=ScaFT, dims=ntuple(+, N), accumulator=sum, weight=1) where {N,T}
myscale = apply_tuple_list((x,y)-> T.(x ./ y), get_scale(size(arr), scale), optional_mat_to_iter(k_max))
myscale2 = apply_tuple_list((x,y)-> T.(-2pi .* x .* y), get_scale(size(arr), scale), optional_mat_to_iter(shift_by))
return cispi.((2 .*Δz) .*
        phase_kz(eltype(arr),size(arr), scale = myscale,  offset=offset, dims = dims, accumulator=accumulator) .+
        2 .* phase_kxy(eltype(arr), size, scale = myscale2, offset=offset, dims = dims, accumulator=accumulator, weight=weight))
end

# These versions of Gaussians interpret sigma as the standard-deviation along each axis separately
function gaussian(::Type{T}, size::NTuple{N, Int}; sigma=1.0, offset=CtrFT, scale=ScaUnit, dims=ntuple(+, N), accumulator=sum, weight=1) where {N,T}
    myscale = apply_tuple_list((x,y)-> T.(x ./( 2 .* y .*y)), get_scale(size, scale), optional_mat_to_iter(sigma))
    return exp_sqr(T, size, scale = myscale, offset=offset, dims = dims, accumulator=accumulator, weight=weight)
end

function gaussian(size::NTuple{N, Int}; sigma=1.0, offset=CtrFT, scale=ScaUnit, dims=ntuple(+, N), accumulator=sum, weight=1) where {N,T}
    myscale = apply_tuple_list((x,y)-> DEFAULT_T.(x ./( 2 .* y .*y)), get_scale(size, scale), optional_mat_to_iter(sigma))
    return exp_sqr(DEFAULT_T, size, scale = myscale, offset=offset, dims = dims, accumulator=accumulator, weight=weight)
end

function gaussian(arr::AbstractArray{T, N}; sigma=1.0, offset=CtrFT, scale=ScaUnit, dims=ntuple(+, N), accumulator=sum, weight=1) where {N,T}
    myscale = apply_tuple_list((x,y)-> T.(x ./( 2 .*y.*y)), get_scale(size(arr), scale), optional_mat_to_iter(sigma))
    return exp_sqr(arr, scale = myscale,  offset=offset, dims = dims, accumulator=accumulator, weight=weight)
end

function normal(::Type{T}, size::NTuple{N, Int}; sigma=1.0, offset=CtrFT, scale=ScaUnit, dims=ntuple(+, N), accumulator=sum, weight=1) where {N,T}
    myscale = apply_tuple_list((x,y)-> T.(x ./( 2 .*y.*y)), get_scale(size, scale), optional_mat_to_iter(sigma))
    return exp_sqr_norm(T, size, scale = myscale, offset=offset, dims = dims, accumulator=accumulator, weight=weight)
end

function normal(sz::NTuple{N, Int}; sigma=1.0, offset=CtrFT, scale=ScaUnit, dims=ntuple(+, N), accumulator=sum, weight=1) where {N}
    myscale = apply_tuple_list((x,y)-> DEFAULT_T.(x ./( 2 .*y.*y)), get_scale(sz, scale), optional_mat_to_iter(sigma))
    return exp_sqr_norm(DEFAULT_T, sz, scale = myscale, offset=offset, dims = dims, accumulator=accumulator, weight=weight)
end

function normal(arr::AbstractArray{T, N}; sigma=1.0, offset=CtrFT, scale=ScaUnit, dims=ntuple(+, N), accumulator=sum, weight=1) where {N,T}
    myscale = apply_tuple_list((x,y)-> T.(x ./( 2 .*y.*y)), get_scale(size(arr), scale), optional_mat_to_iter(sigma))
    return exp_sqr_norm(arr, scale = myscale,  offset=offset, dims = dims, accumulator=accumulator, weight=weight)
end

function box(::Type{T}, sz::NTuple{N, Int}, boxsize; offset=CtrFT, scale=ScaUnit, dims=ntuple(+, N)) where {N,T}
    return convert.(T,idx_max(sz,scale=get_scale(sz, scale) ./ boxsize, offset=offset, dims=dims) .< 0.5)
end

function box(sz::NTuple{N, Int}, boxsize; offset=CtrFT, scale=ScaUnit, dims=ntuple(+, N)) where {N}
    box(DEFAULT_T, sz, boxsize; offset=offset, scale=scale, dims=dims)
end

function box(::Type{TR}, arr::AbstractArray{T, N}, boxsize; offset=CtrFT, scale=ScaUnit, dims=ntuple(+, N)) where {N,T,TR}
    return convert.(T, idx_max(size(arr),scale=get_scale(size(arr), scale) ./ boxsize, offset=offset, dims=dims) .< 0.5)
end

function box(arr::AbstractArray{T, N}, boxsize; offset=CtrFT, scale=ScaUnit, dims=ntuple(+, N)) where {N,T}
    box(eltype(arr), arr, boxsize; offset=offset, scale=scale, dims=dims)
end

function disc(::Type{T}, sz::NTuple{N, Int}, disc_radius; offset=CtrFT, scale=ScaUnit, dims=ntuple(+, N)) where {N,T}
    return convert.(T,rr2(sz,scale=get_scale(sz, scale) ./ disc_radius, offset=offset, dims=dims) .<= 1.0)
end

function disc(sz::NTuple{N, Int}, disc_radius; offset=CtrFT, scale=ScaUnit, dims=ntuple(+, N)) where {N}
    disc(DEFAULT_T, sz, disc_radius; offset=offset, scale=scale, dims=dims)
end

function disc(::Type{TR}, arr::AbstractArray{T, N}, disc_radius; offset=CtrFT, scale=ScaUnit, dims=ntuple(+, N)) where {N,T,TR}
    return convert.(T, rr2(size(arr),scale=get_scale(size(arr), scale) ./ disc_radius, offset=offset, dims=dims) .<= 1.0)
end

function disc(arr::AbstractArray{T, N}, disc_radius; offset=CtrFT, scale=ScaUnit, dims=ntuple(+, N)) where {N,T}
    disc(eltype(arr), arr, disc_radius; offset=offset, scale=scale, dims=dims)
end


text:
module testPortfolio

using Base.Test
using DataFrames
using DataArrays
using Datetime
using TimeData
include(string(Pkg.dir("AssetMgmt"), "/src/AssetMgmt.jl"))

println("\n Running portfolio tests\n")

simVals = rand(1, 8)
simVals = AssetMgmt.makeWeights(simVals)

df = DataFrame(simVals)
pf = AssetMgmt.Portfolio(df)
pf = AssetMgmt.Portfolio(simVals)

simVals = rand(4, 3)
simVals = AssetMgmt.makeWeights(simVals)
pf = AssetMgmt.Portfolio(simVals)

simVals = rand(1, 8)
@test_throws AssetMgmt.Portfolio(simVals)

simVals = rand(1, 8)
simVals = AssetMgmt.makeWeights(simVals)
df = DataFrame(simVals)
df[1, 2] = NA
@test_throws AssetMgmt.Portfolio(df)

end


text:
using EmpiricalCDFs

# Create an empirical cdf with 10^6 normally distributed samples
cdf = EmpiricalCDF(randn(10^6))

# add another sample
push!(cdf, randn())

# add 10^6 more samples
append!(cdf, randn(10^6))

# sort the cdf before using it
sort!(cdf)

# Pick a random sample from the distribution approximated by `cdf`
rand(cdf)

# The normal distribution is symmetric, so this should be about 0.5
cdf(0)

# This should be approximately zero
getinverse(cdf,.5)


text:
mutable struct IncrementBinner{T}
    keep::Int64         # number of values per block
    compression::Int64  # current N-mean
    stage::Int64        # stage * keep == current max length

    # for local averages
    sum::T
    count::Int64

    output::Vector{T}
end

# Overload some basic Base functions
Base.eltype(B::IncrementBinner{T}) where {T} = T
Base.length(B::IncrementBinner) = length(B.output)
Base.ndims(B::IncrementBinner{T}) where {T} = ndims(eltype(B))
Base.isempty(B::IncrementBinner) = length(B) == 0

Base.show(io::IO, B::IncrementBinner{T}) where {T} = print(io, "IncrementBinner{$T}()")

"""
    IncrementBinner([::Type{T}; blocksize = 64])
    IncrementBinner(zero_element::T[; blocksize = 64])

Creates an `IncrementBinner` from a `zero_element` numeric type `T`.

Values pushed to an `IncrementBinner` are averaged in stages. For the first
`blocksize` values pushed no averaging happens. After that `2blocksize` elements
are averaged 2 at a time, then `4blocksize` element 4 at a time, etc. This means
values pushed become progressively more compressed and smoothed.
"""
IncrementBinner(::Type{T} = Float64; kw...) where {T} = IncrementBinner(zero(T); kw...)
function IncrementBinner(x::T; blocksize = 64) where {T <: Union{Number, AbstractArray}}
    # check keyword args
    blocksize <= 0 && throw(ArgumentError("`blocksize` must be finite and positive."))

    # got_timeseries = didn't receive a zero && is a vector
    got_timeseries = count(!iszero, x) > 0 && ndims(T) == 1

    if got_timeseries
        # x = time series
        S = _sum_type_heuristic(eltype(T), ndims(x[1])) # from log
        el = zero(x[1])
    else
        # x = zero_element
        S = _sum_type_heuristic(T, ndims(x))
        el = x
    end

    B = IncrementBinner(blocksize, 1, 1, el, 0, S[])

    got_timeseries && append!(B, x)

    return B
end

function Base.push!(b::IncrementBinner{T}, x::T) where {T}
    if b.count >= b.compression - 1
        push!(b.output, (b.sum + x) / b.compression)
        b.count = 0
    elseif b.count == 0
        b.sum = x
        b.count = 1
    else
        b.sum += x
        b.count += 1
    end
    if length(b.output) == b.stage * b.keep
        b.compression *= 2
        b.stage += 1
    end
    nothing
end

function Base.push!(b::IncrementBinner{T}, x::T) where {T <: AbstractArray}
    if b.count >= b.compression - 1
        push!(b.output, (b.sum .+ x) / b.compression)
        b.count = 0
    elseif b.count == 0
        b.sum .= x
        b.count = 1
    else
        b.sum .+= x
        b.count += 1
    end
    if length(b.output) == b.stage * b.keep
        b.compression *= 2
        b.stage += 1
    end
    nothing
end

function Base.append!(b::IncrementBinner{T}, v::Vector{T}) where {T}
    for x in v
        push!(b, x)
    end
    nothing
end

Base.values(B::IncrementBinner) = B.output

function indices(B::IncrementBinner)
    out = zeros(length(B.output))
    step = 1
    out[1] = 1
    for i in 2:length(B.output)
        if i % B.keep == 1
            out[i] = out[i-1] + 1.5step
            step *= 2
        else
            out[i] = out[i-1] + step
        end
    end
    out
end

max_stars_repo_name: crstnbr/BinningAnalysis


text:
using Parameters, StatsBase, NamedTupleTools

import Distributions: logpdf, loglikelihood

struct Semantic{T1,T2} <: ContinuousUnivariateDistribution
    blc::T1
    parms::T2
end

loglikelihood(d::Semantic, data::Array{<:NamedTuple,1}) = logpdf(d, data)

Semantic(;blc, parms) = Semantic(blc, parms)

function populate_memory(act=0.0)
    chunks = [
        Chunk(object=:shark, attribute=:dangerous, value=:True, act=act),
        Chunk(object=:shark, attribute=:locomotion, value=:swimming, act=act),
        Chunk(object=:shark, attribute=:category, value=:fish, act=act),
        Chunk(object=:salmon, attribute=:edible, value=:True, act=act),
        Chunk(object=:salmon, attribute=:locomotion, value=:swimming, act=act),
        Chunk(object=:salmon, attribute=:category, value=:fish, act=act),
        Chunk(object=:fish, attribute=:breath, value=:gills, act=act),
        Chunk(object=:fish, attribute=:locomotion, value=:swimming, act=act),
        Chunk(object=:fish, attribute=:category, value=:animal, act=act),
        Chunk(object=:animal, attribute=:moves, value=:True, act=act),
        Chunk(object=:animal, attribute=:skin, value=:True, act=act),
        Chunk(object=:canary, attribute=:color, value=:yellow, act=act),
        Chunk(object=:canary, attribute=:sings, value=:True, act=act),
        Chunk(object=:canary, attribute=:category, value=:bird, act=act),
        Chunk(object=:ostritch, attribute=:flies, value=:False, act=act),
        Chunk(object=:ostritch, attribute=:height, value=:tall, act=act),
        Chunk(object=:ostritch, attribute=:category, value=:bird, act=act),
        Chunk(object=:bird, attribute=:wings, value=:True, act=act),
        Chunk(object=:bird, attribute=:locomotion, value=:flying, act=act),
        Chunk(object=:bird, attribute=:category, value=:animal, act=act),
    ]
    return chunks
end

function get_stimuli()
    stim = NamedTuple[]
    push!(stim, (object = :canary, category = :bird, ans = :yes))
    push!(stim, (object = :canary, category = :animal, ans = :yes))
    return vcat(stim...)
end

function logpdf(d::Semantic, data::Array{<:NamedTuple,1})
    LL = computeLL(d.parms, data; blc=d.blc)
    return LL
end

function simulate(parms, stimulus, n_reps; blc)
    # populate declarative memory
    chunks = populate_memory()
    # generate declarative memory object
    memory = Declarative(;memory=chunks)
    # generate ACTR object
    actr = ACTR(;declarative=memory, parms..., blc)
    # the number of correct responses
    k = 0
    # count the number of correct answers, k
    for rep in 1:n_reps
        k += simulate_trial(actr, stimulus)
    end
    # return data that constains stimulus information, number of trials,
    # and correct answers
    return (stimulus..., N = n_reps, k = k)
end

function simulate_trial(actr, stimulus)
    retrieving = true
    # create memory probe or retrieval request
    probe = stimulus
    chunks = actr.declarative.memory
    # k = 1 if answer is "yes", 0 otherwise
    k = 0
    while retrieving
        # generate retrieval probabilities
        p,_ = retrieval_probs(actr; object=probe.object, attribute=:category)
        # sample a chunk index proportional to retrieval probabilities
        idx = sample(1:length(p), weights(p))
        # Last element corresponds to a retrieval failure
        # stop retrieval processes
        if idx == length(p)
            retrieving = false
        # retrieved chunk matches retrieval request, stop retrieving
        # and set k = 1 for "yes" response
        elseif direct_verify(chunks[idx], probe)
            retrieving = false
            k += 1
        # perform another retrieval with category chaining
        # modify the retrieval request based on the retrieved chunk
        elseif chain_category(chunks[idx], probe)
            probe = delete(probe, :object)
            probe = (object = chunks[idx].slots.value, probe...)
        # no chunks match, stop retrieving and respond "no" with k = 0
        else
            retrieving = false
        end
    end
    return k
end

function simulate_exact(parms, stimulus, n_reps; blc)
    chunks = populate_memory()
    memory = Declarative(;memory=chunks)
    actr = ACTR(;declarative=memory, parms..., blc)
    k = 0
    for rep in 1:n_reps
        retrieving = true
        probe = stimulus
        while retrieving
            chunk = retrieve(actr; object=probe.object, attribute=:category)
            if isempty(chunk)
                retrieving = false
            elseif direct_verify(chunk[1], probe)
                retrieving = false
                k += 1
            elseif chain_category(chunk[1], probe)
                probe = delete(probe, :object)
                probe = (object = chunk[1].slots.value, probe...)
            else
                retrieving = false
            end
        end
    end
    return (stimulus..., N = n_reps, k = k)
end

"""
Answer yes via direct verification if retrieved chunk matches
probe on the object slot, the attribute slot equals category and the
value slot matches the value of the probe's category slot
"""
function direct_verify(chunk, probe)
    return match(chunk, object=probe.object,
        value=probe.category, attribute=:category)
end

"""
Chain category if retrieved chunk matches
probe on the object slot, the attribute slot equals category and the
value slot does not match the value of the probe's category slot
"""
function chain_category(chunk, probe)
    return match(chunk, ==, !=, ==, object=probe.object,
        value=probe.category, attribute=:category)
end

function initial_state(blc)
    s0 = zeros(typeof(blc), 5)
    s0[1] = 1
    return s0
end

function computeLL(parms, data; blc)
    act = zero(typeof(blc))
    # populate declarative memory
    chunks = populate_memory(act)
    # create declarative memory object
    memory = Declarative(;memory=chunks)
    # create act-r object
    actr = ACTR(;declarative=memory, parms..., blc)
    # create initial state vector
    s0 = initial_state(blc)
    LL = 0.0
    for d in data
        # create transition matrix
        tmat = transition_matrix(actr, d, blc)
        # compute probability of "yes"
        LL += probability_yes(tmat, s0, d)
    end
    return LL
end

function probability_yes(tmat, s0, d)
    z = s0' * tmat^3; θ = z[4]
    # sometimes θ is nan because of exponentiation of activation
    return isnan(θ) ? (return -Inf) : logpdf(Binomial(d.N, θ), d.k)
end

"""
populatates transition matrix consisting of 5 states:
* `s1`: initial retrieval
* `s2 `: chain category 1
* `s3`: chain category 2
* `s4`: respond yes
* `s5`: respond no
"""
function transition_matrix(actr, stim, blc)
    chunks = actr.declarative.memory
    Nc = length(chunks) + 1
    probe::typeof(stim) = stim
    probe = stim
    N = 5
    # populate transition matrix
    tmat = zeros(typeof(blc), N, N)
    # compute retrieval probabilities, p
    p,_ = retrieval_probs(actr; object=get_object(probe), attribute=:category)
    # find indices of chunks associated with direct verification, category chaining and mismatching conditions
    direct_indices = find_indices(actr, object=get_object(probe), value=get_category(probe))
    chain_indices = find_indices(actr, ==, !=, ==, object=get_object(probe), value=get_category(probe), attribute=:category)
    mismatch_indices = setdiff(1:Nc, direct_indices, chain_indices)
    # use indices to compute probability of category chain, direct verification (yes), and mismatch (no)
    tmat[1,2] = sum(p[chain_indices])
    tmat[1,4] = sum(p[direct_indices])
    tmat[1,5] = sum(p[mismatch_indices])
    # attempt to extract chunk associated with category chaining
    chain_chunk = get_chunks(actr,==,!=,==, object = get_object(probe),
    value = get_category(probe), attribute=:category)
    cnt = 1
    # continue the process above as long as category chaining can be performed.
    while !isempty(chain_chunk)
        cnt += 1
        probe = (object = get_chunk_value(chain_chunk[1]), delete(probe, :object)...)
        p,_ = retrieval_probs(actr; object=get_object(probe), attribute=:category)
        direct_indices = find_indices(actr, object=get_object(probe), value=get_category(probe))
        chain_indices = find_indices(actr, ==, !=, ==, object=get_object(probe), value=get_category(probe), attribute=:category)
        mismatch_indices = setdiff(1:Nc, direct_indices, chain_indices)
        tmat[cnt,2] = sum(p[chain_indices])
        tmat[cnt,4] = sum(p[direct_indices])
        tmat[cnt,5] = sum(p[mismatch_indices])
        chain_chunk = get_chunks(actr,==,!=,==, object = get_object(probe),
        value = get_category(probe), attribute=:category)
    end
    # set self-transitions to 1 if row i sums to 0.0
    map(i -> sum(tmat[i,:]) == 0.0 ? (tmat[i,i] = 1.0) : nothing, 1:size(tmat, 2))
    return tmat
end

"""
Type-stable accessors
"""
get_object(x) = x.object
get_category(x) = x.category
get_chunk_value(x) = x.slots.value

function hit_rate(parms, stimulus, n_reps; blc)
    data = simulate(parms, stimulus, n_reps; blc)
    return data.k / data.N
end

function parse_lisp_data(data)
    new_data = NamedTuple[]
    for row in eachrow(data)
        temp = (object = Symbol(row[1]),category = Symbol(row[2]),
        ans = Symbol(row[3]),N = row[4],k = row[5])
        push!(new_data, temp)
    end
    return vcat(new_data)
end


text:
module PredictableComponents

using CoupledFields
export InputSpace, ModelObj
 export AIF, CVfn, KCCA, gKYA, predict
 export G, parzen

"""
    AIF(G::Matrix, X::Vector, Y::Vector, τ::Int) \n
Auto-Information Function

"""
function AIF(G::Matrix{Float64}, X::T, Y::T, τ::Int) where T<:Vector{Float64}
    n = size(G,1)-τ
    Gx = G[1:n,1:n]
    Ya = Y[(1+τ):end,:]

# Gx needs to be centered so that cor(Uₓ) = Iₚ
    N = fill(1/n, n, n)
    Gx += -N*Gx + Gx*-N + N*Gx*N
    ev, Ux = eig(Symmetric(Gx))
    l = cumsum(flipdim(ev,1))/sum(ev)
    d = sum(l.≤0.99)
    j = n:-1:(n-d+1)
    U = Ux[:,j]
    yhat = U * (U'U \ U'Ya)
    imin = indmin(X[1:n])
    imax = indmax(X[1:n])
    return sign(yhat[imax]-yhat[imin])*cor(Ya, yhat)[1]
end

"""
    CVfn(parm::Matrix, X::Matrix, Y::Matrix, modelfn::Function)  \n
Cross-validation function
"""
function CVfn(parm::T, X::T, Y::T, modelfn::Function) where T<:Matrix{Float64}

    # CV metric parameters?
    # par[3] = τ lags

    # Model parameters
    # par[1] = Kernel parameter
    # par[2] = no. of PCs

    local n::Int64
    local p::Int64
    local τ::Int64
    nrg = size(parm, 1)
    n, p = size(X)
    # Setup
    ta = 1:floor(Int, n/2)
    tb = setdiff(1:n,ta)

    xm1 = mean(X[ta,:],1);     ym1 = mean(Y[ta,:],1)

    X1train = X[ta,:].-xm1;  Y1train = Y[ta,:].-ym1
    X1test = X[tb,:].-xm1;  Y1test = Y[tb,:].-ym1

    lags = [1:floor(Int, n/8);]

    Rsq = zeros(nrg)
    for i in 1:nrg
        percent = round(100*i/nrg)
        @printf " %03d%%" percent
        par = parm[i,:]
        τ = par[3]
        model1 = modelfn(par, X1train, Y1train, lags)
        W1 = model1.W[:,1:1]; A1 = model1.A[:,1:1]
        Rcv = X1test[:, 1:Int(par[2])]*W1
        Tcv = Y1test[:, 1:Int(par[2])]*A1
        Gr = G(Rcv)
        Rsq[i] = KCCA(Gr, Tcv, Int(par[3]))[1]
        println([parm[i:i,:] Rsq[i] ])
    end
    parm = [parm Rsq]
#    println(parm)
    imax = indmax(parm[:,end])
    return parm[imax,:]
end

"""
    G(X::Matrix) \n
Gram matrix of `X`
"""
function G(X::T) where T<:Matrix{Float64}
    kpars = GaussianKP(X)
    G = CoupledFields.Kf([1.0], X, kpars)
    return G
end

"""
    gKPrC(par::Vector, X::Matrix, Y::Matrix, lags::Vector{Int})  \n
gradient Kernel Predictable Components Analysis \n
The vector `par` has 3 elements: `par[1]=σ`, `par[2]=` the number of components to keep, and \n
`par[3]=τ` the lead time to calculate `A` the projection vector of `Y`.
"""
function gKPrC(par::Vector{Float64}, X::T, Y::T, lags::Vector{Int}) where T<:Matrix{Float64}
    local n::Int64
    local p::Int64
#    local lagmax::Int64
    Zx = X[:,1:Int(par[2])]
    Zy = Y[:,1:Int(par[2])]
    n, p = size(Zx)
    kpars = GaussianKP(Zx)
    ∇K = CoupledFields.∇Kf(par, Zx, kpars)
    Gx = CoupledFields.Kf(par, Zx, kpars)
    M = zeros(p,p)
    for τ in lags
        nr = n-τ
        H = Gx[1:nr,1:nr] \ Zy[(τ+1):n,:]
        F = H*H'
        Mi = zeros(p,p)
        for i in 1:nr
            dK = ∇K[i,1:nr,:]
            Mi += dK'*F*dK
        end
       M += Mi./nr
    end
    ev = eigfact(Symmetric(M))
    W = flipdim(ev.vectors, 2)
    Rx = Zx * W
    Gr = G(Rx[:,1:1])
    sv, A = KCCA(Gr, Zy, Int(par[3]))
    Ry = Zy * A
    return ModelObj(W, Rx, A, Ry, flipdim(ev.values,1), par, "gKPrC")
end

"""
    KCCA(G::Matrix, Y::Matrix, τ::Int)  \n
Kernel Canonical Correlation Analysis  \n
Calculates the projection vector of `Y`
"""
function KCCA(G::T, Y::T, τ::Int) where T<:Matrix{Float64}
    n = size(G,1)-τ
    Gx = G[1:n,1:n]
    Ya = Y[(1+τ):end,:]

# Gx needs to be centered so that cor(Uₓ) = Iₚ
    N = fill(1/n, n, n)
    Gx += -N*Gx + Gx*-N + N*Gx*N
    ev, Ux = eig(Symmetric(Gx))
    l = cumsum(flipdim(ev,1))/sum(ev)
    d = sum(l.≤0.99)
    j = n:-1:(n-d+1)
    M = cor(Ux[:,j],Ya)
    U,_,V = svd(M)
    sv = cor(Ux[:,j]*U[:,1], Ya*V[:,1])
    return sv^2, V
end

"""
    gKYA(model::ModelObj, Y::Matrix, dcca::Int) \n

"""
function gKYA(model::ModelObj, Y::Matrix{Float64}, dcca::Int)
    par = model.pars
    Zy = Y[:,1:Int(par[2])]
    A = zeros(Int(par[2]) ,dcca)
    for j in 1:dcca
        Gr = G(model.R[:,j:j])
        A[:,j] = KCCA(Gr, Zy, 0)[2][:,1]
    end
    T = Zy * A
    return ModelObj(model.W, model.R, A, T, model.evals, par, "gKPrC")
end

# kcr
"""
    predict(X::Vector, Y::Vector, τ::Int) \n

"""
function predict(X::T, Y::T, τ::Int) where T<:Vector{Float64}
    Gx = G(X[1:(end-τ),1:1])
    n = size(Gx,1)
    Ya = Y[(1+τ):end,:]

    # Gx needs to be centered so that cor(Uₓ) = Iₚ
    N = fill(1/n, n, n)
    Gx += -N*Gx + Gx*-N + N*Gx*N
    ev, Ux = eig(Symmetric(Gx))
    l = cumsum(flipdim(ev,1))/sum(ev)
    d = sum(l.≤0.99)
    j = n:-1:(n-d+1)
    U = Ux[:,j]
    yhat = U * (U'U \ U'Ya)
    return yhat[:]
end

"""
    parzen(tau::Vector{Int}, M::Int)
Parzen filter.
"""
function parzen(tau::Vector{Int}, M::Int)
    τ = abs.(tau)
    i = τ .>= M/2
    z =   (i).*(2(1-τ/M).^3) +  (.!i).*(1-6(τ/M).^2 + 6(τ/M).^3)
    z *= 1.0/sum(z)
    return z
end

end # module


text:
#!/usr/bin/env julia

# Tests against astropy.

using SkyCoords
using Base.Test

datapath = joinpath(dirname(@__FILE__), "data")
TOL = 0.0001  # tolerance in arcseconds

# Angular separation between two points (angles in radians)
#
#   The angular separation is calculated using the Vincenty formula [1]_,
#   which is slighly more complex and computationally expensive than
#   some alternatives, but is stable at at all distances, including the
#   poles and antipodes.
#
#   [1] http://en.wikipedia.org/wiki/Great-circle_distance
function angsep(lon1, lat1, lon2, lat2)

    sdlon = sin(lon2 - lon1)
    cdlon = cos(lon2 - lon1)
    slat1 = sin(lat1)
    slat2 = sin(lat2)
    clat1 = cos(lat1)
    clat2 = cos(lat2)

    num1 = clat2 * sdlon
    num2 = clat1 * slat2 - slat1 * clat2 * cdlon
    denom = slat1 * slat2 + clat1 * clat2 * cdlon

    atan2(sqrt(num1*num1 + num2*num2), denom)
end

lon(c::GalCoords) = c.l
lat(c::GalCoords) = c.b
lon(c::AbstractSkyCoords) = c.ra
lat(c::AbstractSkyCoords) = c.dec

function angsep{T<:AbstractSkyCoords}(c1::T, c2::T)
    angsep(lon(c1), lat(c1), lon(c2), lat(c2))
end

function angsep{T<:AbstractSkyCoords}(c1::Array{T}, c2::Array{T})
    size(c1) == size(c2) || error("size mismatch")
    result = similar(c1, Float64)
    for i=1:length(c1)
        result[i] = angsep(c1[i], c2[i])
    end
    result
end

rad2arcsec(r) = 3600.*rad2deg(r)

# input coordinates
fname = joinpath(datapath, "input_coords.csv")
indata, inhdr = readcsv(fname; header=true)

for (insys, T) in (("icrs", ICRSCoords), ("fk5j2000", FK5Coords{2000}),
                   ("fk5j1975", FK5Coords{1975}), ("gal", GalCoords))

    c_in = T[T(indata[i, 1], indata[i, 2]) for i=1:size(indata,1)]

    for (outsys, S) in (("icrs", ICRSCoords), ("fk5j2000", FK5Coords{2000}),
                        ("fk5j1975", FK5Coords{1975}), ("gal", GalCoords))
        (outsys == insys) && continue
        c_out = convert(Vector{S}, c_in)

        # Read in reference answers.
        fname = joinpath(datapath, "$(insys)_to_$(outsys).csv")
        refdata, hdr = readcsv(fname; header=true)
        c_ref = S[S(refdata[i, 1], refdata[i, 2]) for i=1:size(refdata,1)]

        # compare
        sep = angsep(c_out, c_ref)
        maxsep = rad2arcsec(maximum(sep))
        meansep = rad2arcsec(mean(sep))
        minsep = rad2arcsec(minimum(sep))
        @printf "%8s --> %8s : max=%6.4f\"  mean=%6.4f\"   min=%6.4f\"\n" insys outsys maxsep meansep minsep
        @test maxsep < TOL
    end
end

println()
println("All tests passed.")


text:
using Test
using TestSetExtensions
using LinearAlgebra
using Qaintessent

##==----------------------------------------------------------------------------------------------------------------------

# adapted from https://github.com/FluxML/Zygote.jl/blob/master/test/gradcheck.jl
function ngradient(f, xs::AbstractArray...)
    grads = zero.(xs)
    for (x, Δ) in zip(xs, grads), i in 1:length(x)
        δ = sqrt(eps())
        tmp = x[i]
        x[i] = tmp - δ/2
        y1 = f(xs...)
        x[i] = tmp + δ/2
        y2 = f(xs...)
        x[i] = tmp
        Δ[i] = (y2-y1)/δ
        if eltype(x) <: Complex
            # derivative with respect to imaginary part
            x[i] = tmp - im*δ/2
            y1 = f(xs...)
            x[i] = tmp + im*δ/2
            y2 = f(xs...)
            x[i] = tmp
            Δ[i] += im*(y2-y1)/δ
        end
    end
    return grads
end

##==----------------------------------------------------------------------------------------------------------------------

@testset ExtendedTestSet "gate gradients density matrix" begin

    @testset "single qubit gates" begin

        # fictitious density matrix
        ρ = DensityMatrix(randn(Float64, 4), 1)

        # fictitious gradients of cost function with respect to output density matrix
        Δ = DensityMatrix(randn(Float64, 4), 1)

        for g in [RxGate, RyGate, RzGate]
            θval = 2π*rand()

            f = (θ)->dot(Δ.v, apply(ρ, CircuitGate((1,), g(θ[]))).v)
            ngrad = ngradient(f, [θval])
            dg = Qaintessent.backward_density(g(θval), reshape(kron(ρ.v, Δ.v), 4, 4))
            @test isapprox(dg.θ, ngrad[1], rtol=1e-5, atol=1e-5)

            fa = (θ)->dot(Δ.v, Qaintessent.apply_mixed_add(ρ, CircuitGate((1,), g(θ[]))).v)
            ngrad = ngradient(fa, [θval])
            dg = Qaintessent.backward_density_mixed_add(g(θval), reshape(kron(ρ.v, Δ.v), 4, 4))
            @test isapprox(dg.θ, ngrad[1], rtol=1e-5, atol=1e-5)

            fs = (θ)->dot(Δ.v, Qaintessent.apply_mixed_sub(ρ, CircuitGate((1,), g(θ[]))).v)
            ngrad = ngradient(fs, [θval])
            dg = Qaintessent.backward_density_mixed_sub(g(θval), reshape(kron(ρ.v, Δ.v), 4, 4))
            @test isapprox(dg.θ, ngrad[1], rtol=1e-5, atol=1e-5)
        end

        for g in [XGate, YGate, ZGate, HadamardGate, SGate, SdagGate, TGate, TdagGate]
            dg = Qaintessent.backward_density(g(), reshape(kron(ρ.v, Δ.v), 4, 4))
            @test isapprox(dg, g())
        end

        begin
            nθval = randn(3)

            f = (nθ)->dot(Δ.v, apply(ρ, CircuitGate((1,), RotationGate(nθ))).v)
            ngrad = ngradient(f, nθval)
            dg = Qaintessent.backward_density(RotationGate(nθval), reshape(kron(ρ.v, Δ.v), 4, 4))
            @test isapprox(dg.nθ, ngrad[1], rtol=1e-5, atol=1e-5)

            fa = (nθ)->dot(Δ.v, Qaintessent.apply_mixed_add(ρ, CircuitGate((1,), RotationGate(nθ))).v)
            ngrad = ngradient(fa, nθval)
            dg = Qaintessent.backward_density_mixed_add(RotationGate(nθval), reshape(kron(ρ.v, Δ.v), 4, 4))
            @test isapprox(dg.nθ, ngrad[1], rtol=1e-5, atol=1e-5)

            fs = (nθ)->dot(Δ.v, Qaintessent.apply_mixed_sub(ρ, CircuitGate((1,), RotationGate(nθ))).v)
            ngrad = ngradient(fs, nθval)
            dg = Qaintessent.backward_density_mixed_sub(RotationGate(nθval), reshape(kron(ρ.v, Δ.v), 4, 4))
            @test isapprox(dg.nθ, ngrad[1], rtol=1e-5, atol=1e-5)

            # special case: zero vector

            ngrad = ngradient(f, zeros(3))
            dg = Qaintessent.backward_density(RotationGate(zeros(3)), reshape(kron(ρ.v, Δ.v), 4, 4))
            @test isapprox(dg.nθ, ngrad[1], rtol=1e-5, atol=1e-5)

            ngrad = ngradient(fa, zeros(3))
            dg = Qaintessent.backward_density_mixed_add(RotationGate(zeros(3)), reshape(kron(ρ.v, Δ.v), 4, 4))
            @test isapprox(dg.nθ, ngrad[1], rtol=1e-5, atol=1e-5)

            ngrad = ngradient(fs, zeros(3))
            dg = Qaintessent.backward_density_mixed_sub(RotationGate(zeros(3)), reshape(kron(ρ.v, Δ.v), 4, 4))
            @test isapprox(dg.nθ, ngrad[1], rtol=1e-5, atol=1e-5)
        end

        begin
            ϕval = 2π*rand()

            f = (ϕ)->dot(Δ.v, apply(ρ, CircuitGate((1,), PhaseShiftGate(ϕ[]))).v)
            ngrad = ngradient(f, [ϕval])
            dg = Qaintessent.backward_density(PhaseShiftGate(ϕval), reshape(kron(ρ.v, Δ.v), 4, 4))
            @test isapprox(dg.ϕ, ngrad[1], rtol=1e-5, atol=1e-5)

            fa = (ϕ)->dot(Δ.v, Qaintessent.apply_mixed_add(ρ, CircuitGate((1,), PhaseShiftGate(ϕ[]))).v)
            ngrad = ngradient(fa, [ϕval])
            dg = Qaintessent.backward_density_mixed_add(PhaseShiftGate(ϕval), reshape(kron(ρ.v, Δ.v), 4, 4))
            @test isapprox(dg.ϕ, ngrad[1], rtol=1e-5, atol=1e-5)

            fs = (ϕ)->dot(Δ.v, Qaintessent.apply_mixed_sub(ρ, CircuitGate((1,), PhaseShiftGate(ϕ[]))).v)
            ngrad = ngradient(fs, [ϕval])
            dg = Qaintessent.backward_density_mixed_sub(PhaseShiftGate(ϕval), reshape(kron(ρ.v, Δ.v), 4, 4))
            @test isapprox(dg.ϕ, ngrad[1], rtol=1e-5, atol=1e-5)
        end
    end

    @testset "entanglement gates" begin

        # fictitious density matrix
        ρ = DensityMatrix(randn(Float64, 16), 2)

        # fictitious gradients of cost function with respect to output density matrix
        Δ = DensityMatrix(randn(Float64, 16), 2)

        for g in [EntanglementXXGate, EntanglementYYGate, EntanglementZZGate]
            θval = 2π*rand()

            f = (θ)->dot(Δ.v, apply(ρ, CircuitGate((1, 2), g(θ[]))).v)
            ngrad = ngradient(f, [θval])
            dg = Qaintessent.backward_density(g(θval), reshape(kron(ρ.v, Δ.v), 16, 16))
            @test isapprox(dg.θ, ngrad[1], rtol=1e-5, atol=1e-5)

            fa = (θ)->dot(Δ.v, Qaintessent.apply_mixed_add(ρ, CircuitGate((1, 2), g(θ[]))).v)
            ngrad = ngradient(fa, [θval])
            dg = Qaintessent.backward_density_mixed_add(g(θval), reshape(kron(ρ.v, Δ.v), 16, 16))
            @test isapprox(dg.θ, ngrad[1], rtol=1e-5, atol=1e-5)

            fs = (θ)->dot(Δ.v, Qaintessent.apply_mixed_sub(ρ, CircuitGate((1, 2), g(θ[]))).v)
            ngrad = ngradient(fs, [θval])
            dg = Qaintessent.backward_density_mixed_sub(g(θval), reshape(kron(ρ.v, Δ.v), 16, 16))
            @test isapprox(dg.θ, ngrad[1], rtol=1e-5, atol=1e-5)
        end
    end

    @testset "controlled gates 1" begin

        # fictitious density matrix
        ρ = DensityMatrix(randn(Float64, 64), 3)

        # fictitious gradients of cost function with respect to output density matrix
        Δ = DensityMatrix(randn(Float64, 64), 3)

        θval = 2π*rand()

        for g in [RxGate, RyGate, RzGate]

            f1(θ) = dot(Δ.v, apply(ρ, CircuitGate((1, 2, 3), ControlledGate{g}(g(θ[]), 2))).v)
            ngrad = ngradient(f1, [θval])
            dg = Qaintessent.backward_density(ControlledGate{g}(g(θval), 2), reshape(kron(ρ.v, Δ.v), 64, 64))
            #@code_warntype(Qaintessent.backward_density(ControlledGate{RyGate}(RyGate(θval), 2), reshape(kron(ρ.v, Δ.v), 64, 64)))
            @test isapprox(dg.U.θ, ngrad[1], rtol=1e-5, atol=1e-5)
            if !isapprox(dg.U.θ, ngrad[1], rtol=1e-5, atol=1e-5)
                println(θval)
                println(norm(dg.U.θ - ngrad[1]) / norm(dg.U.θ))
            end
        end

        for g in [EntanglementXXGate, EntanglementYYGate, EntanglementZZGate]
            f2(θ) = dot(Δ.v, apply(ρ, CircuitGate((1, 2, 3), ControlledGate{g}(g(θ[]), 1))).v)
            ngrad = ngradient(f2, [θval])
            dg = Qaintessent.backward_density(ControlledGate{g}(g(θval), 1), reshape(kron(ρ.v, Δ.v), 64, 64))
            @test isapprox(dg.U.θ, ngrad[1], rtol=1e-5, atol=1e-5)
            if !isapprox(dg.U.θ, ngrad[1], rtol=1e-5, atol=1e-5)
                println(θval)
                println(norm(dg.U.θ - ngrad[1]) / norm(dg.U.θ))
            end
        end

        for g in [XGate, YGate, ZGate, HadamardGate, SGate, SdagGate, TGate, TdagGate]
            dg = Qaintessent.backward_density(ControlledGate{g}(g(), 2), reshape(kron(ρ.v, Δ.v), 64, 64))
            @test isapprox(dg, ControlledGate{g}(g(), 2))
        end
    end

    @testset "controlled gates 2" begin

        # fictitious density matrix
        ρ = DensityMatrix(randn(Float64, 256), 4)

        # fictitious gradients of cost function with respect to output density matrix
        Δ = DensityMatrix(randn(Float64, 256), 4)

        θval = 2π*rand()

        f1(θ) = dot(Δ.v, apply(ρ, CircuitGate((1, 2, 3, 4), ControlledGate{PhaseShiftGate}(PhaseShiftGate(θ[]), 3))).v)
        ngrad = ngradient(f1, [θval])
        dg = Qaintessent.backward_density(ControlledGate{PhaseShiftGate}(PhaseShiftGate(θval), 3), reshape(kron(ρ.v, Δ.v), 256, 256))
        @test isapprox(dg.U.ϕ, ngrad[1], rtol=1e-5, atol=1e-5)

        f2(θ) = dot(Δ.v, apply(ρ, CircuitGate((1, 2, 3, 4), ControlledGate{EntanglementYYGate}(EntanglementYYGate(θ[]), 2))).v)
        ngrad = ngradient(f2, [θval])
        dg = Qaintessent.backward_density(ControlledGate{EntanglementYYGate}(EntanglementYYGate(θval), 2), reshape(kron(ρ.v, Δ.v), 256, 256))
        @test isapprox(dg.U.θ, ngrad[1], rtol=1e-5, atol=1e-5)
    end
end

##==----------------------------------------------------------------------------------------------------------------------

@testset ExtendedTestSet "circuit gradients density matrix" begin

    @testset "individual circuit gate" begin
        # fictitious density matrix
        ρ = DensityMatrix(randn(Float64, 64), 3)

        # fictitious gradients of cost function with respect to output density matrix
        Δ = DensityMatrix(randn(Float64, 64), 3)

        θval = 2π*rand()

        cg(θ) = CircuitGate((3, 1), EntanglementYYGate(θ))
        f(θ) = dot(Δ.v, apply(ρ, cg(θ[])).v)
        ngrad = ngradient(f, [θval])
        dcg = Qaintessent.backward_density(cg(θval), ρ, Δ)
        @test isapprox(dcg.gate.θ, ngrad[1], rtol=1e-5, atol=1e-5)
    end

    @testset "moments" begin
        # fictitious density matrix
        ρ = DensityMatrix(randn(Float64, 256), 4)

        # fictitious gradients of cost function with respect to output density matrix
        Δ = DensityMatrix(0.1*randn(Float64, 256), 4)

        moments(θ, κ, ϕ) = [Moment([circuit_gate(2, RyGate(θ)), circuit_gate(3, 1, EntanglementZZGate(κ))]),
                            Moment(circuit_gate(3, PhaseShiftGate(ϕ), (1, 4)))]

        θval = 2π*rand()
        κval = 2π*rand()
        ϕval = 2π*rand()

        f(rθ, rκ, rϕ) = dot(Δ.v, apply(ρ, moments(rθ[], rκ[], rϕ[])).v)
        ngrad = ngradient(f, [θval], [κval], [ϕval])
        ρs = apply(ρ, moments(θval, κval, ϕval))
        dm, Δ0 = Qaintessent.backward_density(moments(θval, κval, ϕval), ρs, Δ)

        @test isapprox(dm[1].gates[1].gate.θ, ngrad[1], rtol=1e-5, atol=1e-5)
        @test isapprox(dm[1].gates[2].gate.θ, ngrad[2], rtol=1e-5, atol=1e-5)
        @test isapprox(dm[2].gates[1].gate.U.ϕ, ngrad[3], rtol=1e-5, atol=1e-5)
    end

    @testset "full circuit" begin

        N = 4
        A = randn(ComplexF64, 2 ,2)
        U, R = qr(A)
        U = Array(U)
        i = rand(1:N)

        cgc(θ, ϕ, χ, κ, ωn) = CircuitGate[
            circuit_gate(3, HadamardGate()),
            circuit_gate(2, RzGate(θ), (1, 4)),
            circuit_gate(2, TGate()),
            circuit_gate(2, 3, SwapGate()),
            circuit_gate(2, SdagGate(), 1),
            circuit_gate(3, PhaseShiftGate(ϕ)),
            circuit_gate(1, RyGate(χ)),
            circuit_gate(4, HadamardGate(), 2),
            circuit_gate(3, X, 1),
            circuit_gate(2, Z, 4),
            circuit_gate(4, Y),
            circuit_gate(3, 1, EntanglementYYGate(κ)),
            circuit_gate(3, RotationGate(ωn)),
            circuit_gate(4, SGate()),
            circuit_gate(2, TdagGate()),
            circuit_gate(i, MatrixGate(U)),
        ]
        # measurement operators
        meas(M) = MeasurementOperator.([Matrix{Float64}(I, 2^N, 2^N), Hermitian(M)], (Tuple(1:N),))

        # parameter values
        θval = 1.5π
        ϕval  = 0.3
        χval  = √2
        κval  = exp(-0.4)
        ωnval = 0.1π * randn(Float64, 3)
        Mval = randn(ComplexF64, 2^N, 2^N)
        Mval = 0.5*(Mval + adjoint(Mval))

        c = Circuit{N}(cgc(θval, ϕval, χval, κval, ωnval), meas(Mval))

        # input density matrix
        ρ = DensityMatrix(randn(Float64, 4^N), N)
        # fictitious gradients of cost function with respect to circuit output
        Δ = [0.3, -0.2]

        f(rθ, rϕ, rχ, rκ, ωn, M) = dot(Δ, apply(ρ, Circuit{N}(cgc(rθ[], rϕ[], rχ[], rκ[], ωn), meas(M))))
        # numeric gradients
        ngrad = ngradient(f, [θval], [ϕval], [χval], [κval], ωnval, Mval)
        # symmetrize gradient with respect to measurement operator
        ngrad[end][:] = 0.5*(ngrad[end] + adjoint(ngrad[end]))
        dc = Qaintessent.gradients(c, ρ, Δ)[1]

        @test all(isapprox.(ngrad,
        (dc.moments[1][2].gate.U.θ,
         dc.moments[4][2].gate.ϕ,
         dc.moments[5][1].gate.θ,
         dc.moments[7][2].gate.θ,
         dc.moments[8][1].gate.nθ,
         sparse_matrix(dc.meas[2])), rtol=1e-5, atol=1e-5))
    end
end