text:
using TriangleMesh
using Plots
using Triplot
using LinearAlgebra
using SparseArrays
using ForwardDiff
using UnPack
include("./TriFEMUtils.jl") # plotting and misc routines

poly = polygon_unitSquare()
# poly = polygon_regular(5)
# poly = polygon_Lshape()
max_triangle_area = .05/64
mesh = create_mesh(poly, quality_meshing=true,
                   add_switches="penva"*string(max_triangle_area)*"q") # switches

# manufactured solution
uexact(x,y) = exp(sin(1+pi*x)*sin(pi*y))
uBC(x,y) = uexact(x,y)
dudx(x,y) = ForwardDiff.derivative(x->uexact(x,y),x)
dudy(x,y) = ForwardDiff.derivative(y->uexact(x,y),y)
f(x,y) = -(ForwardDiff.derivative(x->dudx(x,y),x) +
           ForwardDiff.derivative(y->dudy(x,y),y)) # -(du^2/dx^2 + du^2/dy^2)

# define Dirichlet and Neumann boundaries
on_Neumann_boundary(x,y) = x ≈ 1 && y > 0 && y < 1 # right face x = 1 (excluding corners)
on_Dirichlet_boundary(x,y) = !on_Neumann_boundary(x,y) # Dirichlet boundary = everywhere else
∇u_dot_n(x,y) = dudx(x,y) # note that n = [1,0] on the Neumann face

# ordering of reference element vertices
# 3
# |`.
# 1--2
reference_vertices = [[-1,1,-1],[-1,-1,1]] # stored as [r,s]
reference_face_indices = [[1,2],[2,3],[1,3]] # ordering of faces for the reference element
reference_elem_info = (;reference_face_indices,reference_vertices)

# define reference basis functions
λ1(r,s) = -(r+s)/2
λ2(r,s) = (1+r)/2
λ3(r,s) = (1+s)/2
λ(r,s) = [λ1.(r,s) λ2.(r,s) λ3.(r,s)]
dλr() = [-.5 .5 0.0]
dλs() = [-.5 0.0 .5]

# x,y = lists of vertices
function compute_geometric_terms(x,y)
    # dx(r,s)/dr = dλ1/dr*x1 + dλ2/dr*x2 + dλ3/dr*x3 = [dλ1/dr dλ2/dr dλ3/dr] * [x1; x2; x3]
    dxdr,dydr = dλr()*x, dλr()*y
    dxds,dyds = dλs()*x, dλs()*y
    G = [dxdr dxds; dydr dyds] # G*r = x, r = inv(G)*x
    J = det(G) # change of variables det(Jacobian) for integration
    drdx, dsdx, drdy, dsdy = inv(G)
    return J, drdx, dsdx, drdy, dsdy
end

# assemble stiffness matrix and RHS vector
function assemble_FE_matrix(mesh)
    VX,VY,EToV = unpack_mesh_info(mesh)
    num_vertices = length(VX)
    num_elements = size(EToV,2) # number of elements = of columns

    A = spzeros(num_vertices, num_vertices)
    b = zeros(num_vertices)
    for e = 1:num_elements # loop through all elements
        ids = EToV[:,e] # vertex ids = local to global index maps
        xv,yv = VX[ids],VY[ids]

        # compute geometric mappings
        J,drdx,dsdx,drdy,dsdy = compute_geometric_terms(xv,yv)
        dλdx = drdx*dλr() + dsdx*dλs() # dλr() = [dλ1/dr dλ2/dr dλ3/dr]
        dλdy = drdy*dλr() + dsdy*dλs() # dλs() = [dλ1/ds dλ2/ds dλ3/ds]
        reference_elem_area = 2.0
        @. A[ids,ids] += J*reference_elem_area*(dλdx'*dλdx + dλdy'*dλdy) # (dλdx'*dλdx)_ij = dλj/dx * dλi/dx

        # midpoint quadrature rule
        x_mid,y_mid = sum(xv)/3, sum(yv)/3 # midpoint rule = averages of vertex locations
        w_mid = 2.0
        b[ids] .+= J*vec(λ(-1/3,-1/3))*(w_mid.*f(x_mid,y_mid)) # λ at midpt = 1/3 constant
    end
    return A,b
end

function impose_Neumann_BCs!(b,mesh,reference_elem_info)

    @unpack reference_face_indices, reference_vertices = reference_elem_info
    VX,VY,EToV = unpack_mesh_info(mesh)
    boundary_indices, boundary_faces = get_boundary_info(reference_face_indices, mesh)
    # contributions from Neumann BCs
    for (f,e) in boundary_faces
        ref_fids = reference_face_indices[f]
        fids = EToV[ref_fids,e]
        xf,yf = VX[fids],VY[fids]
        x_mid,y_mid = sum(xf)/2,sum(yf)/2

        if on_Neumann_boundary(x_mid,y_mid)
            r,s = reference_vertices
            r_mid, s_mid = sum(r[ref_fids])/2, sum(s[ref_fids])/2
            w_mid = 2.0

            ids = EToV[:,e]
            face_length = sqrt((xf[1]-xf[2])^2 + (yf[1]-yf[2])^2)
            b[ids] .+= face_length/2 * vec(λ(r_mid,s_mid)) * w_mid * ∇u_dot_n(x_mid,y_mid) # ∫ ∇u⋅n ϕ_i on each Neumann face
        end
    end
    return b
end
function impose_Dirichlet_BCs!(A,b,mesh,reference_elem_info)

    @unpack reference_face_indices, reference_vertices = reference_elem_info
    VX,VY,EToV = unpack_mesh_info(mesh)
    boundary_indices, boundary_faces = get_boundary_info(reference_face_indices, mesh)

    # impose Dirichlet BCs
    for i in boundary_indices
        xi,yi = VX[i],VY[i]
        if on_Dirichlet_boundary(xi,yi)
            b -= Vector(A[:,i]*uBC(xi,yi))
            A[:,i] .= 0
            A[i,:] .= 0
            A[i,i] = 1.0
        end
    end
    for i in boundary_indices
        xi,yi = VX[i],VY[i]
        if on_Dirichlet_boundary(xi,yi)
            b[i] = uBC(xi,yi)
        end
    end
    return A,b
end

A,b = assemble_FE_matrix(mesh)
b   = impose_Neumann_BCs!(b,mesh,reference_elem_info)
A,b = impose_Dirichlet_BCs!(A,b,mesh,reference_elem_info)
u = A\b
VX,VY,EToV = unpack_mesh_info(mesh)
@show maximum(abs.(u .- uexact.(VX,VY)))
triplot(VX,VY,u .- uexact.(VX,VY),EToV)
# triplot(VX,VY,uexact.(VX,VY),EToV)


text:
using Roots
function Yfind(Q,P,H)

    Hamil(XX,YY,QQ,PP)=( (QQ-XX)^2+(PP-YY)^2 )*( (QQ+XX)^2+(PP+YY)^2 )/((PP^4+2*PP^2*(XX^2-1)+( 1+XX^2)^2 )*(QQ^4+2*QQ^2*(YY^2+1)+(YY^2-1)^2 ))
    Y_find(y)=Hamil(0,y,Q,P)-H
    try
        Y=find_zero(Y_find,.01,maxeval=100,maxfnevals=300,tol=1e-15)
    catch
        Y=zeros(0)
    end
end
# @everywhere function Yfind(Q,P,H)
#     Y_find1(y)=Hamil(0,y,Q,P)-H
#         try
#             Y=find_zeros(Y_find1,0,6, maxeval=100,maxfnevals=300,tol=1e-12)
#          catch
#             Y=zeros(0)
#         end
#     end

ext: jl


text:
# abstract type Matrix{T} <: AbstractMatrix{T} end

"""
# Summary
    Embedding{T} <: AbstractMatrix{T}

# Fields
    X::Matrix{T} # Embedding  (each column is an item of the Embedding)

# Examples
Building an embedding from the desired dimensions:

    Embedding(d::Int, n::Int)

Construct a random d × n Embedding.

---
Building an embedding from a tuple:

    Embedding(t::Tuple{Int,Int})

Construct a random n × d Embedding from a tuple.
If typeof(tuple) == Tuple{Int64}, the embedding
is reshaped such that it is a matrix of 1 × n.

---
Building an embedding from a matrix:

    Embedding(X₀::Matrix{T}) where T <: Real

Construct an Embedding from a matrix of points.

We assume that the points are in a d × n matrix, where
  - d: dimensions
  - n: Number of points

If eltype(X₀) == Int, an AbstractFloat is forced.
"""
mutable struct Embedding{T} <: AbstractMatrix{T}
    X::Matrix{T} # Embedding

    function Embedding(d::Int, n::Int; σ::Float64=0.0001)
        @assert n ≥ d "n ≥ d is required"

        X = σ * randn(d, n)

        new{eltype(X)}(X)
    end

    function Embedding(t::Tuple{Int,Int}; σ::Float64=0.0001)
        @assert t[2] ≥ t[1] "n ≥ d is required"

        X = σ * randn(Float64, t)

        new{Float64}(X)
    end

    function Embedding(t::Tuple{Int}; σ::Float64=0.0001)
        X = σ * randn(Float64, 1, t[1])

        if maximum(X) > 1e-3
          @warn "Norm of X might be too large for initialization. Values should be O(1e-5)."
        end

        new{Float64}(X)
    end

    function Embedding(X₀::AbstractMatrix{T}) where T <: Real
        d, n = size(X₀)
        @assert n ≥ d "n ≥ d is required"

        # X = X₀ * V(d)
        X = X₀ .+ 0.0 # This forces X to be a Matrix{AbstractFloat}

        new{eltype(X)}(X)
    end

    function Embedding(X₀::AbstractVector{T}) where T <: Real
        X = reshape(X₀ .+ 0.0, 1, length(X₀)) # This forces X to be a Matrix{AbstractFloat}

        new{eltype(X)}(X)
    end

end

Base.size(X::Embedding) = size(X.X)
Base.getindex(X::Embedding, inds...) = getindex(X.X, inds...)
function Base.Vector(X::Embedding)
  size(X, 1) == 1 || throw(ArgumentError("X must be of size (1, n)."))

  return dropdims(Matrix(X.X'), dims=2)
end

"""
    ndims(X::Matrix)

Obtain the number of dimensions in the Embedding X. If X is d × n,
ndims(X) returns d.

# Examples
    julia> X = Embedding(2,10);

    julia> ndims(X)
    2
"""
Base.ndims(X::Embedding) = size(X,1)

"""
    nitems(X::Matrix)

Obtain the number of items in the Embedding X. If X is d × n,
nitems(X) returns n.

# Examples
    julia> X = Embedding(2,10);

    julia> nitems(X)
    10
"""
nitems(X::Embedding) = size(X,2)

"""
    V(n::Int)

Compute a normalizing matrix of dimension n. The matrix V centers the rows or columns
of a matrix that is pre or post-multiplied by V.

# Examples

    julia> V(10)
    10×10 Array{Float64,2}:
      0.9  -0.1  -0.1  -0.1  -0.1  -0.1  -0.1  -0.1  -0.1  -0.1
     -0.1   0.9  -0.1  -0.1  -0.1  -0.1  -0.1  -0.1  -0.1  -0.1
     -0.1  -0.1   0.9  -0.1  -0.1  -0.1  -0.1  -0.1  -0.1  -0.1
     -0.1  -0.1  -0.1   0.9  -0.1  -0.1  -0.1  -0.1  -0.1  -0.1
     -0.1  -0.1  -0.1  -0.1   0.9  -0.1  -0.1  -0.1  -0.1  -0.1
     -0.1  -0.1  -0.1  -0.1  -0.1   0.9  -0.1  -0.1  -0.1  -0.1
     -0.1  -0.1  -0.1  -0.1  -0.1  -0.1   0.9  -0.1  -0.1  -0.1
     -0.1  -0.1  -0.1  -0.1  -0.1  -0.1  -0.1   0.9  -0.1  -0.1
     -0.1  -0.1  -0.1  -0.1  -0.1  -0.1  -0.1  -0.1   0.9  -0.1
     -0.1  -0.1  -0.1  -0.1  -0.1  -0.1  -0.1  -0.1  -0.1   0.9

---

    julia> using Statistics: mean

    julia> A = rand(2,2); mean(A, dims=1)
    1×2 Array{Float64,2}:
     0.502544  0.538744

    julia> mean(V(2)*A, dims=1)
    1×2 Array{Float64,2}:
     0.0  0.0

    julia> mean(A*V(2), dims=2)
    2×1 Array{Float64,2}:
     0.0
     0.0

"""
V(n::Int) = Matrix(I, n, n) - ones(n,n)/n

Gram(X::Embedding) = X'X

function L(n::Int, (i,j,k))
  L_t = spzeros(n,n)
  L_t[i,j] = -1
  L_t[i,k] = 1
  L_t[j,i] = -1
  L_t[j,j] = 1
  L_t[k,i] = 1
  L_t[k,k] = -1
  return L_t
end


text:
ENV["GKSwstype"]="100"

using DifferentiableNAS
using Flux
using Flux: throttle, logitcrossentropy, onecold, onehotbatch
using StatsBase: mean
using Parameters
using CUDA
using Distributions
using BSON
using Plots
gr()
include("CIFAR10.jl")
include("training_utils.jl")
include("visualize.jl")

#folder_name = "test/models/osirim/darts_2021-01-13T14:17:49.647"
#folder_name = "test/models/osirim/darts_2021-01-13T14:24:11.711"
folder_name = "test/models/osirim/admm_2021-01-13T17:43:21.59"

trials = Dict("test/models/osirim/darts_2021-01-13T14:17:49.647" => "DARTS_none",
                "test/models/osirim/darts_2021-01-13T14:24:11.711"=> "DARTS_sans_none",
                "test/models/osirim/admm_2021-01-13T17:43:21.59"=> "ADMM_sans_none",
                "test/models/osirim/admm_2021-01-14T16:55:27.481"=> "ADMM_none",
                "test/models/osirim/admm_2021-01-14T17:16:51.535"=> "ADMM_128_RTX",
                "test/models/osirim/admm_2021-01-15T12:06:38.741"=> "ADMM_128_RTX_sans_none_disc",)

for folder_name in keys(trials)

    file_name = string(folder_name, "/histbatch.bson")
    BSON.@load file_name histbatch

    prim = [
        "none",
        "max_pool_3x3",
        "avg_pool_3x3",
        "skip_connect",
        "sep_conv_3x3",
        "sep_conv_5x5",
        #"sep_conv_7x7",
        "dil_conv_3x3",
        "dil_conv_5x5",
        #"conv_7x1_1x7"
    ]

    nodelabels = ["c_{k-2}", "c_{k-1}", "1", "2", "3", "4", "c_{k}"]

    connects = vcat([[(j,i) for j = 1:i-1] for i = 3:6]...)
    if typeof(histbatch)<:histories
        normal_ = histbatch.normal_αs
        reduce_ = histbatch.reduce_αs

        if length(normal_[1][1]) == 7
            prim = prim[2:8]
        end
        n_y_min = minimum([softmax(a[i])[j] for a in histbatch.normal_αs for i in 1:14 for j in 1:length(prim)])
        n_y_max = maximum([softmax(a[i])[j] for a in histbatch.normal_αs for i in 1:14 for j in 1:length(prim)])
        p = Vector(undef, 14)
        for i = 1:14
            p[i] = plot(title = string("Op ",nodelabels[connects[i][1]],"->",nodelabels[connects[i][2]]), ylim=(n_y_min,n_y_max), legend=false)
            for j = 1:length(prim)
                plot!([softmax(a[i])[j] for a in histbatch.normal_αs], label=prim[j])
            end
        end
        plot(p..., layout = (2,7), size = (2200,600));
        savefig(string(folder_name, "/normal_plots.png"))

        r_y_min = minimum([softmax(a[i])[j] for a in histbatch.reduce_αs for i in 1:14 for j in 1:length(prim)])
        r_y_max = maximum([softmax(a[i])[j] for a in histbatch.reduce_αs for i in 1:14 for j in 1:length(prim)])
        p = Vector(undef, 14)
        for i = 1:14
            p[i] = plot(title = string("Op ",nodelabels[connects[i][1]],"->",nodelabels[connects[i][2]]), ylim=(r_y_min,r_y_max), legend=false)
            for j = 1:length(prim)
                plot!([softmax(a[i])[j] for a in histbatch.reduce_αs], label=prim[j])
            end
        end
        plot(p..., layout = (2,7), size = (2200,600));
        savefig(string(folder_name, "/reduce_plots.png"))
    else
        normal_ = histbatch.normal_αs_sm
        reduce_ = histbatch.reduce_αs_sm
        if length(normal_[1][1]) == 7
            prim = prim[2:8]
        else
            n_y_min = minimum([a[i][j] for a in normal_ for i in 1:14 for j in 2:length(prim)])
            n_y_max = maximum([a[i][j] for a in normal_ for i in 1:14 for j in 2:length(prim)])
            p = Vector(undef, 15)
            for i = 1:14
                p[i] = plot(title = string("Op ",nodelabels[connects[i][1]],"->",nodelabels[connects[i][2]]), ylim=(n_y_min,n_y_max), legend=false)
                for j = length(prim):-1:2
                    plot!([a[i][j] for a in normal_], label=prim[j])
                end
            end
            p[15] = plot(showaxis=false, grid = false, xlims=(-1,-1))
            for j = length(prim):-1:1
                plot!([a[1][j] for a in normal_], label=prim[j])
            end
            plot(p..., layout = (3,5), size = (1200,800));
            savefig(string(folder_name, "/", trials[folder_name], "-normal_plots_sans_none.png"))

            r_y_min = minimum([a[i][j] for a in reduce_ for i in 1:14 for j in 1:length(prim)])
            r_y_max = maximum([a[i][j] for a in reduce_ for i in 1:14 for j in 1:length(prim)])
            p = Vector(undef, 15)
            for i = 1:14
                p[i] = plot(title = string("Op ",nodelabels[connects[i][1]],"->",nodelabels[connects[i][2]]), ylim=(r_y_min,r_y_max), legend=false)
                for j = length(prim):-1:2
                    plot!([a[i][j] for a in reduce_], label=prim[j])
                end
            end
            p[15] = plot(showaxis=false, grid = false, xlims=(-1,-1))
            for j = length(prim):-1:2
                plot!([a[1][j] for a in reduce_], label=prim[j])
            end
            plot(p..., layout = (3,5), size = (1200,800));
            savefig(string(folder_name, "/", trials[folder_name], "-reduce_plots_sans_none.png"))
        end
        n_y_min = minimum([a[i][j] for a in normal_ for i in 1:14 for j in 1:length(prim)])
        n_y_max = maximum([a[i][j] for a in normal_ for i in 1:14 for j in 1:length(prim)])
        p = Vector(undef, 15)
        for i = 1:14
            p[i] = plot(title = string("Op ",nodelabels[connects[i][1]],"->",nodelabels[connects[i][2]]), ylim=(n_y_min,n_y_max), legend=false)
            for j = length(prim):-1:1
                plot!([a[i][j] for a in normal_], label=prim[j])
            end
        end
        p[15] = plot(showaxis=false, grid = false, xlims=(-1,-1))
        for j = length(prim):-1:1
            plot!([a[1][j] for a in normal_], label=prim[j])
        end
        plot(p..., layout = (3,5), size = (1200,800));
        savefig(string(folder_name, "/", trials[folder_name], "-normal_plots.png"))

        r_y_min = minimum([a[i][j] for a in reduce_ for i in 1:14 for j in 1:length(prim)])
        r_y_max = maximum([a[i][j] for a in reduce_ for i in 1:14 for j in 1:length(prim)])
        p = Vector(undef, 15)
        for i = 1:14
            p[i] = plot(title = string("Op ",nodelabels[connects[i][1]],"->",nodelabels[connects[i][2]]), ylim=(r_y_min,r_y_max), legend=false)
            for j = length(prim):-1:1
                plot!([a[i][j] for a in reduce_], label=prim[j])
            end
        end
        p[15] = plot(showaxis=false, grid = false, xlims=(-1,-1))
        for j = length(prim):-1:1
            plot!([a[1][j] for a in reduce_], label=prim[j])
        end
        plot(p..., layout = (3,5), size = (1200,800));
        savefig(string(folder_name, "/", trials[folder_name], "-reduce_plots.png"))
        if typeof(histbatch)<:historiessml
            p = Vector(undef, 2)
            p = plot(title = "Losses")
            plot!(histbatch.train_losses, label = "training loss")
            plot!(histbatch.val_losses, label = "validation loss")
            """
            p = Vector(undef, 2)
            p[1] = plot(histbatch.train_losses, label = "training loss")
            p[2] = plot(histbatch.val_losses, label = "validation loss", seriescolor = :red)
            plot(p..., layout = (2,1), size = (600,800), title = "Losses")
            """
            savefig(string(folder_name, "/", trials[folder_name], "-loss.png"))
        end
    end

    visualize(normal_[length(normal_)], string(folder_name, "/", trials[folder_name], "-normal_graph.png"))
    visualize(reduce_[length(reduce_)], string(folder_name, "/", trials[folder_name], "-reduce_graph.png"))
end

max_stars_repo_name: huangyxi/DifferentiableNAS

max_issues_repo_name: huangyxi/DifferentiableNAS

max_forks_repo_name: huangyxi/DifferentiableNAS


text:
module Tinyurl

# Dependencies
using Requests

# Function for shortening the urls, using tinyurl service.
"""
    tinyurl(url::AbstractString)

Shortens the supplied URL.
"""
function tinyurl(url::AbstractString)
 try
  robject = Requests.get("http://tinyurl.com/api-create.php?url=" * url)
  convert(String, robject.data)
 catch Exception
  throw(error("Error! Could not shorten $url"))
 end
end

# Export functions to be accessed!
export tinyurl
end


text:
include(joinpath(pwd(), "examples/implicit_dynamics/utils.jl"))
# include_implicit_dynamics()
include(joinpath(pwd(), "examples/implicit_dynamics/soc_projection_2.jl"))

struct QuadrupedFB{I, T} <: Model{I, T}
      n::Int
      m::Int
      d::Int

      mass          # mass
      inertia       # inertia matrix
      inertia_inv   # inertia matrix inverse
      gravity       # gravity
end

function f(model::QuadrupedFB, z, u, w)
      # states
      x = view(z,1:3)
      r = view(z,4:6)
      v = view(z,7:9)
      ω = view(z,10:12)

      contact_mode = w[12 .+ (1:4)]

      # force in world frame
      F_world = [contact_mode[i] * view(u, (i - 1) * 3 .+ (1:3)) for i = 1:4]

      # torque in body frame
      τ = transpose(MRP(r...)) * sum([cross(w[(i - 1) * 3 .+ (1:3)] - x, F_world[i]) for i = 1:4])

      SVector{12}([v;
                   0.25 * ((1.0 - r' * r) * ω - 2.0 * cross(ω, r) + 2.0*(ω' * r) * r);
                   model.gravity + (1.0 / model.mass) * sum(F_world);
                   model.inertia_inv * (τ - cross(ω, model.inertia * ω))])
end

n, m, d = 12, 12, 16

mass = 1.0
inertia = Diagonal(@SVector[1.0, 1.0, 1.0])
inertia_inv = inv(inertia)

f_max = 10.0
f_min = 0.0
μ_friction = 1.0

# a1
# mass = 108.0 / 9.81
# inertia = 10.0 * Diagonal(@SVector[0.017, 0.057, 0.064])
# inertia_inv = inv(inertia)
# body_height = 0.24

# laikago
# mass = 200.0 / 9.81
# inertia = Diagonal(@SVector[0.07335, 0.25068, 0.25447])
# inertia_inv = inv(inertia)
# body_height = 0.45

SOC_PROJ = true

model = QuadrupedFB{Midpoint, FixedTime}(n, m, d,
                  mass,
                  inertia,
                  inertia_inv,
                  @SVector[0.0, 0.0, -9.81])

function fd(model::QuadrupedFB{Midpoint, FixedTime}, x, u, w, h, t)
	u_proj = SOC_PROJ ? vcat([soc_projection(u[(i - 1) * 3 .+ (1:3)], f_min, f_max, μ_friction) for i = 1:4]...) : u
	rz(z) = fd(model, z, x, u_proj, w, h, t) # implicit midpoint integration
	# x⁺ = newton(rz, copy(x))
    x⁺ = levenberg_marquardt(rz, copy(x))
	return x⁺
end

fd(model, ones(model.n), zeros(model.m), zeros(model.d), h, 1)
function fdx(model::QuadrupedFB{Midpoint, FixedTime}, x, u, w, h, t)
    u_proj = SOC_PROJ ? vcat([soc_projection(u[(i - 1) * 3 .+ (1:3)], f_min, f_max, μ_friction) for i = 1:4]...) : u
    rz(z) = fd(model, z, x, u_proj, w, h, t) # implicit midpoint integration
    # x⁺ = newton(rz, copy(x))
    x⁺ = levenberg_marquardt(rz, copy(x))
    ∇rz = ForwardDiff.jacobian(rz, x⁺)
    rx(z) = fd(model, x⁺, z, u_proj, w, h, t) # implicit midpoint integration
    ∇rx = ForwardDiff.jacobian(rx, x)
    return -1.0 * ∇rz \ ∇rx
end

fdx(model, zeros(model.n), zeros(model.m), zeros(model.d), h, 1)

function fdu(model::QuadrupedFB{Midpoint, FixedTime}, x, u, w, h, t)
    u_proj = SOC_PROJ ? vcat([soc_projection(u[(i - 1) * 3 .+ (1:3)], f_min, f_max, μ_friction) for i = 1:4]...) : u
	u_proj_jac = SOC_PROJ ? cat([soc_projection_jacobian(u[(i - 1) * 3 .+ (1:3)], f_min, f_max, μ_friction) for i = 1:4]..., dims=(1,2)) : I

    rz(z) = fd(model, z, x, u_proj, w, h, t) # implicit midpoint integration
    # x⁺ = newton(rz, copy(x))
    x⁺ = levenberg_marquardt(rz, copy(x))
    ∇rz = ForwardDiff.jacobian(rz, x⁺)
    ru(z) = fd(model, x⁺, x, z, w, h, t) # implicit midpoint integration
    ∇ru = ForwardDiff.jacobian(ru, u_proj)
    return -1.0 * ∇rz \ (∇ru * u_proj_jac)
    # return -1.0 * ∇rz \ ∇ru

	# f(z) = fd(model, x, z, w, h, t)
	# return ForwardDiff.jacobian(f, u)
end

fdu(model, ones(model.n), ones(model.m), ones(model.d), h, 1)

##
# x0 = [q_body; zeros(3); zeros(3)]
# u_gravity = [0.0; 0.0; model.mass * 9.81 / 4.0]
# u_stand = vcat([u_gravity for i = 1:4]...)
# w_stand = [pf1_ref[1]; pf2_ref[2]; pf3_ref[3]; pf4_ref[4]; 1; 1; 1; 1]
#
# # soc_projection(u_gravity)
# T_sim = 101
# x_hist = [x0]
# for t = 1:T_sim-1
#     push!(x_hist, fd(model, x_hist[end], u_stand, w_stand, h, 1))
# end
#
# include(joinpath(pwd(), "models/visualize.jl"))
# vis = Visualizer()
# render(vis)
# visualize!(vis, model, x_hist,
#     [pf1_ref[1] for t = 1:T_sim],
#     [pf2_ref[1] for t = 1:T_sim],
#     [pf3_ref[1] for t = 1:T_sim],
#     [pf4_ref[1] for t = 1:T_sim],
#     Δt = h)
#
# function foot_vis!(vis, i, t;
#     h = 0.01,
#     r = 0.05)
#     if convert(Bool, contact_modes[t][i])
#         setobject!(vis["f$(i)_$(j)_$(t)"],
#             Cylinder(Point(0.0, 0.0, 0.0), Point(0.0, 0.0, h), r),
#             MeshPhongMaterial(color = RGBA(0.0, 1.0, 0.0, 0.01)))
#         settransform!(vis["f$(i)_$(j)_$(t)"], Translation(eval(Symbol("pf$(i)_ref"))[t]))
#     end
#     # setvisible!(vis["f$(i)_$(j)_$(t)"], convert(Bool, contact_modes[t][i]))
# end
#
# function cone_vis!(vis, i, t;
#     h = 0.05,
#     l = h * 1.0,
#     w = h * 1.45)
#
#     pyramid = Pyramid(Point3(0.0, 0.0, 0.0), l, w)
#
#     n = 50
#     if convert(Bool, contact_modes[t][i])
#         for p = 1:n
#             setobject!(vis["pyramid$(i)_$(p)_$(t)"], pyramid,
#                 MeshPhongMaterial(
#                 # color = RGBA(1,153/255,51/255, 1.0))
#                 color = RGBA(51/255,1,1, 0.01))
#                 )
#             settransform!(vis["pyramid$(i)_$(p)_$(t)"],
#                 compose(Translation(eval(Symbol("pf$(i)_ref"))[t][1], eval(Symbol("pf$(i)_ref"))[t][2], l),
#                     LinearMap(RotX(π) * RotZ(π * p / n))),
#                     )
#             setvisible!(vis["pyramid$(i)_$(p)_$(t)"], convert(Bool, contact_modes[t][i]))
#         end
#     end
# end
#
# for i = 1:4
#     for t = 1:T
#         foot_vis!(vis, i, t)
#         cone_vis!(vis, i, t)
#     end
# end

max_issues_repo_name: thowell/DirectMotionPlanning


text:
# module for estimating pose by monte carlo localization
# based on particle filter
# particle num is decided by kld sampling
# parameter nn: range std on straight movement
# parameter no: range std on rotation movement
# parameter on: direction std on straight movement
# parameter oo: direction std on rotation movement

using Distributions, LinearAlgebra, StatsBase

include(joinpath(split(@__FILE__, "src")[1], "src/model/map/map.jl"))
include(joinpath(split(@__FILE__, "src")[1], "src/localization/particle_filter/particle/particle.jl"))
include(joinpath(split(@__FILE__, "src")[1], "src/common/state_transition/state_transition.jl"))
include(joinpath(split(@__FILE__, "src")[1], "src/common/observation_function/observation_function.jl"))

mutable struct KldMcl
  particles
  max_num
  widths
  epsilon
  delta
  bin_num
  motion_noise_rate_pdf
  map
  dist_dev
  dir_dev
  max_likelihood_particle
  estimated_pose

  function KldMcl(init_pose::Array, max_num::Int64;
                  motion_noise_stds::Dict=Dict("nn"=>0.20, "no"=>0.001, "on"=>0.11, "oo"=>0.20),
                  env_map=nothing, dist_dev_rate=0.14, dir_dev=0.05,
                  widths=[0.2, 0.2, pi/18], epsilon=0.1, delta=0.01)
    self = new()
    self.particles = [Particle(init_pose, 1.0)] # initial particle num is one
    self.map = env_map
    self.dist_dev = dist_dev_rate
    self.dir_dev = dir_dev
    self.max_likelihood_particle = self.particles[1]
    self.estimated_pose = self.max_likelihood_particle.pose

    # parameters for kld sampling
    self.widths = widths # width of x,y,theta's bin
    self.max_num = max_num # maximum particle num
    self.epsilon = epsilon # ϵ
    self.delta = delta # δ
    self.bin_num = 0 # num of bin

    v = motion_noise_stds
    c = diagm(0 => [v["nn"]^2, v["no"]^2, v["on"]^2, v["oo"]^2])
    self.motion_noise_rate_pdf = MvNormal(c)
    return self
  end
end

function motion_update(self::KldMcl, speed, yaw_rate, time_interval)
  ws = [(if p.weight <= 0.0 1e-100 else p.weight end) for p in self.particles]

  new_particles = []
  bins = Set()
  for i in 1:self.max_num
    chosen_p = sample(self.particles, Weights(ws)) # choose one particle randomly
    p = deepcopy(chosen_p)
    motion_update(p, speed, yaw_rate, time_interval, self.motion_noise_rate_pdf)
    push!(bins, [Int(floor(e)) for e in p.pose./self.widths])
    push!(new_particles, p)

    if length(bins) > 1
      self.bin_num = length(bins)
    else
      self.bin_num = 2
    end

    if length(new_particles) > ceil(quantile(Chisq(self.bin_num-1), 1-self.delta)/(2*self.epsilon))
      break
    end
  end

  self.particles = new_particles
  for i in 1:length(self.particles) # normalization
    self.particles[i].weight = 1.0/length(self.particles)
  end
end

function set_max_likelihood_pose(self::KldMcl)
  max_index = argmax([p.weight for p in self.particles])
  self.max_likelihood_particle = self.particles[max_index]
  self.estimated_pose = self.max_likelihood_particle.pose
end

function observation_update(self::KldMcl, observation)
  for p in self.particles
    observation_update(p, observation, self.map, self.dist_dev, self.dir_dev)
  end
  set_max_likelihood_pose(self)
end

function draw!(self::KldMcl)
  k = 0.5 # scale for length of arrows
  # all of particles
  px = [p.pose[1] for p in self.particles]
  py = [p.pose[2] for p in self.particles]
  vx = [cos(p.pose[3])*k for p in self.particles]
  vy = [sin(p.pose[3])*k for p in self.particles]
  quiver!(px, py, quiver=(vx, vy), aspect_ratio=true, color="blue")
  # maximum likelihood particle
  mx = [self.estimated_pose[1]]
  my = [self.estimated_pose[2]]
  mvx = [cos(self.estimated_pose[3])*k]
  mvy = [sin(self.estimated_pose[3])*k]
  quiver!(mx, my, quiver=(mvx, mvy), aspect_ratio=true, color="red")
  annotate!(-4.5, -4.5, text("particle:$(length(self.particles)), bin:$(self.bin_num)", :left, :black))
end

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy


text:
using PtFEM, Test

data = Dict(
  # Plane(ndim, nst, nxe, nye, nip, direction, finite_element(nod, nodof), axisymmetric)
  :struc_el => Plane(2, 3, 2, 2, 1, :x, Triangle(3, 2), false),
  :properties => [1.0e6 0.3;],
  :x_coords => [0.0,  0.5,  1.0],
  :y_coords => [0.0,  -0.5,  -1.0],
  :support => [
    (1, [0 1]),
    (4, [0 1]),
    (7, [0 0]),
    (8, [1 0]),
    (9, [1 0])
    ],
  :loaded_nodes => [
    (1, [0.0 -1//4]),
    (2, [0.0 -1//2]),
    (3, [0.0 -1//4])
    ]
)

@time m, dis_df, fm_df = p51(data)
println()

dis_df |> display
println()

println("\nThe integration point (nip = $(data[:struc_el].nip)) stresses are:")
fm_df |> display
println()


text:
using HorizonSideRobots

function main(r::Robot)
    num_vert = moves!(r, Sud)
    num_hor = moves!(r, West)
    print(num_vert)
    print(num_hor)

    height_box = moves!(r, Nord)
    print(height_box)

    z_down(height_box)
    vygol()
    home(num_vert,num_hor)
end

function home(num_vert,num_hor)
    moves!(r, Nord, num_vert)
    moves!(r, Ost, num_hor)
end

function z_down(height_box)
    for i in 1:height_box
        putmarker!(r)
        putmarkers!(r, Ost)
        move!(r,Sud)

        moves!(r,West)
    end
    putmarker!(r)
    putmarkers!(r, Ost)
end

function vygol()
    moves!(r, Sud)
    moves!(r, West)
end

function putmarkers!(r::Robot, side::HorizonSide)
    while isborder(r,side)==false
        move!(r,side)
        putmarker!(r)
    end
end

r=Robot(5,5, animate=true)

max_stars_repo_name: Savko-fokin/mirea-progs

max_issues_repo_name: Savko-fokin/mirea-progs

max_forks_repo_name: Savko-fokin/mirea-progs


text:
__precompile__()

module SpatialFields

import ForwardDiff
import MultiPoly
import Base: convert
import DataStructures: OrderedDict
import StaticArrays: SVector, MVector

export InterpolatingSurface,
    HermiteRadialField,
    PolynomialScalarField,
    PolynomialVectorField

derivative(f::Function) = x -> ForwardDiff.derivative(f, x)
gradient(f::Function) = x -> ForwardDiff.gradient(f, x)

abstract AbstractScalarField{N} <: Function
abstract AbstractVectorField{N} <: Function

type ScalarField{N, F <: Function} <: AbstractScalarField{N}
    func::F
end
ScalarField{F <: Function}(N::Integer, f::F) = ScalarField{N, F}(f)
(field::ScalarField{N}){N}(x::AbstractVector) = field.func(x)

type VectorField{N, F <: Function} <: AbstractVectorField{N}
    func::F
end
VectorField{F <: Function}(N::Integer, f::F) = VectorField{N, F}(f)
(field::VectorField{N}){N}(x::AbstractVector) = field.func(x)

gradient{N}(field::AbstractScalarField{N}) = VectorField(N, x -> ForwardDiff.gradient(field, x))

include("radial_functions.jl")
include("hrbf.jl")
include("interpolating.jl")
include("polynomial.jl")

end # module


text:
using ReversibleSeismic, CUDA
using Test, Pkg

@testset "simulation" begin
    include("simulation.jl")
end

@testset "reversible" begin
    include("reversible.jl")
end

@testset "bennett" begin
    include("bennett.jl")
end

@testset "treeverse" begin
    include("treeverse.jl")
end

function isinstalled(target)
    deps = Pkg.dependencies()
    for (uuid, dep) in deps
        dep.is_direct_dep || continue
        dep.name == target && return true
    end
    return false
end

if CUDA.functional()
    @testset "cuda" begin
        include("cuda.jl")
    end
    @testset "reversible_gpu" begin
        include("reversible_gpu.jl")
    end
end


text:
"""
    Linear Integer Model Scoring System

A LIM scoring system for classification which can present itself
and its predictions nicely.

TODO: Merge this with the existing SLIM scoring systems, for now
I just copied in the LASSO-specific code
"""
struct LIMScoringSystem{S} <: InterpretableClassificationModel
    maxscore::Int               # Max score per rule
    scores::AbstractVector{Int} # Actual scores per rule
    treshold::Int
    colnames::Vector{S}
    classfeaturename::S
end

LIMScoringSystem(maxscore::Int, scores::AbstractVector{Int}, treshold::Int) =
    LIMScoringSystem{String}(maxscore, scores, treshold, String["F$i" for i in 1:length(scores)], "<<unknown>>")

zerooneconvert(v::Number) = (v < 0) ? 0 : 1

function classify(ss::LIMScoringSystem, X::AbstractMatrix{N}) where {N<:Number}
    zerooneconvert.(ss.treshold .+ X * ss.scores)
end

StatsBase.predict(m::LIMScoringSystem, X::AbstractMatrix) = classify(m, X)

# Calculate integer score in range from -L to L, given beta value and absolute max beta.
rescale(beta, absmaxbeta, L) = round(Int, L*beta/absmaxbeta)

function rescale_lasso_coefficients(betas::AbstractVector{Float64}, beta0 = 0; L::Int = 10)
    amb = maximum(abs.(betas))
    rescaledbetas = rescale.(betas, amb, L)
    return rescaledbetas, rescale(beta0, amb, L)
end

# A Lasso-based Linear Integer Model optimizer:
struct LassoLIMOptimizer <: InterpretableModelOptimizer
    withintercept::Bool
    maxnummodels::Int
    maxscore::Int
    LassoLIMOptimizer(; maxscore::Int = 10, withintercept = true, maxnummodels = 30) =
        new(withintercept, maxnummodels, maxscore)
end

function optimize(::Type{<:LIMScoringSystem}, o::LassoLIMOptimizer,
    X::AbstractMatrix, y::AbstractVector)
    lp = fit(LassoPath, X, y, Binomial(), LogitLink(); intercept = o.withintercept)
    nmodels = min(o.maxnummodels, size(lp.coefs, 2))
    wi = o.withintercept
    coefs = unique(map(i ->
        rescale_lasso_coefficients(lp.coefs[:, i], (wi ? lp.b0[i] : 0); L = o.maxscore),
        2:nmodels))
    LIMScoringSystem[LIMScoringSystem(o.maxscore, c...) for c in coefs]
end


text:
@testset "minimum_error" begin

    samples = map(x->gaussian_mixture(x, P₁ = 0.5, μ₁ = 64, σ₁ = 15, P₂ = 0.5, μ₂ = 192, σ₂ = 15), 0:255)
    edges = 0:255
    t = find_threshold(MinimumError(), round.(Int,samples * 10^10), edges)
    @test t == 128
    t = find_threshold(MinimumError(), samples, edges)
    @test t == 128

    # Replicates the result of Figure 2 in the original paper of Kittler and Illingworth.
    samples = map(x->gaussian_mixture(x, P₁ = 0.5, μ₁ = 50, σ₁ = 4, P₂ = 0.5, μ₂ = 150, σ₂ = 30), 0:255)
    edges = 0:255
    t = find_threshold(MinimumError(), ceil.(Int,samples * 10^10), edges)
    @test t == 64
    t = find_threshold(MinimumError(), samples, edges)
    @test t == 64

    img = testimage("cameraman")
    # Uses the underlying UInt8 type when counting the intensity frequencies.
    edges, counts = build_histogram(reinterpret.(gray.(img)), 0:1:255)
    t = find_threshold(MinimumError(),counts[1:end], 0:1:255)
    @test t == 25
end


text:
macro scachefast(type, common_cache_dir, ex)
	Base.depwarn("`scachefast` is deprecated, use `scache` instead.", :scachefast)

	:(@scache $(esc(type)) $(esc(common_cache_dir)) $(esc(ex)))
end


text:
## Extrapolations are currently just constant

function sde_extrapolant!(out,Θ,integrator::DEIntegrator,idxs,deriv::Type)
  if idxs === nothing
    recursivecopy!(out,integrator.u)
  else
    out[idxs] .= @view integrator.u[idxs]
  end
end

function sde_extrapolant(Θ,integrator::DEIntegrator,idxs,deriv::Type)
  if idxs === nothing
    return integrator.u
  else
    return integrator.u[idxs]
  end
end

function sde_interpolant!(out,Θ,integrator::DEIntegrator,idxs,deriv::Type)
  sde_interpolant!(out,Θ,integrator.dt,integrator.uprev,integrator.u,idxs,deriv)
end

function sde_interpolant(Θ,integrator::DEIntegrator,idxs,deriv::Type{T}) where T
  sde_interpolant(Θ,integrator.dt,integrator.uprev,integrator.u,idxs,deriv)
end

@muladd function sde_interpolant(Θ,dt,u0,u1,idxs::Nothing,deriv::Type{Val{0}})
  @. (1-Θ)*u0 + Θ*u1
end

@muladd function sde_interpolant(Θ,dt,u0,u1,idxs,deriv::Type{Val{0}})
  @. (1-Θ)*u0[idxs] + Θ*u1[idxs]
end

function sde_interpolant(Θ,dt,u0,u1,idxs::Nothing,deriv::Type{Val{1}})
  @. (u1-u0)/dt
end

function sde_interpolant(Θ,dt,u0,u1,idxs,deriv::Type{Val{1}})
  @. (u1[idxs]-u0[idxs])/dt
end

@muladd function sde_interpolant!(out,Θ,dt,u0,u1,idxs,deriv::Type{Val{0}})
  Θm1 = (1-Θ)
  if idxs === nothing
    @. out = Θm1*u0 + Θ*u1
  else
    @views @. out = Θm1*u0[idxs] + Θ*u1[idxs]
  end
end

function sde_interpolant!(out,Θ,dt,u0,u1,idxs,deriv::Type{Val{1}})
  if idxs === nothing
    @. out = (u1-u0)/dt
  else
    @views @. out = (u1[idxs]-u0[idxs])/dt
  end
end

@inline function current_interpolant(t::Number,integrator::DEIntegrator,idxs,deriv)
  Θ = (t-integrator.tprev)/integrator.dt
  sde_interpolant(Θ,integrator,idxs,deriv)
end

@inline function current_interpolant(t,integrator::DEIntegrator,idxs,deriv)
  Θ = (t.-integrator.tprev)./integrator.dt
  [sde_interpolant(ϕ,integrator,idxs,deriv) for ϕ in Θ]
end

@inline function current_interpolant!(val,t::Number,integrator::DEIntegrator,idxs,deriv)
  Θ = (t-integrator.tprev)/integrator.dt
  sde_interpolant!(val,Θ,integrator,idxs,deriv)
end

@inline function current_interpolant!(val,t,integrator::DEIntegrator,idxs,deriv)
  Θ = (t.-integrator.tprev)./integrator.dt
  [sde_interpolant!(val,ϕ,integrator,idxs,deriv) for ϕ in Θ]
end

@inline function current_extrapolant(t::Number,integrator::DEIntegrator,idxs=nothing,deriv=Val{0})
  Θ = (t-integrator.tprev)/(integrator.t-integrator.tprev)
  sde_extrapolant(Θ,integrator,idxs,deriv)
end

@inline function current_extrapolant!(val,t::Number,integrator::DEIntegrator,idxs=nothing,deriv=Val{0})
  Θ = (t-integrator.tprev)/(integrator.t-integrator.tprev)
  sde_extrapolant!(val,Θ,integrator,idxs,deriv)
end

@inline function current_extrapolant(t::AbstractArray,integrator::DEIntegrator,idxs=nothing,deriv=Val{0})
  Θ = (t.-integrator.tprev)./(integrator.t-integrator.tprev)
  [sde_extrapolant(ϕ,integrator,idxs,deriv) for ϕ in Θ]
end

@inline function current_extrapolant!(val,t,integrator::DEIntegrator,idxs=nothingderiv=Val{0})
  Θ = (t.-integrator.tprev)./(integrator.t-integrator.tprev)
  [sde_extrapolant!(val,ϕ,integrator,idxs,deriv) for ϕ in Θ]
end

"""
sde_interpolation(tvals,ts,timeseries,ks)

Get the value at tvals where the solution is known at the
times ts (sorted), with values timeseries and derivatives ks
"""
@inline function sde_interpolation(tvals,id,idxs,deriv,p,continuity::Symbol=:left)
  @unpack ts,timeseries = id
  tdir = sign(ts[end]-ts[1])
  idx = sortperm(tvals,rev=tdir<0)
  tdir*tvals[idx[end]] > tdir*ts[end] && error("Solution interpolation cannot extrapolate past the final timepoint. Either solve on a longer timespan or use the local extrapolation from the integrator interface.")
  tdir*tvals[idx[1]] < tdir*ts[1] && error("Solution interpolation cannot extrapolate before the first timepoint. Either start solving earlier or use the local extrapolation from the integrator interface.")
  i = 2 # Start the search thinking it's between ts[1] and ts[2]
  if typeof(idxs) <: Number
    vals = Vector{eltype(first(timeseries))}(undef,length(tvals))
  elseif typeof(idxs) <: AbstractArray
     vals = Vector{Array{eltype(first(timeseries)),ndims(idxs)}}(undef,length(tvals))
  else
    vals = Vector{eltype(timeseries)}(undef,length(tvals))
  end
  @inbounds for j in idx
    t = tvals[j]
    i = searchsortedfirst(@view(ts[i:end]),t,rev=tdir<0)+i-1 # It's in the interval ts[i-1] to ts[i]
    if ts[i] == t
      lasti = lastindex(ts)
      k = continuity == :right && i+1 <= lasti && ts[i+1] == t ? i+1 : i
      if idxs === nothing
        vals[j] = timeseries[k]
      else
        vals[j] = timeseries[k][idxs]
      end
    elseif ts[i-1] == t # Can happen if it's the first value!
      if idxs === nothing
        vals[j] = timeseries[i-1]
      else
        vals[j] = timeseries[i-1][idxs]
      end
    else
      dt = ts[i] - ts[i-1]
      Θ = (t-ts[i-1])/dt
      vals[j] = sde_interpolant(Θ,dt,timeseries[i-1],timeseries[i],idxs,deriv)
    end
  end
  DiffEqArray(vals,tvals)
end

"""
sde_interpolation(tval::Number,ts,timeseries,ks)

Get the value at tval where the solution is known at the
times ts (sorted), with values timeseries and derivatives ks
"""
@inline function sde_interpolation(tval::Number,id,idxs,deriv,p,continuity::Symbol=:left)
  @unpack ts,timeseries = id
  tdir = sign(ts[end]-ts[1])
  tdir*tval > tdir*ts[end] && error("Solution interpolation cannot extrapolate past the final timepoint. Either solve on a longer timespan or use the local extrapolation from the integrator interface.")
  tdir*tval < tdir*ts[1] && error("Solution interpolation cannot extrapolate before the first timepoint. Either start solving earlier or use the local extrapolation from the integrator interface.")
  @inbounds i = searchsortedfirst(ts,tval,rev=tdir<0) # It's in the interval ts[i-1] to ts[i]
  @inbounds if ts[i] == tval
    lasti = lastindex(ts)
    k = continuity == :right && i+1 <= lasti && ts[i+1] == tval ? i+1 : i
    if idxs === nothing
      val = timeseries[k]
    else
      val = timeseries[k][idxs]
    end
  elseif ts[i-1] == tval # Can happen if it's the first value!
    if idxs === nothing
      val = timeseries[i-1]
    else
      val = timeseries[i-1][idxs]
    end
  else
    dt = ts[i] - ts[i-1]
    Θ = (tval-ts[i-1])/dt
    val = sde_interpolant(Θ,dt,timeseries[i-1],timeseries[i],idxs,deriv)
  end
  val
end

@inline function sde_interpolation!(out,tval::Number,id,idxs,deriv,p,continuity::Symbol=:left)
  @unpack ts,timeseries = id
  tdir = sign(ts[end]-ts[1])
  tdir*tval > tdir*ts[end] && error("Solution interpolation cannot extrapolate past the final timepoint. Either solve on a longer timespan or use the local extrapolation from the integrator interface.")
  tdir*tval < tdir*ts[1] && error("Solution interpolation cannot extrapolate before the first timepoint. Either start solving earlier or use the local extrapolation from the integrator interface.")
  @inbounds i = searchsortedfirst(ts,tval,rev=tdir<0) # It's in the interval ts[i-1] to ts[i]
  @inbounds if ts[i] == tval
    lasti = lastindex(ts)
    k = continuity == :right && i+1 <= lasti && ts[i+1] == tval ? i+1 : i
    if idxs === nothing
      copyto!(out,timeseries[k])
    else
      copyto!(out,timeseries[k][idxs])
    end
  elseif ts[i-1] == tval # Can happen if it's the first value!
    if idxs === nothing
      copyto!(out,timeseries[i-1])
    else
      copyto!(out,timeseries[i-1][idxs])
    end
  else
    dt = ts[i] - ts[i-1]
    Θ = (tval-ts[i-1])/dt
    sde_interpolant!(out,Θ,dt,timeseries[i-1],timeseries[i],idxs,deriv)
  end
end

@inline function sde_interpolation!(vals,tvals,id,idxs,deriv,p,continuity::Symbol=:left)
  @unpack ts,timeseries = id
  tdir = sign(ts[end]-ts[1])
  idx = sortperm(tvals,rev=tdir<0)
  tdir*tvals[idx[end]] > tdir*ts[end] && error("Solution interpolation cannot extrapolate past the final timepoint. Either solve on a longer timespan or use the local extrapolation from the integrator interface.")
  tdir*tvals[idx[1]] < tdir*ts[1] && error("Solution interpolation cannot extrapolate before the first timepoint. Either start solving earlier or use the local extrapolation from the integrator interface.")
  i = 2 # Start the search thinking it's between ts[1] and ts[2]
  @inbounds for j in idx
    t = tvals[j]
    i = searchsortedfirst(@view(ts[i:end]),t,rev=tdir<0)+i-1 # It's in the interval ts[i-1] to ts[i]
    if ts[i] == t
      lasti = lastindex(ts)
      k = continuity == :right && i+1 <= lasti && ts[i+1] == t ? i+1 : i
      if idxs === nothing
        vals[j] = timeseries[k]
      else
        vals[j] = timeseries[k][idxs]
      end
    elseif ts[i-1] == t # Can happen if it's the first value!
      if idxs === nothing
        vals[j] = timeseries[i-1]
      else
        vals[j] = timeseries[i-1][idxs]
      end
    else
      dt = ts[i] - ts[i-1]
      Θ = (t-ts[i-1])/dt
      if eltype(timeseries) <: AbstractArray
        sde_interpolant!(vals[j],Θ,dt,timeseries[i-1],timeseries[i],idxs,deriv)
      else
        vals[j] = sde_interpolant(Θ,dt,timeseries[i-1],timeseries[i],idxs,deriv)
      end
    end
  end
end


text:
###
### (::Tensor) ⊗ (::Tensor)
###

using TensorCore

function TensorCore.tensor(t::Tensor{<:Any,N,ipos}, s::Tensor{<:Any,M,jpos}) where {N,ipos,M,jpos}
    ts = tensor(parent(t), parent(s))
    ijpos = IndexPos{N+M}((ipos.x | jpos.x << N))
    Tensor{eltype(ts), N+M, ijpos, typeof(ts)}(ts)
end


text:
import Test: @test, @testset, @test_throws

@testset "Vector plasma" begin
    # scalar modal model, vector modal at 0 degrees and vector modal at 45 degrees should be identical
    using Luna
    import LinearAlgebra: norm
    a = 9e-6
    gas = :Ar
    pres = 5
    flength = 5e-2
    τfwhm = 30e-15
    λ0 = 1500e-9
    energy = 1.7e-6
    grid = Grid.RealGrid(flength, λ0, (600e-9, 3000e-9), 1e-12)
    energyfun, energyfunω = Fields.energyfuncs(grid)
    dens0 = PhysData.density(gas, pres)
    densityfun(z) = dens0
    ionpot = PhysData.ionisation_potential(gas)
    ionrate = Ionisation.ionrate_fun!_ADK(ionpot)
    # scalar modal
    modes = (
        Capillary.MarcatiliMode(a, gas, pres, n=1, m=1, kind=:HE, ϕ=0.0, loss=false),
    )
    nmodes = length(modes)
    plasma = Nonlinear.PlasmaCumtrapz(grid.to, grid.to,
                                    ionrate, ionpot)
    responses = (Nonlinear.Kerr_field(PhysData.γ3_gas(gas)),
                plasma)
    inputs = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy)
    Eω, transform, FT = Luna.setup(grid, densityfun, responses, inputs, modes,
                                :y; full=false)
    linop = LinearOps.make_const_linop(grid, modes, λ0)
    statsfun = Stats.default(grid, Eω, modes, linop, transform; gas=gas)
    outscalar = Output.MemoryOutput(0, grid.zmax, 2, statsfun)
    Luna.run(Eω, grid, linop, transform, FT, outscalar)
    # vector linear 0 degrees
    plasma = Nonlinear.PlasmaCumtrapz(grid.to, Array{Float64}(undef, length(grid.to), 2),
                                      ionrate, ionpot)
    responses = (Nonlinear.Kerr_field(PhysData.γ3_gas(gas)),
             plasma)
    Eω, transform, FT = Luna.setup(grid, densityfun, responses, inputs, modes,
                                   :xy; full=false)
    outvector = Output.MemoryOutput(0, grid.zmax, 2, statsfun)
    Luna.run(Eω, grid, linop, transform, FT, outvector)

    Iωs = abs2.(outscalar.data["Eω"][:,1,1])
    Iωv = abs2.(outvector.data["Eω"][:,1,1])
    @test Iωs ≈ Iωv

    modes = (
        Capillary.MarcatiliMode(a, gas, pres, n=1, m=1, kind=:HE, ϕ=0.0, loss=false),
        Capillary.MarcatiliMode(a, gas, pres, n=1, m=1, kind=:HE, ϕ=π/2, loss=false),
    )
    nmodes = length(modes)
    inputs = ((mode=1, fields=(Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy/2),)),
              (mode=2, fields=(Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy/2),)))
    Eω, transform, FT = Luna.setup(grid, densityfun, responses, inputs, modes,
                                   :xy; full=false)
    linop = LinearOps.make_const_linop(grid, modes, λ0)
    statsfun = Stats.default(grid, Eω, modes, linop, transform; gas=gas)
    outvector45 = Output.MemoryOutput(0, grid.zmax, 2, statsfun)
    Luna.run(Eω, grid, linop, transform, FT, outvector45)
    Iωv45 = abs2.(outvector45.data["Eω"][:,1,1]) .+ abs2.(outvector45.data["Eω"][:,2,1])
    @test Iωs ≈ Iωv45
end


text:
@testset "utils/saveload.jl" begin

    xlen = 3
    zlen = 2
    keep = 2

    enc = ConditionalMvNormal(Dense(xlen, zlen))
    dec = ConditionalMvNormal(Dense(zlen, xlen))
    model = VAE(zlen, enc, dec)

    nowarn_logger = SimpleLogger(stdout, Logging.Error)
    model_dir   = mktempdir()
    @debug "  model_dir: $model_dir"

    @debug "  Testing `save_checkpoint`"
    model_ckpt = joinpath(model_dir, "ckpt.bson")
    history = MVHistory()
    push!(history, :loss, 1, 1)
    with_logger(nowarn_logger) do
        save_checkpoint(model_ckpt, model, history, keep=keep)
        save_checkpoint(model_ckpt, model, history, keep=keep)
        save_checkpoint(model_ckpt, model, history, keep=keep)
    end

    @test isfile(model_ckpt)
    @test length(readdir(dirname(model_ckpt))) == keep

    @debug "  Testing `load_checkpoint`"
    loaded_model, history = with_logger(nowarn_logger) do
        load_checkpoint(model_ckpt)
    end
    @test model.encoder.mapping.W == loaded_model.encoder.mapping.W

    opt = ADAM()
    lossf(x) = -elbo(model, x, β=1e-3)
    data = [(randn(Float32, xlen, 10),)]
    Flux.train!(lossf, params(model), data, opt)
    params_trained = get_params(model)

    with_logger(nowarn_logger) do
        save_checkpoint(model_ckpt, model, history, keep=keep)
    end

    loaded_model, history = with_logger(nowarn_logger) do
        load_checkpoint(model_ckpt)
    end

    @test length(params(model)) == length(params(loaded_model))
    @test !any(param_change(params_trained, loaded_model)) # did the params change?
    @test size(mean(loaded_model.encoder, randn(Float32, xlen))) == (zlen,)
end


text:
module TensorFlowDiffEq

using Reexport
@reexport using DiffEqBase

using TensorFlow, Compat
import DiffEqBase: solve, interpolation

# Abstract Types
@compat abstract type TensorFlowAlgorithm <: AbstractODEAlgorithm end

# ODE Algorithms
immutable odetf{Opt<:train.Optimizer} <: TensorFlowAlgorithm
  hl_widths::Vector{Int}
  optimizer::Opt
end

"""
    odetf
Specify the hyper-parameters used by the TensorFlow solver.
Including the optimisation hyper-paramters

for example
```julia
odetf(
    256,
    # use a neural network with 256 neurons in one-hiddlen layer in trail solution
    optimizer=train.MomentumOptimizer(0.5, 0.9)
    # Fit the network using a momentum based optimizser, learning rate of 0.5, momentum of 0.9
)
```

To make a deep network, pass a `Vector` to `hl_width`, given the width of each hiddden layer
Eg a 3 layer network
```julia
odetf(hl_width = [64, 128, 64])
```
"""
function odetf(;hl_width=64, optimizer=train.AdamOptimizer())
    odetf(hl_width, optimizer)
end

odetf(hl_width::Integer, optimizer) = odetf([hl_width], optimizer)

export odetf

immutable TensorFlowInterpolation{uElType}
    sess::Session
    u::Tensor{uElType}
end

grad_node(interp, deriv::Type{Val{0}}) = interp.u

function grad_node{N}(interp, deriv::Type{Val{N}})
    N >= 0 || throw(DomainError())
    #GOLDPLATE: this code could be enhanced, a lot. Right now it creates new nodes in the graph every time it is called. Which is a minor waste. The fix makes the code pretty hard to read and so I do not think it is worth it.
    as_default(interp.sess.graph) do
        prev = grad_node(interp, Val{N-1})
        tt = interp.sess.graph["tt"] #load the `tt` node from the graph
        grads(prev, tt)
    end
end

function (id::TensorFlowInterpolation){N}(tvals, idxs, deriv::Type{Val{N}})
    gn = grad_node(id, deriv)
    vals = [run(id.sess, gn, Dict(id.sess.graph["t"]=>[t]))' for t in tvals]
    #PREM-OPT: the indexing could be moved inside the network, to avoid even calculating gradients for columns are are not using, by first slicing `interp.u` inside the `grad_node` function.
    idxs == nothing ? _ret = vals : _ret = [val[:,idxs] for val in vals]
    typeof(idxs) <: Number || (idxs == nothing && length(first(_ret)) == 1) ? (ret = [first(_r) for _r in _ret]) : (ret = _ret)
    ret
end

(id::TensorFlowInterpolation){N}(tval::Number, idxs, deriv::Type{Val{N}}) = first(id([tval], idxs, deriv))

function (id::TensorFlowInterpolation){N}(v, tvals, idxs, deriv::Type{Val{N}})
    # In-place version, noting that truely inplace operations between julia and tensorflow are actually impossible
    v[:] = @view id(tvals, idxs, deriv)[:]
end

"Calculate the gradient per column of us, with regards to ts"
function grads(us, ts)
    outdim = get_shape(us, 2)
    hcat(map(1:outdim) do u_ii
        gradients(us[:,u_ii], ts)
    end...)
end

function solve(
    prob::AbstractODEProblem,
    alg::TensorFlowAlgorithm,
    timeseries = [], ts = [], ks = [];
    verbose=true, dt = nothing, maxiters = Int(1e4),
    progress_steps = 100, dense = true,
    timeseries_errors = true,
    dense_errors = false,
    kwargs...)

    u0 = prob.u0
    tspan = prob.tspan

    if dt == nothing
      error("dt must be set.")
    end

    sess=Session(Graph())
    f = prob.f

    uElType = eltype(u0)
    tType = typeof(tspan[1])
    t0 = tspan[1] # t0 is a constant as part of the u(t)
                  # so that even if interpolating to a different domain t0 remains from the training domain.
    outdim = length(u0)
    t_obs = collect(tspan[1]:tType(dt):tspan[2])

    @tf begin #Automatically name nodes based on RHS
        t = placeholder(tType, shape=[-1])
        tt = expand_dims(t, 2) #make it a matrix

        # u_trial trial network definition
        z = tt
        for (ii, hl_width) in enumerate(alg.hl_widths)
            width_below = get_shape(z, 2)
            w = get_variable("w_$ii", [width_below, hl_width], uElType)
            b = get_variable("b_$ii", [hl_width], uElType)
            z = nn.sigmoid(z*w + b; name="z_$ii")
        end
        width_below = get_shape(z, 2)
        w_out = get_variable([width_below, outdim], uElType)

        u = u0 + (tt .- t0).*z*w_out
        du_dt = grads(u, tt)
        deq_rhs = f(tt,u)

        loss = reduce_mean((du_dt - deq_rhs).^2)
        opt = train.minimize(alg.optimizer, loss)
    end

    run(sess, global_variables_initializer())

    if (verbose)
        @show run(sess, size(u), Dict(t=>t_obs))
        @show run(sess, size(du_dt), Dict(t=>t_obs))
        @show run(sess, size(deq_rhs), Dict(t=>t_obs))
    end

    for ii in 1:maxiters
        _, loss_o = run(sess, [opt, loss], Dict(t=>t_obs))
        if verbose && ii%progress_steps == 1
            println(loss_o)
        end
    end

    verbose && println("get final estimates with u_net(t)")
    u_net, du_dt_net = run(sess, [u, du_dt], Dict(t=>t_obs))
    verbose && println("preparing results")

    if typeof(u0) <: AbstractArray
        timeseries = Vector{typeof(u0)}(0)
        if dense
            du = Vector{typeof(u0)}(0)
        end
    for i=1:size(u_net, 1)
        push!(timeseries, reshape(view(u_net, i, :)', size(u0)))
        if dense
            push!(du, reshape(view(du_dt_net, i, :)', size(u0)))
        end
    end
    else
        timeseries = vec(u_net)
        if dense
          du = vec(du_dt_net)
        end
    end

    if !dense
      du = []
    end

    build_solution(prob,alg,t_obs,timeseries,
               dense = dense, du = du,
               interp = TensorFlowInterpolation(sess, u),
               timeseries_errors = timeseries_errors,
               dense_errors = dense_errors,
               retcode = :Success)

  end

end # module


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

"""
Either `GridSagLinear` or `GridSagBicubic` - determines the interpolation between sample points in the grid for a [`GridSagSurface`](@ref).
"""
@enum GridSagInterpolation GridSagLinear GridSagBicubic
export GridSagInterpolation, GridSagLinear, GridSagBicubic

"""
    GridSagSurface{T,N,S<:Union{ZernikeSurface{T,N},ChebyshevSurface{T,N}},Nu,Nv} <: ParametricSurface{T,N}

Either a Zernike (circular) or Chebyshev (rectangular) surface with grid sag height added to the base sag.
The surface shape is determined by either a linear or a bicubic spline interpolation of the `Nu×Nv` grid of sag values, set by the `interpolation` argument taking either `GridSagLinear` or `GridSagBicubic`.

Each entry in the grid is a vector of the form ``[z, \\frac{\\partial z}{\\partial x}, \\frac{\\partial z}{\\partial y}, \\frac{\\partial^2 z}{\\partial x \\partial y}]``.
The first data item corresponds to the lower left corner of the surface, that is, the corner defined by the -u and -v limit.
Each point that follows is read across the face of the surface from left to right moving upwards.
If zero is given for the partials (and using bicubic interpolation) then the partials will be approximated using finite differences.

The sag grid can be decentered from the surface in uv space, if so the surface may become wild outside of the area over which the grid is defined.
It is advised to clip the surface to the valid area using CSG operations in this case.

A surface can also be generated from a `.GRD` file by passing in the filename as the first and only positional argument. In this case the surface will be rectangular with optional radius and conic.

See docs for [`ZernikeSurface`](@ref) and [`ChebyshevSurface`](@ref) for details of the base surface.

```julia
GridSagSurface(basesurface::Union{ZernikeSurface{T,N},ChebyshevSurface{T,N}}, sag_grid::AbstractArray{T,3}; interpolation = GridSagBicubic, decenteruv = (0, 0))
GridSagSurface{T}(filename::String; radius = Inf, conic = 0, interpolation = GridSagBicubic)
```
"""
struct GridSagSurface{T,N,S<:Union{ZernikeSurface{T,N},ChebyshevSurface{T,N}},Nu,Nv} <: ParametricSurface{T,N}
    basesurface::S
    gridsag::SVector{Nu,SVector{Nv,SVector{4,T}}}
    interpolation::GridSagInterpolation
    decenteruv::Tuple{T,T}

    function GridSagSurface(basesurface::S, sag_grid::AbstractArray{T,2}; interpolation::GridSagInterpolation = GridSagBicubic, decenteruv::Tuple{T,T} = (zero(T), zero(T))) where {T<:Real,N,S<:Union{ZernikeSurface{T,N},ChebyshevSurface{T,N}}}
        grid = []
        Nv, Nu, d = size(sag_grid)
        @assert d == 4
        for r in 1:Nu
            row = []
            for c in 1:Nv
                push!(row, SVector{4,T}(sag_grid[c, r], 0, 0, 0))
            end
            push!(grid, row)
        end
        grid = SVector{Nu,SVector{Nv,SVector{4,T}}}(grid)
        return GridSagSurface(basesurface, grid, interpolation = interpolation, decenteruv = decenteruv)
    end

    function GridSagSurface(basesurface::S, sag_grid::AbstractArray{T,3}; interpolation::GridSagInterpolation = GridSagBicubic, decenteruv::Tuple{T,T} = (zero(T), zero(T))) where {T<:Real,N,S<:Union{ZernikeSurface{T,N},ChebyshevSurface{T,N}}}
        grid = []
        Nv, Nu, d = size(sag_grid)
        @assert d == 4
        for r in 1:Nu
            row = []
            for c in 1:Nv
                push!(row, SVector{4,T}(sag_grid[c, r, :]))
            end
            push!(grid, row)
        end
        grid = SVector{Nu,SVector{Nv,SVector{4,T}}}(grid)
        return GridSagSurface(basesurface, grid, interpolation = interpolation, decenteruv = decenteruv)
    end

    function GridSagSurface(basesurface::S, sag_grid::AbstractArray{<:AbstractVector{T},2}; interpolation::GridSagInterpolation = GridSagBicubic, decenteruv::Tuple{T,T} = (zero(T), zero(T))) where {T<:Real,N,S<:Union{ZernikeSurface{T,N},ChebyshevSurface{T,N}}}
        grid = []
        Nv, Nu = size(sag_grid)
        @assert length(sag_grid[1, 1]) == 4
        for r in 1:Nu
            push!(grid, sag_grid[:, r])
        end
        grid = SVector{Nu,SVector{Nv,SVector{4,T}}}(grid)
        return GridSagSurface(basesurface, grid, interpolation = interpolation, decenteruv = decenteruv)
    end

    function GridSagSurface(basesurface::S, sag_grid::SVector{Nu,SVector{Nv,SVector{4,T}}}; interpolation::GridSagInterpolation = GridSagBicubic, decenteruv::Tuple{T,T} = (zero(T), zero(T))) where {T<:Real,N,S<:Union{ZernikeSurface{T,N},ChebyshevSurface{T,N}},Nu,Nv}
        # TODO can probably simplify this?
        if interpolation === GridSagLinear
            @warn "Linear interpolation doesn't ensure C1 continuity, may cause problems..."
        end
        Δx = 2 * halfsizeu(basesurface) / (Nu - 1)
        Δy = 2 * halfsizev(basesurface) / (Nv - 1)
        kx = 1 / 2Δx
        ky = 1 / 2Δy
        kxy = kx * ky
        dxdκ = 2 * halfsizeu(basesurface) / (Nu - 1)
        dydγ = 2 * halfsizev(basesurface) / (Nv - 1)
        grid = MVector{Nu,MVector{Nv,SVector{4,T}}}(sag_grid)
        # if no derivs are provided (and we need them, i.e. bicubic) then approximate them here
        calcderivs = interpolation === GridSagBicubic && all(all.(iszero.(getindex.(v, [2:4])) for v in sag_grid))
        for u in 1:Nu
            for v in 1:Nv
                if calcderivs
                    # use finite differences if missing derivatives
                    z = sag_grid[u][v][1]
                    # at the edges just assume the derivative continues unchanged from that between the adjacent and current points
                    # in this case we need to half Δx and/or Δy hence umod and vmod below
                    u₊₁ = u + 1
                    u₋₁ = u - 1
                    umod = 1.0
                    if u === Nu
                        u₊₁ = u
                        umod = 2
                    elseif u === 1
                        u₋₁ = 1
                        umod = 2
                    end
                    v₊₁ = v + 1
                    v₋₁ = v - 1
                    vmod = 1.0
                    if v === Nv
                        v₊₁ = v
                        vmod = 2
                    elseif v === 1
                        v₋₁ = 1
                        vmod = 2
                    end
                    dzdx = kx * (sag_grid[u₊₁][v][1] - sag_grid[u₋₁][v][1]) * umod
                    dzdy = ky * (sag_grid[u][v₊₁][1] - sag_grid[u][v₋₁][1]) * vmod
                    dzdxdy = kxy * ((sag_grid[u₊₁][v₊₁][1] - sag_grid[u₊₁][v₋₁][1]) - (sag_grid[u₋₁][v₊₁][1] - sag_grid[u₋₁][v₋₁][1])) * umod * vmod
                else
                    z, dzdx, dzdy, dzdxdy = sag_grid[u][v]
                end
                # transform the derivatives from xy to κγ (i.e [0,1] in each patch)
                dzdκ = dzdx * dxdκ
                dzdγ = dzdy * dydγ
                dzdκdγ = dzdxdy * dxdκ * dydγ
                grid[u][v] = SVector{4,T}(z, dzdκ, dzdγ, dzdκdγ)
            end
        end
        return new{T,N,S,Nu,Nv}(basesurface, SVector{Nu,SVector{Nv,SVector{4,T}}}(grid), interpolation, decenteruv)
    end

    function GridSagSurface{T}(filename::String; radius::T = typemax(T), conic::T = zero(T), interpolation::GridSagInterpolation = GridSagBicubic) where {T<:Real}
        @assert !isnan(radius) && !isnan(conic)
        @assert isfile(filename)
        if !isvalid(readuntil(filename, " "))
            f = open(filename, enc"UTF-16LE", "r")
        else
            f = open(filename, "r")
        end
        spec = readline(f)
        spec = replace(spec, "\ufeff" => "")
        specs = split(spec)
        Nv = parse(Int, specs[1])
        Nu = parse(Int, specs[2])
        Δx = parse(Float64, specs[3])
        Δy = parse(Float64, specs[4])
        units = parse(Int, specs[5])
        if units == 0
            unitmod = 1
        elseif units == 1
            unitmod = 10
        elseif units == 2
            unitmod = 25.4
        elseif units == 3
            unitmod = 1000
        else
            throw(ErrorException("Invalid units"))
        end
        halfsizeu = (Nu - 1) * Δx / 2 / unitmod
        halfsizev = (Nv - 1) * Δy / 2 / unitmod
        xdec = parse(Float64, specs[6]) / unitmod
        ydec = parse(Float64, specs[7]) / unitmod
        decenteruv = (xdec / halfsizeu, ydec / halfsizev)
        grid = Array{Float64,3}(undef, Nv, Nu, 4)
        thisrowcount = 1
        thisrownum = 1
        for l in readlines(f)
            vs = [parse(Float64, x) for x in split(l)]
            if vs[5] == 1
                z, dzdx, dzdy, d2zdxdy = zeros(4) # nodata
            else
                z, dzdx, dzdy, d2zdxdy = vs[1:4]
            end
            # flip the y direction as GRD goes from -x, +y right and down while we go from -x, -y right and up
            grid[Nv - thisrownum + 1, thisrowcount, :] = [z / unitmod, dzdx / unitmod, dzdy / unitmod, d2zdxdy / unitmod^2]
            if thisrowcount == Nu
                thisrownum += 1
                thisrowcount = 1
            else
                thisrowcount += 1
            end
        end
        close(f)
        return GridSagSurface(ChebyshevSurface(halfsizeu, halfsizev, nothing, radius = radius, conic = conic), grid, interpolation = interpolation, decenteruv = decenteruv)
    end
end
export GridSagSurface

uvrange(::Type{GridSagSurface{T,N,S,Nu,Nv}}) where {T<:Real,N,S<:Union{ZernikeSurface{T,N},ChebyshevSurface{T,N}},Nu,Nv} = uvrange(S)

@inline function gridsag(s::GridSagSurface{T,3,S,Nu,Nv}, ui::Int, vi::Int, i::Int = 0) where {T<:Real,S,Nu,Nv}
    @assert 0 <= ui < Nu
    @assert 0 <= vi < Nv
    @assert i <= 4
    if i <= 0
        return @inbounds s.gridsag[ui + 1][vi + 1]
    else
        return @inbounds s.gridsag[ui + 1][vi + 1][i]
    end
end

@inline function bicubicα(surf::GridSagSurface{T,3,S,Nu,Nv}, uli::Int, uui::Int, vli::Int, vui::Int) where {T<:Real,S,Nu,Nv}
    f00, fu00, fv00, fuv00 = gridsag(surf, uli, vli)
    f01, fu01, fv01, fuv01 = gridsag(surf, uli, vui)
    f11, fu11, fv11, fuv11 = gridsag(surf, uui, vui)
    f10, fu10, fv10, fuv10 = gridsag(surf, uui, vli)
    return SMatrix{4,4,T,16}(1, 0, -3, 2, 0, 0, 3, -2, 0, 1, -2, 1, 0, 0, -1, 1) * SMatrix{4,4,T,16}(f00, f10, fu00, fu10, f01, f11, fu01, fu11, fv00, fv10, fuv00, fuv10, fv01, fv11, fuv01, fuv11) * SMatrix{4,4,T,16}(1, 0, 0, 0, 0, 0, 1, 0, -3, 3, -2, -1, 2, -2, 1, 1)
end

@inline function bicubicζ(surf::GridSagSurface{T,3,S,Nu,Nv}, uli::Int, uui::Int, vli::Int, vui::Int, κ::T, γ::T) where {T<:Real,S,Nu,Nv}
    return (SMatrix{1,4,T,4}(1, κ, κ^2, κ^3) * bicubicα(surf, uli, uui, vli, vui) * SMatrix{4,1,T,4}(1, γ, γ^2, γ^3))[1]
end

@inline function bicubicdζ(surf::GridSagSurface{T,3,S,Nu,Nv}, uli::Int, uui::Int, vli::Int, vui::Int, κ::T, γ::T) where {T<:Real,S,Nu,Nv}
    α = bicubicα(surf, uli, uui, vli, vui)
    dκ = (SMatrix{1,4,T,4}(0, 1, 2κ, 3κ^2) * α * SMatrix{4,1,T,4}(1, γ, γ^2, γ^3))[1]
    dγ = (SMatrix{1,4,T,4}(1, κ, κ^2, κ^3) * α * SMatrix{4,1,T,4}(0, 1, 2γ, 3γ^2))[1]
    return dκ, dγ
end

@inline function normalizedcoords(Nu::Int, Nv::Int, u::T, v::T, decenter::Tuple{T,T}) where {T<:Real}
    u = u - decenter[1]
    v = v - decenter[2]
    û = (u + 1) / 2
    v̂ = (v + 1) / 2
    uli = max(min(floor(Int, û * (Nu - 1)), Nu - 2), 0)
    uui = min(uli + 1, Nu - 1)
    vli = max(min(floor(Int, v̂ * (Nv - 1)), Nv - 2), 0)
    vui = min(vli + 1, Nv - 1)
    κ = û * (Nu - 1) - uli
    γ = v̂ * (Nv - 1) - vli
    return uli, uui, vli, vui, κ, γ
end

@inline function ζ(surf::GridSagSurface{T,3,S,Nu,Nv}, u::T, v::T) where {T<:Real,S,Nu,Nv}
    uli, uui, vli, vui, κ, γ = normalizedcoords(Nu, Nv, u, v, surf.decenteruv)
    if surf.interpolation === GridSagLinear
        sagll = gridsag(surf, uli, vli, 1)
        saglu = gridsag(surf, uli, vui, 1)
        sagul = gridsag(surf, uui, vli, 1)
        saguu = gridsag(surf, uui, vui, 1)
        z = κ * (γ * saguu + (1 - γ) * sagul) + (1 - κ) * (γ * saglu + (1 - γ) * sagll)
    elseif surf.interpolation === GridSagBicubic
        z = bicubicζ(surf, uli, uui, vli, vui, κ, γ)
    end
    return z
end

@inline function dζ(surf::GridSagSurface{T,3,S,Nu,Nv}, u::T, v::T) where {T<:Real,S,Nu,Nv}
    uli, uui, vli, vui, κ, γ = normalizedcoords(Nu, Nv, u, v, surf.decenteruv)
    # calculate sag derivative at point
    if surf.interpolation === GridSagLinear
        sagll = gridsag(surf, uli, vli, 1)
        saglu = gridsag(surf, uli, vui, 1)
        sagul = gridsag(surf, uui, vli, 1)
        saguu = gridsag(surf, uui, vui, 1)
        dζdκ = (γ * saguu + (1 - γ) * sagul) - (γ * saglu + (1 - γ) * sagll)
        dζdγ = (κ * saguu + (1 - κ) * saglu) - (κ * sagul + (1 - κ) * sagll)
        # NOT C1 CONTINUOUS!!!
    elseif surf.interpolation === GridSagBicubic
        dζdκ, dζdγ = bicubicdζ(surf, uli, uui, vli, vui, κ, γ)
        # C1 continuous, but not C2 continuous (at grid boundaries)
    end
    # change variables
    dζdu = dζdκ * (Nu - 1) / 2
    dζdv = dζdγ * (Nv - 1) / 2
    return dζdu, dζdv
end

function point(surf::GridSagSurface{T,3,S}, ρ::T, ϕ::T)::SVector{3,T} where {T<:Real,S<:ZernikeSurface{T,3}}
    x, y, z = point(surf.basesurface, ρ, ϕ)
    u = ρ * cos(ϕ)
    v = ρ * sin(ϕ)
    return SVector{3,T}(x, y, z + ζ(surf, u, v))
end

function point(surf::GridSagSurface{T,3,S}, u::T, v::T)::SVector{3,T} where {T<:Real,S<:ChebyshevSurface{T,3}}
    x, y, z = point(surf.basesurface, u, v)
    return SVector{3,T}(x, y, z + ζ(surf, u, v))
end

function partials(surf::GridSagSurface{T,3,S}, ρ::T, ϕ::T)::Tuple{SVector{3,T},SVector{3,T}} where {T<:Real,S<:ZernikeSurface{T,3}}
    pρ, pϕ = partials(surf.basesurface, ρ, ϕ)
    u = ρ * cos(ϕ)
    v = ρ * sin(ϕ)
    dζdu, dζdv = dζ(surf, u, v)
    dζdρ = dζdu * cos(ϕ) + dζdv * sin(ϕ)
    dζdϕ = dζdv * u - dζdu * v
    return SVector{3,T}(pρ[1], pρ[2], pρ[3] + dζdρ), SVector{3,T}(pϕ[1], pϕ[2], pϕ[3] + dζdϕ)
end

function partials(surf::GridSagSurface{T,3,S}, u::T, v::T)::Tuple{SVector{3,T},SVector{3,T}} where {T<:Real,S<:ChebyshevSurface{T,3}}
    pu, pv = partials(surf.basesurface, u, v)
    du, dv = dζ(surf, u, v)
    return SVector{3,T}(pu[1], pu[2], pu[3] + du), SVector{3,T}(pv[1], pv[2], pv[3] + dv)
end

function normal(z::GridSagSurface{T,N}, u::T, v::T)::SVector{N,T} where {T<:Real,N}
    du, dv = partials(z, u, v)
    return normalize(cross(du, dv))
end

uv(s::GridSagSurface{T,3}, p::SVector{3,T}) where {T<:Real} = uv(s.basesurface, p)

function onsurface(surf::GridSagSurface{T,3,S}, p::SVector{3,T}) where {T<:Real,S}
    u, v = uv(surf, p)
    if abs(u) > one(T) || (abs(v) > one(T) && S <: ChebyshevSurface{T,3})
        return false
    else
        # not sure it's really necessary/correct to include the cylinder here
        surfpoint = point(surf, u, v)
        return @inbounds samepoint(p[3], surfpoint[3]) # || (ρ == 1 && p[3] < z)
    end
end

function inside(surf::GridSagSurface{T,3,S}, p::SVector{3,T}) where {T<:Real,S}
    u, v = uv(surf, p)
    if abs(u) > one(T) || (abs(v) > one(T) && S <: ChebyshevSurface{T,3})
        return false
    else
        surfpoint = point(surf, u, v)
        return @inbounds p[3] < surfpoint[3]
    end
end

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

# Assumes the ray has been transformed into the canonical coordinate frame which has the vertical axis passing through (0,0,0) and aligned with the z axis.
function surfaceintersection(surf::AcceleratedParametricSurface{T,3,GridSagSurface{T,3,S,Nu,Nv}}, r::AbstractRay{T,3}) where {T<:Real,S,Nu,Nv}
    boundingintvl = surfaceintersection(boundingobj(surf.surface.basesurface), r)
    if boundingintvl isa EmptyInterval{T}
        return EmptyInterval(T)
    else
        if doesintersect(surf.triangles_bbox, r) || inside(surf.triangles_bbox, origin(r))
            surfint = triangulatedintersection(surf, r)
            if !(surfint isa EmptyInterval{T})
                return intervalintersection(boundingintvl, surfint)
            end
        end
        # hasn't hit the surface
        if lower(boundingintvl) isa RayOrigin{T} && upper(boundingintvl) isa Infinity{T}
            if inside(surf.surface, origin(r))
                return Interval(RayOrigin(T), Infinity(T))
            else
                return EmptyInterval(T)
            end
            # otherwise check that the intersection is underneath the surface
        else
            p = point(closestintersection(boundingintvl, false))
            ρ, ϕ = uv(surf, p)
            surfpoint = point(surf.surface, ρ, ϕ)
            if @inbounds p[3] < surfpoint[3]
                return boundingintvl # TODO!! UV (and interface) issues?
            else
                return EmptyInterval(T)
            end
        end
    end
end

function AcceleratedParametricSurface(surf::GridSagSurface{T,N,Z}, numsamples::Int = 17; interface::NullOrFresnel{T} = NullInterface(T)) where {T<:Real,N,Z<:ZernikeSurface{T,N}}
    # Zernike uses ρ, ϕ uv space so need to modify extension of triangulation
    a = AcceleratedParametricSurface(surf, triangulate(surf, numsamples, true, false, true, false), interface = interface)
    emptytrianglepool!(T)
    return a
end

function BoundingBox(surf::GridSagSurface{T,3}) where {T<:Real}
    bbox = BoundingBox(surf.basesurface)
    maxsag = maximum(maximum.(getindex.(surf.gridsag[i], 1) for i in 1:length(surf.gridsag)))
    minsag = minimum(minimum.(getindex.(surf.gridsag[i], 1) for i in 1:length(surf.gridsag)))
    return BoundingBox(bbox.xmin, bbox.xmax, bbox.ymin, bbox.ymax, bbox.zmin + minsag, bbox.zmax + maxsag)
end


text:
# ------------ #
# Original API #
# ------------ #

# lindef.m -- DONE
function lindef(breaks::AbstractVector, evennum::Int=0)
    p = LinParams(breaks, evennum)
    return (
        length(p.breaks),
        minimum(p.breaks),
        maximum(p.breaks),
        Any[p.breaks, evennum],
        p
    )
end

# linnode.m -- DONE
linnode(breaks, evennum) = nodes(LinParams(breaks, evennum))

# lindop.m
function lindop(breaks, evennum, order)
    D, params = derivative_op(LinParams(breaks, evennum), order)
    n, a, b = length(params.breaks), params.breaks[1], params.breaks[end]
    D, n, a, b, params
end

# linbase.m -- DONE
function linbase(breaks, evennum=0, x=breaks, order=0)
    B = evalbase(LinParams(breaks, evennum), x, order)
    B, x
end


text:
using Interpolations

export interpolate_data

# """
#     interpolate_data(ntwk::NetworkData{T}, freq) where {T<:NetworkParams}
# Interpolate `NetworkData` at a new freqeuncy range `freq`.
# """
for p in (:Sparams, :Yparams, :Zparams, :ABCDparams)
    @eval interpolate_data(ntwk::NetworkData{S, ($p){T}},
        freq::AbstractVector{R}) where {S<:Real, T<:Real, R<:Real} = begin
            nPort = ntwk.nPort
            d = []
            U = promote_type(T, S)
            for i in 1:nPort
                e = []
                for j in 1:nPort
                    Sij_itp = interpolate((U.(ntwk.frequency), ),
                        Complex{U}.(ntwk[(i, j), :]), Gridded(Linear()))
            # itp = interpolate(ntwk[(i, j), :], BSpline(Cubic(Line())), OnGrid())

            # frng = ntwk.frequency[1]:(ntwk.frequency[end]-ntwk.frequency[1])/(ntwk.nPoint-1):ntwk.frequency[end]
            # sitp = scale(itp, frng)
#              sitp = scale(itp, ntwk.frequency)
                    push!(e, Sij_itp[freq])
                    # push!(e, [sitp[f] for f in freq])
                end
                push!(d, e)
            end
            params = [($p)([d[i][j][n] for i in 1:nPort, j in 1:nPort]) for n in 1:length(freq)]
            NetworkData(ntwk.ports, freq, params)
        end
end

# """
#     interpolate_data(cdata::CircuitData{T}, freq) where {T<:CircuitParams}
# Interpolate `CircuitData` at a new freqeuncy range `freq`.
# """
for p in (:Impedance, :Admittance)
    @eval interpolate_data(cdata::CircuitData{S, ($p){T}},
        freq::AbstractVector{R}) where {S<:Real, T<:Real, R<:Real} = begin
            U = promote_type(T, S)
            itp = interpolate((U.(cdata.frequency), ), Complex{U}.(cdata[:]), Gridded(Linear()))
            CircuitData(freq, ($p)(itp[freq]))
        end
end


text:
# order 2
h = [5//12, 7//6, 11//12, 1, 1]
e = [1//1, 0, 0, 0, 0]
H = Diagonal(h)

q1 = [0 7//12 -1//12 0 0]
q2 = [-q1[2] 0 7//12 0 0]
q3 = [-q1[3] -q2[3] 0 1//2 0]
q4 = [0 0 -1//2 0 1//2]
q5 = [0 0 0 -1//2 0]
Q = vcat(q1, q2, q3, q4, q5)

D1 = H \ (Q - 1//2*e*e'); display(D1)

x=1//1:length(h)
ks = 0:5
for k in ks
    u=x.^k
    println("k = ", k)
    println(D1*u - k*x.^(k-1))
end

# order 4
h = [511//1600, 3971//2880, 929//1440, 587//480, 2659//2880, 14551//14400, 1, 1, 1]
e = [1//1, 0, 0, 0, 0, 0, 0, 0, 0]
H = Diagonal(h)

q1 = [0 436061//680400 -312391//5443200 -100109//907200 54107//5443200 263//15552 0 0 0]
q2 = [-q1[2] 0 238831//544320 290357//1088640 -125//5184 -220357//5443200 0 0 0]
q3 = [-q1[3] -q2[3] 0 3467//7776 -74827//1088640 3859//907200 0 0 0]
q4 = [-q1[4] -q2[4] -q3[4] 0 342799//544320 -149959//5443200 0 0 0]
q5 = [-q1[5] -q2[5] -q3[5] -q4[5] 0 428789//680400 -1//12 0 0]
q6 = [-q1[6] -q2[6] -q3[6] -q4[6] -q5[6] 0 2//3 -1//12 0]
q7 = [0 0 0 0 1//12 -2//3 0 2//3 -1//12]
q8 = [0 0 0 0 0 1//12 -2//3 0 2//3]
q9 = [0 0 0 0 0 0 1//12 -2//3 0]
Q = vcat(q1, q2, q3, q4, q5, q6, q7, q8, q9)

D1 = H \ (Q - 1//2*e*e'); display(D1)

x=1//1:length(h)
ks = 0:5
for k in ks
    u=x.^k
    println("k = ", k)
    println(D1*u - k*x.^(k-1))
end

# order 6
h = [7497391//25401600, 1106227//725760, 105317//403200, 260179//145152, 303631//725760,
    513973//403200, 671171//725760, 25631999//25401600, 1, 1, 1, 1]
e = [1//1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
H = Diagonal(h)

q1 = [0 14085349//21168000 -205189//8064000 -7583609//38102400 -89573//4233600 9//80 -39279943//1524096000 -13//2000 0 0 0 0]
q2 = [-q1[2] 0 56626691//304819200 97579603//152409600 1214239//16934400 -246914711//762048000 23380423//304819200 3//200 0 0 0 0]
q3 = [-q1[3] -q2[3] 0 22302157//101606400 -6048877//60963840 9737351//169344000 -508309//25401600 780937//304819200 0 0 0 0]
q4 = [-q1[4] -q2[4] -q3[4] 0 8120941//30481920 3467971//6350400 -42123859//304819200 -2075677//152409600 0 0 0 0]
q5 = [-q1[5] -q2[5] -q3[5] -q4[5] 0 29205107//152409600 15263287//304819200 -81127//3386880 0 0 0 0]
q6 = [-q1[6] -q2[6] -q3[6] -q4[6] -q5[6] 0 321139631//508032000 -49594423//762048000 1//60 0 0 0]
q7 = [-q1[7] -q2[7] -q3[7] -q4[7] -q5[7] -q6[7] 0 30841499//43545600 -3//20 1//60 0 0]
q8 = [-q1[8] -q2[8] -q3[8] -q4[8] -q5[8] -q6[8] -q7[8] 0 3//4 -3//20 1//60 0]
q9 = [0 0 0 0 0 -1//60 3//20 -3//4 0 3//4 -3//20 1//60]
q10 = [0 0 0 0 0 0 -1//60 3//20 -3//4 0 3//4 -3//20]
q11 = [0 0 0 0 0 0 0 -1//60 3//20 -3//4 0 3//4]
q12 = [0 0 0 0 0 0 0 0 -1//60 3//20 -3//4 0]
Q = vcat(q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12)

D1 = H \ (Q - 1//2*e*e'); display(D1)

x=1//1:length(h)
ks = 0:7
for k in ks
    u=x.^k
    println("k = ", k)
    println(D1*u - k*x.^(k-1))
end


text: Dict{String,Any}("kilo" => Dict{String,Any}("value" => "1000","type" => "integer"))

max_stars_repo_name: SmiteWindows/julia

max_issues_repo_name: SmiteWindows/julia

max_forks_repo_name: SmiteWindows/julia


text:
using TensorCore
using LinearAlgebra
using Test

@testset "Ambiguities" begin
    if VERSION >= v"1.6-"
        @test isempty(detect_ambiguities(TensorCore))
    else
        @test isempty(detect_ambiguities(TensorCore, Base, Core, LinearAlgebra))
    end
end

@testset "tensor and hadamard" begin
    for T in (Int, Float32, Float64, BigFloat)
        a = [T[1, 2], T[-3, 7]]
        b = [T[5, 11], T[-13, 17]]
        @test map(⋅, a, b) == map(dot, a, b) == [27, 158]
        @test map(⊙, a, b) == map(hadamard, a, b) == [a[1].*b[1], a[2].*b[2]]
        @test map(⊗, a, b) == map(tensor, a, b) == [a[1]*transpose(b[1]), a[2]*transpose(b[2])]
        @test hadamard!(fill(typemax(Int), 2), T[1, 2], T[-3, 7]) == [-3, 14]
        @test tensor!(fill(typemax(Int), 2, 2), T[1, 2], T[-3, 7]) == [-3 7; -6 14]
    end

    @test_throws DimensionMismatch [1,2] ⊙ [3]
    @test_throws DimensionMismatch hadamard!([0, 0, 0], [1,2], [-3,7])
    @test_throws DimensionMismatch hadamard!([0, 0], [1,2], [-3])
    @test_throws DimensionMismatch hadamard!([0, 0], [1], [-3,7])
    @test_throws DimensionMismatch tensor!(Matrix{Int}(undef, 2, 2), [1], [-3,7])
    @test_throws DimensionMismatch tensor!(Matrix{Int}(undef, 2, 2), [1,2], [-3])

    u, v = [2+2im, 3+5im], [1-3im, 7+3im]
    @test u ⋅ v == conj(u[1])*v[1] + conj(u[2])*v[2]
    @test u ⊙ v == [u[1]*v[1], u[2]*v[2]]
    @test u ⊗ v == [u[1]*v[1] u[1]*v[2]; u[2]*v[1] u[2]*v[2]]
    @test hadamard(u, v) == u ⊙ v
    @test tensor(u, v)   == u ⊗ v
    dest = similar(u)
    @test hadamard!(dest, u, v) == u ⊙ v
    dest = Matrix{Complex{Int}}(undef, 2, 2)
    @test tensor!(dest, u, v) == u ⊗ v

    for (A, B, b) in (([1 2; 3 4], [5 6; 7 8], [5,6]),
                      ([1+0.8im 2+0.7im; 3+0.6im 4+0.5im],
                       [5+0.4im 6+0.3im; 7+0.2im 8+0.1im],
                       [5+0.6im,6+0.3im]))
        @test A ⊗ b == cat(A*b[1], A*b[2]; dims=3)
        @test A ⊗ B == cat(cat(A*B[1,1], A*B[2,1]; dims=3),
                           cat(A*B[1,2], A*B[2,2]; dims=3); dims=4)
    end

    A, B = reshape(1:27, 3, 3, 3), reshape(1:4, 2, 2)
    @test A ⊗ B == [a*b for a in A, b in B]

    # Adjoint/transpose is a dual vector, not an AbstractMatrix
    v = [1,2]
    @test_throws ErrorException v ⊗ v'
    @test_throws ErrorException v ⊗ transpose(v)
    @test_throws ErrorException v' ⊗ v
    @test_throws ErrorException transpose(v) ⊗ v
    @test_throws ErrorException v' ⊗ v'
    @test_throws ErrorException transpose(v) ⊗ transpose(v)
    @test_throws ErrorException v' ⊗ transpose(v)
    @test_throws ErrorException transpose(v) ⊗ v'
    @test_throws ErrorException A ⊗ v'
    @test_throws ErrorException A ⊗ transpose(v)

    # Docs comparison to `kron`
    v, w = [1,2,3], [5,7]
    @test kron(v,w) == vec(w ⊗ v)
    @test w ⊗ v == reshape(kron(v,w), (length(w), length(v)))
end

@testset "arrays of arrays" begin
    @testset "pseudo-complex" begin

        # This is a matrix representation of the complex numbers:
        jm = [0 1; -1 0]
        om = [1 0; 0 1]
        @test jm * jm == -om
        @test jm * om == om * jm == jm
        @test om * om == om

        @test jm' == -jm

        mreal(x::Matrix) = tr(x)/2
        jmag(x::Matrix) = tr(-x * jm)/2
        @test mreal((2om+3jm)^2) == real((2+3im)^2)
        @test jmag((2om+3jm)^2) == imag((2+3im)^2)

        v = [im, 2, 3+im]
        vm = [jm, 2om, 3om + jm]

        w = [4 + im, 5, 6+7im]
        wm = [4om + jm, 5om, 6om+7jm]
        @test real.(w) == mreal.(wm)
        @test imag.(w) == jmag.(wm)

        # hadamard & tensor
        @test real.(v ⊙ w) == mreal.(vm ⊙ wm)
        @test imag.(v ⊙ w) == jmag.(vm ⊙ wm)

        @test real.(v ⊗ w) == mreal.(vm ⊗ wm)
        @test imag.(v ⊗ w) == jmag.(vm ⊗ wm)

        # kron
        @test mreal.(kron(vm, wm)) == real.(kron(v, w))
        @test jmag.(kron(vm, wm)) == imag.(kron(v, w))

        @test vec(kron(v', w')') == kron(v, w)
        @test vec(kron(vm', wm')') == kron(vm, wm)

        # times, cross & dot
        outervw = v * w'
        @test outervw == mreal.(vm * wm') + im * jmag.(vm * wm')

        vcrossw = v × w
        @test vcrossw == mreal.(vm × wm) + im * jmag.(vm × wm)

        dvw = dot(v, w) # 36 + 11im
        @test dvw == v' * w
        @test dvw == sum(v[i]' * w[i] for i in 1:3)

        m_star = vm' * wm  # recursive adjoint is important here!
        if !(VERSION < v"1.5-")
            @test dvw == mreal(m_star) + im * jmag(m_star) # fails on 1.4
        end

        m_wrong = first(permutedims(vm) * wm)
        @test mreal(m_wrong) + im * jmag(m_wrong) == 20 + 31im

        m_sum = sum(vm[i]' * wm[i] for i in 1:3)
        @test dvw == mreal(m_sum) + im * jmag(m_sum)

        @test_broken dvw == dot(vm, wm) # 72, wtf?
        sum(dot(vm[i], wm[i]) for i in 1:3)
        # I think this behaviour was introduced here:
        # https://github.com/JuliaLang/julia/pull/27401
        # I don't actually see an argument for this, and there also seem to be no tests
        # of this in LinearAlgebra (although I have not changed & run them to check).

        # dot with kron
        kdot = kron(v, v)' * kron(w, w)
        kmdot = kron(vm, vm)' * kron(wm, wm)
        if !(VERSION < v"1.5-")
            @test kdot == mreal(kmdot) + im * jmag(kmdot)
        end

    end
    @testset "block matrix & vector" begin

        # A matrix of matrices, and a vector of vectors:
        M = [rand(Int8,2,2).+3im for i in 1:3, j in 1:4]
        N = [rand(Int8,2,2).+5im for i in 1:3, j in 1:4]

        V = [rand(Int8,2).+7im for i in 1:3]
        W = [rand(Int8,2).+11im for j in 1:4]

        # ordinary LinearAlgebra
        @test M * N' == (N * M')'

        @test M * W == (W' * M')' # recursive adjoint is important here!
        @test V' * M == (M' * V)'

        @test V' * V == dot(V,V) # here the behaviour of dot is not wrong
        if VERSION >= v"1.4"
            @test V' * M * W == dot(V,M,W)
        end

        # kron & TensorCore
        @test M ⊙ M == M .* M != map(⊙, M, M)

        @test first(M ⊗ N) == first(M) * first(N)
        @test size(M ⊗ N) == (3,4,3,4)

        @test_throws MethodError kron(W, W) # what should this do, if anything?
        @test_throws MethodError W ⊗ V

    end

end

@testset "boxdot" begin

    # Matrices and vectors
    A = [1 2+im; 3 4im]
    B = [5im 6; 7+im 8]
    c = [1, 2+3im]
    d = [3im, 4-5im]

    @test A ⊡ B == A * B
    @test A ⊡ c == A ⊡ c
    @test c ⊡ A == vec(transpose(c) * A)
    @test c ⊡ d == sum(c .* d)

    @test A' ⊡ B == A' * B
    @test A ⊡ B' == A * B'
    @test A' ⊡ B' == A' * B'

    # Dual vectors
    @test c' ⊡ d == dot(c, d)
    @test c ⊡ d' == conj(dot(c, d))
    @test c' ⊡ d' == conj(dot(c', d))

    @test transpose(c) ⊡ d == sum(c .* d)
    @test c ⊡ transpose(d) == sum(c .* d)
    @test transpose(c) ⊡ transpose(d) == sum(c .* d)
    @test transpose(c) ⊡ adjoint(d) == sum(c .* conj(d))
    @test adjoint(c) ⊡ transpose(d) == dot(c,d)

    @test A ⊡ c' isa Adjoint
    @test A ⊡ c' == (c ⊡ A')'
    @test c' ⊡ A isa Adjoint
    @test c' ⊡ A == (A' * c)'

    @test B' ⊡ c' isa Adjoint
    @test B' ⊡ c' == (c ⊡ B)'
    @test c' ⊡ B' isa Adjoint
    @test c' ⊡ B' == (B * c)'

    @test B' ⊡ transpose(c) isa Transpose
    @test B ⊡ transpose(c) == transpose(c ⊡ transpose(B))
    @test transpose(c) ⊡ B' isa Transpose
    @test transpose(c) ⊡ B == transpose(transpose(B) * c)

    # Higher-dimensional arrays
    E3 = cat(A, -B, dims=3)
    F4 = cat(E3, conj(E3 .+ 1), dims=4)
    E3adjoint = conj(permutedims(E3, (3,2,1)))
    E3lazy = PermutedDimsArray(permutedims(E3, (3,2,1)), (3,2,1))
    F4lazy = PermutedDimsArray(permutedims(F4, (4,3,2,1)), (4,3,2,1))
    @test E3lazy == E3
    @test F4lazy == F4

    @test E3 ⊡ A == reshape(reshape(E3, 4,2) * A, 2,2,2)
    @test size(A ⊡ E3) == (2, 2, 2)
    @test size(B ⊡ F4) == (2, 2, 2, 2)
    @test size(E3 ⊡ F4) == (2, 2, 2, 2, 2)
    @test E3lazy ⊡ A == E3 ⊡ A
    @test E3lazy ⊡ F4lazy == E3 ⊡ F4
    @test F4lazy ⊡ E3lazy == F4 ⊡ E3

    @test c ⊡ E3 == reshape(transpose(c) * reshape(E3, 2,4), 2,2)
    @test size(F4 ⊡ c) == (2, 2, 2)

    @test c' ⊡ E3 isa Matrix
    @test c' ⊡ E3 == conj(c) ⊡ E3
    @test c' ⊡ E3 == (E3adjoint ⊡ c)'
    @test transpose(c) ⊡ E3 == c ⊡ E3
    @test c ⊡ E3lazy == c ⊡ E3

    @test E3 ⊡ c' isa Matrix
    @test E3 ⊡ c' == E3 ⊡ conj(c)
    @test E3 ⊡ c' == (c ⊡ E3adjoint)'
    @test E3 ⊡ transpose(c) == E3 ⊡ c
    @test E3 ⊡ transpose(c) == transpose(c ⊡ PermutedDimsArray(E3, (3,2,1)))

    # Errors
    @test_throws DimensionMismatch ones(3) ⊡ ones(4)
    @test_throws DimensionMismatch ones(3) ⊡ ones(4)'
    @test_throws DimensionMismatch ones(3)' ⊡ ones(4)
    @test_throws DimensionMismatch ones(3)' ⊡ ones(4)'
    @test_throws DimensionMismatch ones(2,3) ⊡ ones(4)
    @test_throws DimensionMismatch ones(2,3) ⊡ ones(4)'
    @test_throws DimensionMismatch ones(2,3) ⊡ ones(4,5)

    # In-place
    @test boxdot!(similar(c), A, c) == A * c
    if VERSION >= v"1.3"
        @test boxdot!(similar(c), A, c, 100) == A * c * 100
        @test boxdot!(copy(c), B, d, 100, -5) == B * d * 100 .- 5 .* c
    end

    @test boxdot!(similar(c), A, c') == A * conj(c)
    @test boxdot!(similar(c,1), c, d') == [sum(c .* conj(d))]

    @test boxdot!(similar(c)', c', A) == c' * A
    @test boxdot!(similar(c,1,2), c', A) == c' * A

    @test boxdot!(similar(c,1), c', d) == [dot(c, d)]
end

@testset "_adjoint" begin
    A = [1 2+im; 3 4im]
    E3 = cat(A, -A, dims=3)
    F4 = cat(E3, conj(E3 .+ 1), dims=4)
    M3 = [rand(ComplexF32, 2,2) for _ in 1:3, _ in 1:4, _ in 1:5];

    @test TensorCore._adjoint(E3) == conj(permutedims(E3, (3,2,1)))
    @test TensorCore._transpose(F4) == permutedims(F4, (4,3,2,1))
    @test TensorCore._adjoint(A) == A'
    @test TensorCore._transpose(A) == transpose(A)

    @test TensorCore._adjoint(M3)[5,4,3] == adjoint(M3[3,4,5])
    @test TensorCore._transpose(M3)[1,2,3] == transpose(M3[3,2,1])

    @test TensorCore._adjoint(E3 ⊗ F4) == TensorCore._adjoint(F4) ⊗ TensorCore._adjoint(E3)
    @test TensorCore._adjoint(E3 ⊡ F4) == TensorCore._adjoint(F4) ⊡ TensorCore._adjoint(E3)
    @test TensorCore._adjoint(M3 ⊡ TensorCore._adjoint(M3)) ≈ M3 ⊡ TensorCore._adjoint(M3)
end


text:
# Copyright 2019 Tobias Frilling
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

"""
    readUMX(filename::String, directory=pwd())

Read the contents of a `*.umx` and return a `Matrix{Float64}`
"""
function readUMX(filename::String, directory = pwd())
    filename = prepare_path(filename, "umx", directory)
    result = Matrix{Float64}(undef, (0,0))

    open(filename, "r") do f
        skipStarting(f, ['#', '%'])
        result = readdlm(f, Float64, skipblanks = true)
    end

    return result
end

max_stars_repo_name: ckafi/DataIO

max_issues_repo_name: ckafi/DataIO

max_forks_repo_name: ckafi/DataIO


text:
function globalids(sig, slots)
  used = Set{Int}()
  ret = Int[]
  for T in sig
    i = findfirst(i -> slots[i] == T && !(i in used), 1:length(slots))
    i == nothing && (i = length(push!(slots, T)))
    push!(used, i)
    push!(ret, i-1)
  end
  return ret
end

function multivalue!(mod, f, xs, slots)
  for (i, x) in enumerate(xs)
    if x isa Return && length(f.result) > 1
      splice!(xs, i:0, reverse(SetGlobal.(globalids(f.result, slots))))
      return
    elseif x isa Call
      rs = func(mod, x.name).result
      length(rs) > 1 || continue
      splice!(xs, i+1:0, GetGlobal.(globalids(rs, slots)))
    elseif x isa Union{Block,If,Loop}
      multivalue_shim!(mod, f, x, slots)
    end
  end
end

multivalue_shim!(mod, f, c::Block, slots) = multivalue!(mod, f, c.body, slots)
multivalue_shim!(mod, f, c::Loop, slots)  = multivalue!(mod, f, c.body, slots)

multivalue_shim!(mod, f, c::If, slots) =
  (multivalue!(mod, f, c.t, slots); multivalue!(mod, f, c.t, slots))

function multivalue_shim!(mod::Module)
  slots = WType[]
  for f in mod.funcs
    multivalue_shim!(mod, f, f.body, slots)
  end
  for f in mod.funcs
    length(f.result) > 1 && empty!(f.result)
  end
  append!(mod.globals, Global.(slots))
  return mod
end


text:
#!/usr/bin/env julia
using MFCC: mfcc
using LinearAlgebra: norm
using DataStructures
using Clustering
using ArgParse
using WAV
using Crayons

struct SeamPoint
    startidx::Int
    endidx::Int
    loss::Float32
end

import Base.isless
isless(a::SeamPoint, b::SeamPoint) = a.loss < b.loss

function restride(input::AbstractArray{T, 2}) where T
    out = Vector{T}(undef, length(input))
    m, n = size(input)

    for idx = 1:m, jdx = 1:n
        @inbounds out[(idx-1) * n + jdx] = input[idx, jdx]
    end

    return out
end

# the number of packs per one second
const s = 32

# the number of MFC coefficients
const numcep = 30

# length of the pack (measured in samples)
const τ = 1

function extract(cepstrum::Vector{Float32}, lowerbound::Int, upperbound::Int, numpacks::Int)::Vector{SeamPoint}
    capacity = 128
    heap = BinaryMaxHeap{SeamPoint}([SeamPoint(1, 1, Inf)])

    for a = 1:numpacks - lowerbound - τ, b = a + lowerbound : min(a + upperbound - τ, numpacks - τ)
        pack₁ = @view cepstrum[(a - 1) * numcep + 1 : (a + τ) * numcep]
        pack₂ = @view cepstrum[(b - 1) * numcep + 1 : (b + τ) * numcep]

        loss = norm(pack₁ - pack₂, 1)

        if loss < top(heap).loss
            length(heap) > capacity && pop!(heap)

            push!(heap, SeamPoint(a, b, loss))
        end
    end

    return extract_all!(heap)
end

function main(args)
    rules = ArgParseSettings("Extract precise loops from the audio track")

    @add_arg_table! rules begin
        "audio-path"
        arg_type = String
        required = true
        help = "path to the source audio file"

        "--min-seconds", "--min"
        nargs = '?'
        arg_type = Int
        default = 10
        help = "minimum length of the loop"

        "--max-seconds", "--max"
        nargs = '?'
        arg_type = Int
        default = 300
        help = "minimum length of the loop"

        "--quantity", "-n"
        nargs = '?'
        arg_type = Int
        default = 1
        help = "the number of loops"

        "--quiet", "-q"
        action = :store_true
        help = "silence the output"
    end

    settings = parse_args(args, rules)

    !settings["quiet"] && println(crayon"red", "reading the file...")
    source, sr = wavread(settings["audio-path"])
    source = source[:, 1]

    !settings["quiet"] && println(crayon"yellow", "computing mfcc...")
    cepstrum = Float32.(first(mfcc(source, sr; wintime=1/s, steptime=1/s, numcep=numcep)))

    # since coefficients are strided, restride gives ~2x speedup
    cepstrum = restride(cepstrum)

    # number of samples in a single pack
    numsamples = Int(sr) ÷ s

    # length of the track (measured in the number of packs)
    numpacks = length(source) ÷ numsamples

    !settings["quiet"] && println(crayon"green", "extracting seaming points...")
    seams = extract(cepstrum, settings["min-seconds"] * s, settings["max-seconds"] * s, numpacks)

    matrix = Array{Float32}(undef, 2, length(seams))

    for idx = 1:length(seams)
        matrix[1, idx] = seams[idx].startidx
        matrix[2, idx] = seams[idx].endidx
    end

    !settings["quiet"] && println(crayon"blue", "trimming results...")
    indices = kmeans(matrix, settings["quantity"]) |> assignments

    sourcename = split(settings["audio-path"], '/') |> last |> s -> split(s, '.') |> first

    for idx = 1:settings["quantity"]
        seam = minimum(seams[findall(isequal(idx), indices)])

        loop = source[seam.startidx * numsamples : seam.endidx * numsamples]

        loopname = if settings["quantity"] == 1
            "$sourcename-loop.wav"
        else
            "$sourcename-loop-$idx.wav"
        end

        wavwrite(loop, loopname, Fs=sr, nbits=16, compression=WAVE_FORMAT_PCM)
        !settings["quiet"] && println(crayon"magenta", "congratulations, now run\nmpv $loopname --loop")
    end
end

main(ARGS)

max_stars_repo_name: ogoremeni/liederkreis

max_issues_repo_name: ogoremeni/liederkreis

max_forks_repo_name: ogoremeni/liederkreis


text:
#=
Author: Shuangchi He / Yulv
Email: yulvchi@qq.com
Date: 2022-02-15 12:32:43
LastEditors: Shuangchi He
LastEditTime: 2022-02-25 21:45:58
Description: Statistical Analysis for Pearson Correlation and Bland-Altman Agreement
=#
using Polynomials
using PyPlot
using Statistics

function Pearson_Correlation_Bland_Altman_Agreement(X, Y)
    # Pearson Correlation
    linearfit=fit(X, Y, 1)  # Linear fit.
    k = linearfit[1]
    b = linearfit[0]

    scatter(X, Y, label="Num: $(length(X))", color="r", s=8)
    plot(sort(X), sort(X * k .+ b), color="b", linestyle="--",
        label="linear_fit: y = $(round(k, digits=5)) x + $(round(b, digits=5))\npearson coefficient: $(round(cor(X, Y), digits=5))")
    plot([findmin([findmin(X)[1], findmin(Y)[1]])[1], findmax([findmax(X)[1], findmax(Y)[1]])[1]], [findmin([findmin(X)[1], findmin(Y)[1]])[1], findmax([findmax(X)[1], findmax(Y)[1]])[1]],
        color="k", linestyle="--", label="Equal Line")

    xlabel("Measurement_predict (mm)")
    ylabel("Measurement_GT (mm)")
    title("Pearson Correlation")
    legend(loc=0)
    savefig("Measurement_Pearson_Correlation.png")
    close()

    # Bland-Altman Agreement
    data_mean = (X + Y) / 2.
    data_diff = X - Y
    mean_diff = mean(data_diff)
    std_diff = std(data_diff)  # Julia std defaults to correcting for bias in sample variance by dividing by N-1.

    linearfit2=fit(data_mean, data_diff, 1)  # Linear fit.
    k2 = linearfit2[1]
    b2 = linearfit2[0]

    scatter(data_mean, data_diff, label="Diff num: $(length(data_mean))", color="r", s=8)
    plot(sort(data_mean), sort(mean_diff * ones(length(data_mean))), color="k", linestyle="--",
        label="Mean_Diff: $(round(mean_diff, digits=5))")
    plot(sort(data_mean), sort((mean_diff .+ 1.96 * std_diff) * ones(length(data_mean))), color="b", linestyle="--",
        label="Mean_Diff + 1.96 Std_Diff: $(round(mean_diff + 1.96 * std_diff, digits=5))")
    plot(sort(data_mean), sort((mean_diff .- 1.96 * std_diff) * ones(length(data_mean))), color="b", linestyle="-.",
        label="Mean_Diff - 1.96 Std_Diff: $(round(mean_diff - 1.96 * std_diff, digits=5))")
    plot(sort(data_mean), sort(data_mean * k2 .+ b2), color="y", linestyle="--",
        label="linear_fit: y = $(round(k2, digits=5)) x + $(round(b2, digits=5))")

    xlabel("Mean of Measurement_predict and Measurement_GT (mm)")
    ylabel("Diff between Measurement_predict and Measurement_GT (mm)")
    title("Bland-Altman Agreement")
    legend(loc=0)
    savefig("Measurement_Bland-Altman_Agreement.png")
end

X = [0.125, 0.95, 0.55, 0.60, 0.78, 0.46, 0.88, 0.50, 0.93, 0.35, 0.975, 0.725, 0.285, 0.166, 0.666, 0.888, 0.233]
Y = [0.127, 0.97, 0.53, 0.57, 0.72, 0.49, 0.91, 0.52, 0.90, 0.37, 0.982, 0.718, 0.277, 0.175, 0.666, 0.88, 0.2333]
Pearson_Correlation_Bland_Altman_Agreement(X, Y)