text:
#MDDatasets: Base broadcast support
#-------------------------------------------------------------------------------

#==Nomenclature:
A 1-argument function is one where 1 argument in particular dictates the
dimensonality of the operation.

A 2-argument function is one where 2 arguments dictate the dimensionality of
the operation.
==#

#==Type definitions
===============================================================================#
abstract type CastType end
abstract type CastTypeRed <: CastType end #Identifies a reducing function

#Identifies a function cast with 1 argument of TCAST:
struct CastType1{TCAST, POS} <: CastType; end
#Identifies a reducing function cast with 1 argument of TCAST, returning collection(Number):
struct CastTypeRed1{TCAST, POS} <: CastTypeRed; end
#Identifies a function cast with 2 argument of TCAST1/2:
struct CastType2{TCAST1, POS1, TCAST2, POS2} <: CastType; end
#Identifies a reducing function cast with 2 argument of TCAST1/2, returning collection(Number):
struct CastTypeRed2{TCAST1, POS1, TCAST2, POS2} <: CastTypeRed; end #Reducing function (DataF1->Number)

#Constructors:
CastType(::Type{T}, pos::Int) where T = CastType1{T, pos}()
CastType(::Type{T1}, pos1::Int, ::Type{T2}, pos2::Int) where {T1,T2} =
	CastType2{T1, pos1, T2, pos2}()

CastTypeRed(::Type{T}, pos::Int) where T = CastTypeRed1{T, pos}()
CastTypeRed(::Type{T1}, pos1::Int, ::Type{T2}, pos2::Int) where {T1,T2} =
	CastTypeRed2{T1, pos1, T2, pos2}()

#==Constants
===============================================================================#
#Cast on function capable of operating directly on base types (Number):
const CAST_BASEOP1 = CastType(Number, 1)
const CAST_BASEOPRED1 = CastTypeRed(Number, 1)

#Cast on function capable of operating directly on base types (Number, Number):
const CAST_BASEOP2 = CastType(Number, 1, Number, 2)

#Cast on function capable of operating only on DataF1:
const CAST_MD1 = CastType(DataF1, 1)
const CAST_MDRED1 = CastTypeRed(DataF1, 1)

#Cast on function capable of operating only on DataF1:
const CAST_MD2 = CastType(DataF1, 1, DataF1, 2)
const CAST_MDRED2 = CastTypeRed(DataF1, 1, DataF1, 2)

#==result_type: Figure out result type for a function call with given arguments
===============================================================================#
#=NOTES:
 - Advantage: Helps with readability.
 - Disadvantage: Might add extra level on call stack for simple operations (hope inlines).

 - result_type never returns concrete DataF1{TX,TY}.  Only abstract "DataF1".

 - (Future) special functions might behave differently.
   ie: The default return types specified here are not necessarily correct.
   Methods should be overwritten for special functions.
   Sadly: can't currently dispatch on function.  Need to think about this.
=#
#-------------------------------------------------------------------------------

result_type(fn::Function, ::Type{T1}) where T1 = T1 #In=out... wrong with functions like sin(Int)=>Float64
result_type(fn::Function, ::Type{T1}) where T1<:DataF1 = DataF1

result_type(fn::Function, ::Type{T1}, ::Type{T2}) where {T1<:Number, T2<:Number} = promote_type(T1, T2)
result_type(fn::Function, ::Type{T1}, ::Type{T2}) where {T1<:DataF1, T2<:DataF1} = DataF1
result_type(fn::Function, ::Type{T1}, ::Type{T2}) where {T1<:DataF1, T2<:Number} = DataF1
result_type(fn::Function, ::Type{T1}, ::Type{T2}) where {T1<:Number, T2<:DataF1} = DataF1

result_type(CT::CastTypeRed, fn::Function, ::Type{DataF1{TX,TY}}) where {TX,TY} = TY

#==eltype/valtype: Figure out result type for a function call with given arguments
===============================================================================#

#y-element type of an operation with a DataF1 (result: DataF1):
Base.valtype(fn::Function, d1::Number, d2::DataF1) =
	result_type(fn, typeof(d1), eltype(d2.y))
Base.valtype(fn::Function, d1::DataF1, d2::Number) =
	result_type(fn, eltype(d1.y), typeof(d2))
Base.valtype(fn::Function, d1::DataF1, d2::DataF1) =
	result_type(fn, eltype(d1.y), eltype(d2.y))

#=="apply" function
===============================================================================#

#"apply" between a DataF1 & basic scalar (Number)
function apply(fn::Function, d1::DataF1, d2::Number, args...; kwargs...)
	RT = valtype(fn, d1, d2)
	y = RT[fn(yi, d2, args...; kwargs...) for yi in d1.y]
	return DataF1(d1.x, y)
end
function apply(fn::Function, d1::Number, d2::DataF1, args...; kwargs...)
	RT = valtype(fn, d1, d2)
	y = RT[fn(d1, yi, args...; kwargs...) for yi in d2.y]
	return DataF1(d2.x, y)
end

#Relay to Base. broadcasting system:
#-------------------------------------------------------------------------------
broadcastBase(fn::Function, args...; kwargs...) =
	Base.materialize(Base.broadcasted(fn, args...; kwargs...))

#Generic broadcastMD functions (avoids name collisions):
#NOTE: "core" means most elemental version of the function
#-------------------------------------------------------------------------------
#Trap to provide more useful message:
function broadcastMD(ct::CastType, fn::Function, args...; kwargs...)
	msg = "Cast type not supported for call to $fn: $ct"
	throw(ArgumentError(msg))
end

#fn(DataF1) - core: fn(Number)
broadcastMD(::CastType1{Number,1}, fn::Function, d::Number, args...; kwargs...) =
	fn(d, args...; kwargs...)
broadcastMD(::CastType1{Number,1}, fn::Function, d::DataF1, args...; kwargs...) =
	DataF1(d.x, broadcastBase(fn, d.y, args...; kwargs...))

#Reducing fn(DataF1) - core: fn(Number)
broadcastMD(::CastTypeRed1{Number,1}, fn::Function, d::DataF1, args...; kwargs...) =
	fn(d.y, args...; kwargs...)

#fn(DataF1) - core: fn(DataF1)
broadcastMD(::CastType1{DataF1,1}, fn::Function, d::DataF1, args...; kwargs...) =
	fn(d, args...; kwargs...)

#Reducing fn(DataF1) - core: fn(DataF1)
broadcastMD(::CastTypeRed1{DataF1,1}, fn::Function, d::DataF1, args...; kwargs...) =
	fn(d, args...; kwargs...)

#fn(Number, Number) - core: fn(Number, Number)
broadcastMD(CT::CastType2{Number,1,Number,2}, fn::Function, d1::Number, d2::Number, args...; kwargs...) =
	fn(d1, d2, args...; kwargs...)

#fn(DataF1, DataF1) - core: fn(Number, Number)
broadcastMD(CT::CastType2{Number,1,Number,2}, fn::Function, d1::DataF1, d2::DataF1, args...; kwargs...) =
	apply(fn, d1, d2, args...; kwargs...)

#fn(DataF1, Number) - core: fn(Number, Number)
broadcastMD(CT::CastType2{Number,1,Number,2}, fn::Function, d1::DataF1, d2::Number, args...; kwargs...) =
	apply(fn, d1, d2, args...; kwargs...)

#fn(Number, DataF1) - core: fn(Number, Number)
broadcastMD(CT::CastType2{Number,1,Number,2}, fn::Function, d1::Number, d2::DataF1, args...; kwargs...) =
	apply(fn, d1, d2, args...; kwargs...)

#fn(DataF1, DataF1) - core: fn(DataF1, DataF1)
broadcastMD(CT::CastType2{DataF1,1,DataF1,2}, fn::Function, d1::DataF1, d2::DataF1, args...; kwargs...) =
	fn(d1, d2, args...; kwargs...)

#Last line


text:
# Base16 Encoder
# ==============

struct Base16Encoder <: Codec
    table::CodeTable16
    state::State
end

Base16Encoder(table::CodeTable16) = Base16Encoder(table, State())

"""
    Base16Encoder(;lowercase=false)

Create a base16 encoding codec.

Arguments
- `lowercase`: use [0-9a-f] instead of [0-9A-F].
"""
function Base16Encoder(;lowercase::Bool=false)
    if lowercase
        table = BASE16_LOWER
    else
        table = BASE16_UPPER
    end
    return Base16Encoder(table)
end

const Base16EncoderStream{S} = TranscodingStream{Base16Encoder,S} where S<:IO

"""
    Base16EncoderStream(stream::IO; kwargs...)

Create a base16 encoding stream (see `Base16Encoder` for `kwargs`).
"""
function Base16EncoderStream(stream::IO; kwargs...)
    return TranscodingStream(Base16Encoder(;kwargs...), stream)
end

macro encode16(i, j1, j2)
    quote
        a = input[$(i)]
        output[$(j1)] = encode(table, a >> 4)
        output[$(j2)] = encode(table, a & 0x0f)
    end |> esc
end

function TranscodingStreams.startproc(
        codec :: Base16Encoder,
        state :: Symbol,
        error :: Error)
    start!(codec.state)
    return :ok
end

function TranscodingStreams.process(
        codec  :: Base16Encoder,
        input  :: Memory,
        output :: Memory,
        error  :: Error)
    table = codec.table
    state = codec.state

    # Check if we can encode data.
    if !is_running(state)
        error[] = ArgumentError("encoding is already finished")
        return 0, 0, :error
    elseif input.size == 0
        finish!(state)
        return 0, 0, :end
    elseif output.size < 2
        # Need more output space.
        return 0, 0, :ok
    end

    # Encode the data.
    i = j = 0
    k::Int = min(fld(input.size, 4), fld(output.size, 8))
    @inbounds while k > 0  # ≡ i + 4 ≤ input.size && j + 8 ≤ output.size
        # unrolled loop
        @encode16 i+1 j+1 j+2
        @encode16 i+2 j+3 j+4
        @encode16 i+3 j+5 j+6
        @encode16 i+4 j+7 j+8
        i += 4
        j += 8
        k -= 1
    end
    @inbounds while i + 1 ≤ input.size && j + 2 ≤ output.size
        @encode16 i+1 j+1 j+2
        i += 1
        j += 2
    end
    return i, j, :ok
end

ext: jl


text:
# generate a tridiagonal matrix
function tridiagm(a, b, c, n)
       dd, du = ones(n), ones(n - 1)
       b*diagm(0 => dd) + a*diagm(-1 => du) + c*diagm(1 => du)
end

SUITE["Balls"] = BenchmarkGroup()
SUITE["Linear map"] = BenchmarkGroup()
begin
    for set_type in (Ball1, Ball2, BallInf)
        for n in (1, 10, 100, 1000)
            X = set_type(ones(n), 0.1)
            d = [i/n for i in 1:n]
            SUITE["Balls"][string(set_type), n] = @benchmarkable σ($d, $X)

            A = tridiagm(1, -2, 0.5, n)
            Y = A * X
            SUITE["Linear map"][string(set_type), n] = @benchmarkable σ($d, $Y)
        end
    end
end


text:
function solve_first(input)
	ret = 0
	lines = collect(Iterators.takewhile(line -> line != "", input))
	while length(lines) > 0
		ret += lines |> join |> unique |> length
		lines = collect(Iterators.takewhile(line -> line != "", input))
	end
	ret
end

function solve_second(input)
	ret = 0
	lines = collect(Iterators.takewhile(line -> line != "", input))
	while length(lines) > 0
		ret +=  count(char -> all(line -> contains(line, char), lines), collect(lines[1]))
		lines = collect(Iterators.takewhile(line -> line != "", input))
	end
	ret
end

@assert solve_first(eachline("input/test/test_day6.txt")) == 11 ["test case 1 failed"]
@assert solve_second(eachline("input/test/test_day6.txt")) == 6 ["test case 2 failed"]

println("First star: " * string(solve_first(eachline("input/day6.txt"))))
println("Second star: " * string(solve_second(eachline("input/day6.txt"))))


text:
function accumulative_backpropagation!(net, train_datum::TestDatum,
    dCx_dw_buffer, dCx_db_buffer, dCx_da_buffer, Bj_buffer)

    #feedforward
    feedforward!(net, train_datum);

    #BackPropagate
    update_Bj_buffer!(net, Bj_buffer, dCx_da_buffer, train_datum);
    accumulate_dCx_dw_buffer!(net, dCx_dw_buffer, Bj_buffer);
    accumulate_dCx_db_buffer!(net, dCx_db_buffer, Bj_buffer);

end

function update_Bj_buffer!(net, Bj_buffer, dCx_da_buffer, train_datum)

    #Reset Buffers
    for l in 2:net.depth
        Bj_buffer[l] .= 0.0;
        dCx_da_buffer[l] .= 0.0;
    end

    for l in net.depth:-1:2
        if l == net.depth
            for j in 1:outlength(net)
                dCx_da_buffer[net.depth][j] = single_cost_prime(output(net)[j], output(train_datum)[j]);
                Bj_buffer[net.depth][j] = act_function_prime(net.z[net.depth][j]) * dCx_da_buffer[net.depth][j];
            end
        else
            nextl = l + 1;
            for k in 1:net.lsizes[l]
                for j in 1:net.lsizes[nextl]
                    dCx_da_buffer[l][k] += net.w[nextl][j,k] * Bj_buffer[nextl][j];
                end
                Bj_buffer[l][k] = act_function_prime(net.z[l][k]) * dCx_da_buffer[l][k];
            end
        end
    end
end

function accumulate_dCx_db_buffer!(net, dCx_db_buffer, Bj_buffer)
    for l in 2:net.depth
        for j in 1:net.lsizes[l]
            dCx_db_buffer[l][j] += Bj_buffer[l][j];
        end
    end
end

function accumulate_dCx_dw_buffer!(net, dCx_dw_buffer, Bj_buffer)
    for l in 2:net.depth
        prevl = l-1;
        for j in 1:net.lsizes[l]
            for k in 1:net.lsizes[prevl]
                dCx_dw_buffer[l][j,k] += net.a[prevl][k] * Bj_buffer[l][j];
            end
        end
    end
end


text:
module FEAssembly

using SparseArrays
using LinearAlgebra
import Base.Filesystem

using ..FEStructure
using ..FESparse
using ..LoadCase

export assemble!,assemble_cable_link!,
clear_result!

mutable struct LCNode
    case
    children::Vector{LCNode}
end

mutable struct Assembly
    structure::Structure
    lcset::LoadCaseSet

    node_count::Int
    beam_count::Int
    quad_count::Int
    tria_count::Int

    nDOF::Int
    nfreeDOF::Int
    restrainedDOFs::Vector{Int}

    lc_tree::LCNode
    working_path::String
end

function DFSfind(lcnode::LCNode,lc::String)
    if lcnode.case !="root" && lcnode.case.id==lc
        return lcnode
    elseif isempty(lcnode.children)
        return nothing
    else
        for child in lcnode.children
            res=DFSfind(child,lc)
            if res isa LCNode
                return res
            end
        end
    end
end

function disperse(A::Matrix,i::Vector{Int},N::Int)
    m,n=size(A)
    I = repeat(i,outer=n)
    J = repeat(i,inner=m)
    V = reshape(A,m*n)
    SparseMatrixCOO(N,N,I,J,V)
end

function disperse(A::Vector,i::Vector{Int},N::Int)
    m=length(A)
    I = i
    J = [1 for i in 1:m]
    V = A
    to_array(SparseMatrixCOO(N,1,I,J,V))[:,1]
end

function disperse_idx(node::Node)
    i=node.hid
    collect(6*i-5:6*i)
end

function disperse_idx(elm::Beam)
    i = elm.node1.hid
    j = elm.node2.hid
    [6i-5:6i;6j-5:6j]
end

function disperse_idx(elm::Quad)
    i = elm.node1.hid
    j = elm.node2.hid
    k = elm.node3.hid
    l = elm.node4.hid
    [6i-5:6i;6j-5:6j;6k-5:6k;6l-5:6l]
end

function disperse_idx(elm::Tria)
    i = elm.node1.hid
    j = elm.node2.hid
    k = elm.node3.hid
    [6i-5:6i;6j-5:6j;6k-5:6k]
end

assembleK(elm,N)=disperse(elm.T'*integrateK(elm)*elm.T,disperse_idx(elm),N)
assembleM(elm,N)=disperse(elm.T'*integrateM(elm)*elm.T,disperse_idx(elm),N)

"""
    assemble!(structure,lcset;path=pwd())
集成结构与工况
# 参数
- `structure::Structure`: Structure类型实例
- `lcset::LoadCaseSet`: LoadCaseSet类型实例
# 关键字参数
- `mass_source::String`:`weight` or `loadcases`
- `mass_cases::Array{String}`: 质量荷载工况
- `mass_case_factors::Array{Number}`:工况叠加系数
- `path`: 工作路径
# 返回
- `assembly`: Assembly对象
"""
function assemble!(structure,lcset;mass_source="weight",mass_cases=[],mass_cases_factors=[],path=pwd())
    nDOF=length(structure.nodes)*6
    restrainedDOFs=[]
    node_count=length(structure.nodes)
    beam_count=length(structure.beams)
    quad_count=length(structure.quads)
    tria_count=length(structure.trias)

    K=spzeros_coo(nDOF,nDOF)
    K̄=spzeros(nDOF,nDOF)
    M=spzeros_coo(nDOF,nDOF)
    M̄=spzeros(nDOF,nDOF)
    C=sparse(0.02*I,nDOF,nDOF)
    C̄=sparse(0.02*I,nDOF,nDOF)

    lc_tree=LCNode("root",[])
    lcs=merge(lcset.statics,lcset.modals,lcset.bucklings,lcset.time_histories,lcset.response_spectrums)

    while !isempty(lcs)
        for k in keys(lcs)
            if lcs[k].plc==""
                lcnode=LCNode(pop!(lcs,k),[])
                push!(lc_tree.children,lcnode)#root
            else
                lcnode=DFSfind(lc_tree,lcs[k].plc)
                if lcnode isa LCNode
                    child=LCNode(pop!(lcs,k),[])
                    push!(lcnode.children,child)
                else
                    continue
                end
            end
        end
    end

    for node in values(structure.nodes)
        idx=node.hid*6-6
        for r in node.restraints
            idx+=1
            if r
                push!(restrainedDOFs,idx)
            end
        end
    end

    sort!(restrainedDOFs)
    mask=[(i in restrainedDOFs) ? false : true for i in 1:nDOF]

    # for elm in values(structure.beams)
    #     i = elm.node1.hid
    #     j = elm.node2.hid
    #     T=sparse(elm.T)
    #
    #     I=collect(1:12)
    #     J=[6i-5:6i;6j-5:6j]
    #     G=sparse(I,J,1.0,12,nDOF)
    #
    #     Kᵉ=integrateK!(elm)
    #     Mᵉ=integrateM!(elm)
    #
    #     A=T*G
    #
    #     rDOF=findall(x->x==true,elm.release)
    #     if length(rDOF)!=0
    #         K̄ᵉ,P̄ᵉ=FEStructure.FEBeam.static_condensation(Array(elm.Kᵉ),elm.Pᵉ,rDOF)
    #         K̄ᵉ=sparse(K̄ᵉ)
    #         K+=A'*K̄ᵉ*A
    #     else
    #         K+=A'*Kᵉ*A
    #     end
    #     M+=A'*Mᵉ*A
    # end

@time begin
    if node_count>0
        K+=reduce(+,assembleK.(values(structure.nodes),nDOF))
        M+=reduce(+,assembleM.(values(structure.nodes),nDOF))
    end
end

@time begin
    if beam_count>0
        K+=reduce(+,assembleK.(values(structure.beams),nDOF))
        M+=reduce(+,assembleM.(values(structure.beams),nDOF))
    end
end

@time begin
    if quad_count>0
        K+=reduce(+,assembleK.(values(structure.quads),nDOF))
        M+=reduce(+,assembleM.(values(structure.quads),nDOF))
    end
end

@time begin
    if tria_count>0
        K+=reduce(+,assembleK.(values(structure.trias),nDOF))
        M+=reduce(+,assembleM.(values(structure.trias),nDOF))
    end
end

    if structure.damp=="constant"
        C=sparse(structure.ζ₁*I,nDOF,nDOF)
    elseif structure.damp=="rayleigh"
        C=structure.ζ₁*K+structure.ζ₂*M
    else
        throw("Structural damping error!")
    end

    structure.K=to_csc(K)
    structure.C=C

    #calculate gravity
    _gset=LoadCaseSet()
    add_static_case!(_gset,"__Gravity__",1.0)
    Pg=zeros(nDOF)
    for elm in values(structure.beams)
        #load on repeat id will be added together
        f=[0,0,-9.81*elm.material.ρ*elm.section.A,0,0,0,
           0,0,-9.81*elm.material.ρ*elm.section.A,0,0,0]
        T=sparse(elm.T)
        fᵉ=T*f
        fi1,fi2,fi3,mi1,mi2,mi3,fi1,fi2,fi3,mi1,mi2,mi3=fᵉ
        add_beam_distributed!(_gset,"__Gravity__",elm.id,fi1,fi2,fi3,mi1,mi2,mi3,fi1,fi2,fi3,mi1,mi2,mi3)
    end
    #add quad / tria gravity
    for load in values(_gset.statics["__Gravity__"].beam_forces)
        elm=structure.beams[load.id]
        Kᵉ=integrateK(elm)
        Pᵉ=integrateP(elm,load)
        rDOF=findall(x->x==true,elm.release)
        K̃ᵉ,P̃ᵉ=FEStructure.static_condensation(Kᵉ,Pᵉ,rDOF)
        idx=disperse_idx(elm)
        Pg+=disperse(reshape(elm.T'*P̃ᵉ,12),idx,nDOF)
    end

    cases=merge(lcset.statics,lcset.bucklings,lcset.time_histories)
    for l in keys(cases)
        loadcase=cases[l]
        P=zeros(nDOF)
        if l in keys(lcset.statics)
            P.+=Pg*loadcase.gfactor
        end
        for load in values(loadcase.beam_forces)
            elm=structure.beams[load.id]
            Kᵉ=integrateK(elm)
            Pᵉ=integrateP(elm,load)
            rDOF=findall(x->x==true,elm.release)
            K̃ᵉ,P̃ᵉ=FEStructure.static_condensation(Kᵉ,Pᵉ,rDOF)
            idx=disperse_idx(elm)
            P+=disperse(reshape(elm.T'*P̃ᵉ,12),idx,nDOF)
        end
        for load in values(loadcase.nodal_forces)
            node=structure.nodes[load.id]
            Pⁿ=reshape(node.T'*load.val,6)
            idx=disperse_idx(node)
            P+=disperse(Pⁿ,idx,nDOF)
        end
        cases[l].P=P
    end

    if mass_source=="weight"
        structure.M=to_csc(M)
    elseif mass_source=="loadcases"
        M=zero(M)
        for (lc,fac) in (mass_cases, mass_case_factors)
            M+=Diagonal(lcset.statics[lc].P)*fac
        end
        structure.M=sparse(M)
    else
        throw("mass_source should be only weight or loadcases")
    end
    return Assembly(structure,lcset,node_count,beam_count,quad_count,tria_count,
                    nDOF,nDOF-length(restrainedDOFs),restrainedDOFs,lc_tree,path)
end

function clear_result!(assembly)
    nDOF=assembly.nDOF
    assembly.structure.K=spzeros(nDOF,nDOF)
    assembly.structure.K̄=spzeros(nDOF,nDOF)
    assembly.structure.M=spzeros(nDOF,nDOF)
    assembly.structure.M̄=spzeros(nDOF,nDOF)
    assembly.structure.C=spzeros(nDOF,nDOF)
    assembly.structure.C̄=spzeros(nDOF,nDOF)
    for i in keys(assembly.lcset.statics)
        assembly.lcset.statics[i].P=zeros(nDOF)
        assembly.lcset.statics[i].P̄=zeros(nDOF)
    end
    assembly.restrainedDOFs=[]
    result_dir=joinpath(assembly.working_path,".analysis")
    rm(result_dir,recursive=true)
    return
end

end


text:
#SpiceData: A pure Julia SPICE data reader
#-------------------------------------------------------------------------------
__precompile__(true)

#=
TAGS:
	#WANTCONST, HIDEWARN_0.7
=#

module SpiceData

include("base.jl")
include("show.jl")

#==DataReader object: public members
================================================================================
	.sweepname
	.signalnames
	.sweep #Sweep values
==#

#==Exported symbols
===============================================================================#

#==Un-"exported" symbols
================================================================================
	_open(filepath::String)::DataReader
==#

#==Other interface tools (symbols not exported to avoid collisions):
================================================================================
#Already in base:
	Base.names(reader::DataReader)
	Base.read(reader::DataReader, signum::Int)
	Base.read(reader::DataReader, signame::String)
==#

end #module


text:
#
#--------------------------------------#
export reluHat
#--------------------------------------#
function reluHat(xg::AbstractVector  # grid points
                ,ug::AbstractVector  # grid values
                ,x ::AbstractArray)  # eval points, size [1,N]

    NN = reluHatFunc(xg)
    return reluHat(NN,ug,x)
end
#--------------------------------------#
function reluHat(NN::Function        # reluHatFunc
                ,ug::AbstractVector  # grid values
                ,x ::AbstractArray)  # eval points

    return (ug' * NN(x)) # size [1,N]
end
#--------------------------------------#
export reluHatFunc
#--------------------------------------#
"""
 Generate FEM Hat functions
"""
function reluHatFunc(x0::Number,x1::Number,x2::Number)
    return reluHatFunc([x0,x1,x2])
end
#--------------------------------------#
function reluHatFunc(xg::AbstractVector)

    n  = length(xg)
    h  = diff(xg)
    hi = 1.0 ./ h

    RR,b2 = Zygote.ignore() do
        RR = restrictUpper(n) + restrictLower(n)
        RR = RR'
        b2 = zero(xg); b2[1] = 1
        return RR,b2
    end

    hh = -RR*hi
    W2 = Tdg(hi,hh,hi)

    NN = Chain(x -> relu.(x .- xg) # fix FEM nodes
              ,u -> (W2*u .+ b2)   # form FEM hat functions (linear)
              )

    return (x -> NN(x))

#   return femPOU(n,x->x,xg[2:end-1],xg[1],xg[end])
end
#--------------------------------------#
export addEndPts
#--------------------------------------#
function addEndPts(xe::AbstractArray
                  ,x0::Number,x1::Number)

    n = length(xe) + 2
    A,R = Zygote.ignore() do
        A = zeros(n,2)
        A[1,1] = 1
        A[n,2] = 1
        R = restrict(n)
        return A,R
    end

    xg = A * [x0,x1] + R'*xe

    return xg
end
#--------------------------------------#


text:
using Documenter, MixedSubdivisions

makedocs(
    sitename = "MixedSubdivisions.jl",
    pages = [
        "MixedSubdivisions" => "index.md",
    ]
)

deploydocs(
    repo   = "github.com/saschatimme/MixedSubdivisions.jl.git",
)


text:
# test module for prob_stats

module TestProbStats
    using Test

    # target modules
    include(joinpath(split(@__FILE__, "test")[1], "src/prob_stats/sensor_data/print_sensor_data.jl"))
    include(joinpath(split(@__FILE__, "test")[1], "src/prob_stats/freq_dist/histogram/draw_histogram.jl"))
    include(joinpath(split(@__FILE__, "test")[1], "src/prob_stats/freq_dist/histogram/histogram_mean.jl"))
    include(joinpath(split(@__FILE__, "test")[1], "src/prob_stats/freq_dist/variance/calc_variance.jl"))
    include(joinpath(split(@__FILE__, "test")[1], "src/prob_stats/freq_dist/std_dev/calc_std_dev.jl"))
    include(joinpath(split(@__FILE__, "test")[1], "src/prob_stats/prob_dist/prob_dist.jl"))
    include(joinpath(split(@__FILE__, "test")[1], "src/prob_stats/prob_model/gauss_dist_model/gauss_dist_model.jl"))
    include(joinpath(split(@__FILE__, "test")[1], "src/prob_stats/prob_model/gauss_prob_dist/gauss_prob_dist.jl"))
    include(joinpath(split(@__FILE__, "test")[1], "src/prob_stats/prob_model/gauss_prob_dist_pdf/gauss_prob_dist_pdf.jl"))
    include(joinpath(split(@__FILE__, "test")[1], "src/prob_stats/prob_model/gauss_prob_dist_cdf/gauss_prob_dist_cdf.jl"))
    include(joinpath(split(@__FILE__, "test")[1], "src/prob_stats/prob_model/bar_prob_dist_cdf/bar_prob_dist_cdf.jl"))
    include(joinpath(split(@__FILE__, "test")[1], "src/prob_stats/prob_model/expected_value/calc_expected_value.jl"))
    include(joinpath(split(@__FILE__, "test")[1], "src/prob_stats/complex_dist/histogram_600mm/draw_hist_600.jl"))
    include(joinpath(split(@__FILE__, "test")[1], "src/prob_stats/complex_dist/time_series_600mm/draw_time_series_600.jl"))
    include(joinpath(split(@__FILE__, "test")[1], "src/prob_stats/complex_dist/group_by_hour/group_by_hour.jl"))
    include(joinpath(split(@__FILE__, "test")[1], "src/prob_stats/complex_dist/hist_grp_by_hour/hist_grp_by_hour.jl"))
    include(joinpath(split(@__FILE__, "test")[1], "src/prob_stats/complex_dist/group_hour_heatmap/group_hour_heatmap.jl"))
    include(joinpath(split(@__FILE__, "test")[1], "src/prob_stats/complex_dist/prob_sum_hour/prob_sum_hour.jl"))
    include(joinpath(split(@__FILE__, "test")[1], "src/prob_stats/complex_dist/prob_sum_lidar/prob_sum_lidar.jl"))
    include(joinpath(split(@__FILE__, "test")[1], "src/prob_stats/complex_dist/cond_z_t_bar/cond_z_t_bar.jl"))
    include(joinpath(split(@__FILE__, "test")[1], "src/prob_stats/complex_dist/probs_bayes/probs_bayes.jl"))
    include(joinpath(split(@__FILE__, "test")[1], "src/prob_stats/complex_dist/bayes_estimation_630/bayes_estimation_630.jl"))
    include(joinpath(split(@__FILE__, "test")[1], "src/prob_stats/complex_dist/bayes_est_at_5/bayes_est_at_5.jl"))
    include(joinpath(split(@__FILE__, "test")[1], "src/prob_stats/complex_dist/bayes_est_at_11/bayes_est_at_11.jl"))
    include(joinpath(split(@__FILE__, "test")[1], "src/prob_stats/multi_dim_gauss_dist/covariance/calc_covariance.jl"))
    include(joinpath(split(@__FILE__, "test")[1], "src/prob_stats/multi_dim_gauss_dist/contour_pdf/contour_pdf.jl"))
    include(joinpath(split(@__FILE__, "test")[1], "src/prob_stats/multi_dim_gauss_dist/contour_pdf/contour_pdf_plus20.jl"))
    include(joinpath(split(@__FILE__, "test")[1], "src/prob_stats/multi_dim_gauss_dist/contour_pdf/contour_pdf_200.jl"))
    include(joinpath(split(@__FILE__, "test")[1], "src/prob_stats/multi_dim_gauss_dist/multiple_gauss_dist/multiple_gauss_dist.jl"))
    include(joinpath(split(@__FILE__, "test")[1], "src/prob_stats/multi_dim_gauss_dist/eigen/calc_plot_eigen.jl"))
    include(joinpath(split(@__FILE__, "test")[1], "src/prob_stats/multi_dim_gauss_dist/cov_matrix/calc_cov_mat.jl"))

    function main()
        @testset "ProbStats" begin
            @testset "SensorData" begin
                @test_nowarn PrintSensorData.print_data_200()
                @test_nowarn PrintSensorData.print_data_600()
            end
            @testset "FreqDist" begin
                @test_nowarn DrawHistogram.main(true)
                @test_nowarn HistogramMean.main(true)
                @test_nowarn CalcVariance.main()
                @test_nowarn CalcStdDev.main()
            end
            @testset "ProbDist" begin
                @test_nowarn ProbDist.main(true)
            end
            @testset "ProbModel" begin
                @testset "GaussDistModel" begin
                    @test GaussDistModel.p(190) == 2.064932530633184e-5
                    @test GaussDistModel.p(210) == 0.08231272044708592
                    @test GaussDistModel.p(230) == 1.2364903952298874e-5
                    @test_nowarn GaussDistModel.main(true)
                end
                @testset "GaussProbDist" begin
                    @test GaussProbDist.prob(193) == 0.0002254354061377026
                    @test GaussProbDist.prob(210) == 0.08187587024811953
                    @test GaussProbDist.prob(229) == 3.113715670488921e-5
                    @test_nowarn GaussProbDist.main(true)
                end
                @testset "GaussProbDistPdf" begin
                    @test_nowarn GaussProbDistPdf.main(true)
                end
                @testset "GaussProbDistCdf" begin
                    @test_nowarn GaussProbDistCdf.main(true)
                end
                @testset "BarProbDistCdf" begin
                    @test_nowarn BarProbDistCdf.main(true)
                end
                @testset "CalcExpValue" begin
                    @test_nowarn CalcExpValue.main()
                end
            end
            @testset "ComplexDist" begin
                @testset "DrawHist600" begin
                    @test_nowarn DrawHist600.main(true)
                end
                @testset "DrawTimeSeries600" begin
                    @test_nowarn DrawTimeSeries600.main(true)
                end
                @testset "GroupByHour" begin
                    @test_nowarn GroupByHour.main(true)
                end
                @testset "HistGrpByHour" begin
                    @test_nowarn HistGrpByHour.main(true)
                end
                @testset "GrpHrHeatmap" begin
                    @test_nowarn GrpHrHeatmap.main(true)
                end
                @testset "ProbSumHour" begin
                    @test ProbSumHour.main(true) == 1.0
                end
                @testset "ProbSumLidar" begin
                    @test ProbSumLidar.main(true) == 1.0
                end
                @testset "CondZtBar" begin
                    @test_nowarn CondZtBar.main(true)
                end
                @testset "ProbsBayes" begin
                    @test_nowarn ProbsBayes.main()
                end
                @testset "BayesEst630" begin
                    @test_nowarn BayesEst630.main(true)
                end
                @testset "BayesEstAt5" begin
                    @test_nowarn BayesEstAt5.main(true)
                end
                @testset "BayesEstAt11" begin
                    @test_nowarn BayesEstAt11.main(true)
                end
            end
            @testset "MultiDimGuassDist" begin
                @testset "CalcCovariance" begin
                    @test_nowarn CalcCovariance.main()
                end
                @testset "MultiGaussDist" begin
                    @test_nowarn MultiGaussDist.main(true)
                end
                @testset "CalcPlotEigen" begin
                    @test_nowarn CalcPlotEigen.main(true)
                end
                @testset "CalcCovMat" begin
                    @test_nowarn CalcCovMat.main()
                end
            end
        end
    end
end

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy


text:
using PyPlot
# plt[:style][:use]("ggplot")
plt[:rcParams]["font.family"][1] = "serif"
# include("../src/DSR.jl")
using DSR

function test()

    T, zta = 1.0, 0.02
    xi, alpha = 0.008, 0.2
    kpa = 1.0
    phi = 0.005
    omg = 2.0*pi/T
    dt = 0.005
    n = 4001

    # t = linspace(0,(n-1)*dt,n)
    # p = (t) -> sin(2.0*π/T*t) # 正弦激励

    s0 = DSR.SDF(omg, zta)
    # s  = DSR.SDF_MYD_BW(omg, zta, kpa, phi, 0.02)
    # s.n = 10
    # s  = DSR.SDF_MYD_EP(omg, zta, kpa, phi, 0.05)
    # s  = DSR.SDF_FD(omg, zta, phi)
    # s  = DSR.SDF_VD(omg, zta, xi, alpha)

    mu, xi, kpa = 0.050, 0.02, 0.085
    s  = DSR.SDF_PID_VD(omg, zta, mu, xi, kpa)
    # s  = DSR.SDF_MAXWELL_VD(omg, zta, xi, kpa, alpha)
    # s  = DSR.SDF_MAXWELL_VD(omg, zta, xi*1000.0^(1.0-alpha), kpa, alpha)

    # r = DSR.response(s, p, t)

    a_g = DSR.readtxt("EQ-S-1.txt")
    dt_g = 0.005
    ndiv = 1
    dt = dt_g/ndiv
    nstep = convert(Int64,round(length(a_g)*ndiv*1.1))

    r0 = DSR.response(s0, a_g, dt_g, dt, nstep)
    r  = DSR.response(s,  a_g, dt_g, dt, nstep)

    figure(1,(9,7.5))
    subplot(3,1,1)
    plot(r0[:,1],r0[:,2])
    plot(r[:,1], r[:,2])
    grid(true)

    # subplot(3,1,2)
    # plot(t, r[:,2])
    # plot(t,r0[:,2])
    # grid(true)
    # subplot(3,1,3)
    # plot(t, r[:,3])
    # plot(t,r0[:,3])
    # grid(true)

    subplot2grid((3,1),(1,0),rowspan=2,colspan=1)
    plot(r[:,2],r[:,5])
    # plot(r[:,2],r[:,6])
    plot(r[:,7],r[:,8])
    grid(true)

    # show()
end

test();

max_stars_repo_name: Panchatantra/DSR

max_issues_repo_name: Panchatantra/DSR

max_forks_repo_name: Panchatantra/DSR


text:
# Define a sector for ungraded vector spaces
struct Trivial <: Sector end
Base.iterate(::SectorValues{Trivial}, i = false) = return i ? nothing : (Trivial(), true)
Base.IteratorSize(::Type{SectorValues{Trivial}}) = HasLength()
Base.length(::SectorValues{Trivial}) = 1
Base.getindex(::SectorValues{Trivial}, i::Int) = i == 1 ? Trivial() :
    throw(BoundsError(values(Trivial), i))
findindex(::SectorValues{Trivial}, c::Trivial) = 1
Base.show(io::IO, ::Trivial) = print(io, "Trivial()")

Base.one(::Type{Trivial}) = Trivial()
Base.conj(::Trivial) = Trivial()
Base.isreal(::Type{Trivial}) = true
Base.isless(::Trivial, ::Trivial) = false

FusionStyle(::Type{Trivial}) = UniqueFusion()
⊗(::Trivial, ::Trivial) = (Trivial(),)
Nsymbol(::Trivial, ::Trivial, ::Trivial) = true
Fsymbol(::Trivial, ::Trivial, ::Trivial, ::Trivial, ::Trivial, ::Trivial) = 1

BraidingStyle(::Type{Trivial}) = Bosonic()
Rsymbol(::Trivial, ::Trivial, ::Trivial) = 1


text:
__precompile__()

module DifferentialEquations

  using Reexport

  @reexport using DiffEqBase
  @reexport using DiffEqNoiseProcess
  @reexport using RecursiveArrayTools

  @reexport using SteadyStateDiffEq
  @reexport using StochasticDiffEq
  @reexport using OrdinaryDiffEq
  @reexport using BoundaryValueDiffEq
  using Sundials
  @reexport using DelayDiffEq

  @reexport using DiffEqCallbacks
  @reexport using DiffEqMonteCarlo
  @reexport using DiffEqJump

  @reexport using DiffEqFinancial
  @reexport using MultiScaleArrays
  @reexport using DiffEqPhysics

  @reexport using DimensionalPlotRecipes

  using LinearAlgebra

  import DiffEqBase: solve

  include("default_solve.jl")
  include("default_arg_parsing.jl")
  include("ode_default_alg.jl")
  include("sde_default_alg.jl")
  include("dae_default_alg.jl")
  include("dde_default_alg.jl")
  include("discrete_default_alg.jl")
  include("rode_default_alg.jl")
  include("steady_state_default_alg.jl")

end # module


text:
s = "Hello,How,Are,You,Today"
a = split(s, ",")
t = join(a, ".")

println("The string \"", s, "\"")
println("Splits into ", a)
println("Reconstitutes to \"", t, "\"")

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
"""
	grdcontour(cmd0::String="", arg1=[]; kwargs...)

Reads a 2-D grid file or a GMTgrid type and produces a contour map by tracing each
contour through the grid.

Full option list at [`pscontour`](http://gmt.soest.hawaii.edu/doc/latest/pscontour.html)

Parameters
----------

- $(GMT.opt_J)
- **A** : **annot** : -- Str or Number --       Flags = [-|[+]annot_int][labelinfo]

    *annot_int* is annotation interval in data units; it is ignored if contour levels are given in a file.
    [`-A`](http://gmt.soest.hawaii.edu/doc/latest/grdcontour.html#a)
- $(GMT.opt_B)
- **C** : **cont** : **contours** : **levels** : -- Str or Number --

    Contours to be drawn.
    [`-C`](http://gmt.soest.hawaii.edu/doc/latest/grdcontour.html#c)
- **D** : **dump** : -- Str --

    Dump contours as data line segments; no plotting takes place.
    [`-D`](http://gmt.soest.hawaii.edu/doc/latest/grdcontour.html#d)
- **F** : **force** : -- Str or [] --

    Force dumped contours to be oriented so that higher z-values are to the left (-Fl [Default]) or right.
    [`-F`](http://gmt.soest.hawaii.edu/doc/latest/grdcontour.html#f)
- **G** : **labels** : -- Str --

    Controls the placement of labels along the quoted lines.
    [`-G`](http://gmt.soest.hawaii.edu/doc/latest/grdcontour.html#g)
- $(GMT.opt_Jz)
- **L** : **range** : -- Str --

    Limit range: Do not draw contours for data values below low or above high.
    [`-L`](http://gmt.soest.hawaii.edu/doc/latest/grdcontour.html#l)
- **N** : **fill** : -- Bool or [] --

    Fill the area between contours using the discrete color table given by cpt.
    [`-N`](http://gmt.soest.hawaii.edu/doc/latest/grdcontour.html#n)
- $(GMT.opt_P)
- **Q** : **cut** : -- Str or Number --

    Do not draw contours with less than cut number of points.
    [`-Q`](http://gmt.soest.hawaii.edu/doc/latest/grdcontour.html#q)
- **S** : **smooth** : -- Number --

    Used to resample the contour lines at roughly every (gridbox_size/smoothfactor) interval.
    [`-S`](http://gmt.soest.hawaii.edu/doc/latest/grdcontour.html#s)
- **T** : **ticks** : -- Str --

    Draw tick marks pointing in the downward direction every *gap* along the innermost closed contours.
    [`-T`](http://gmt.soest.hawaii.edu/doc/latest/grdcontour.html#t)
- $(GMT.opt_R)
- $(GMT.opt_U)
- $(GMT.opt_V)
- **W** : **pen** : -- Str or Number --

    Sets the attributes for the particular line.
    [`-W`](http://gmt.soest.hawaii.edu/doc/latest/grdcontour.html#w)
- $(GMT.opt_X)
- $(GMT.opt_Y)
- **Z** : **scale** : -- Str --

    Use to subtract shift from the data and multiply the results by factor before contouring starts.
    [`-Z`](http://gmt.soest.hawaii.edu/doc/latest/grdcontour.html#z)
- $(GMT.opt_bo)
- $(GMT.opt_do)
- $(GMT.opt_e)
- $(GMT.opt_f)
- $(GMT.opt_h)
- $(GMT.opt_p)
- $(GMT.opt_t)
"""
function grdcontour(cmd0::String="", arg1=[], arg2=[]; K=false, O=false, first=true, kwargs...)

	length(kwargs) == 0 && return monolitic("grdcontour", cmd0, arg1)	# Speedy mode

	d = KW(kwargs)
	output, opt_T, fname_ext = fname_out(d)		# OUTPUT may have been an extension only

	cmd  = ""
    cmd, opt_B, = parse_BJR(d, cmd, "", O, " -JX12c/0")
	cmd  = parse_UVXY(cmd, d)
	cmd, = parse_bo(cmd, d)
	cmd, = parse_e(cmd, d)
	cmd, = parse_f(cmd, d)
	cmd, = parse_h(cmd, d)
	cmd, = parse_p(cmd, d)
	cmd, = parse_t(cmd, d)
	cmd  = parse_params(cmd, d)

	cmd, K, O, opt_B = set_KO(cmd, opt_B, first, K, O)		# Set the K O dance
	cmd, got_fname, arg1 = find_data(d, cmd0, cmd, 1, arg1)	# Find how data was transmitted
	if (isa(arg1, Array{<:Number}))		arg1 = mat2grid(arg1)	end

	N_used = got_fname == 0 ? 1 : 0		# To know whether a cpt will go to arg1 or arg2
	cmd, arg1, arg2, = add_opt_cpt(d, cmd, [:C :color :cmap], 'C', N_used, arg1, arg2)

	cmd = add_opt(cmd, 'A', d, [:A :annot])
	if (!occursin(" -C", cmd))			# Otherwise ignore an eventual :cont because we already have it
		cmd = add_opt(cmd, 'C', d, [:cont :contours :levels])
	end

	for symb in [:N :fill]
		if (haskey(d, symb))
			if (isa(d[symb], GMTcpt))
				if (!isempty_(arg2))	# Already have one cpt in arg2, replace it by new one
					arg2 = []
				end
				cmd, arg1, arg2, = add_opt_cpt(d, cmd, [:N :fill], 'N', N_used, arg1, arg2)
			else
				cmd = cmd * " -N"
			end
			break
		end
	end

	cmd = add_opt(cmd, 'D', d, [:D :dump])
	cmd = add_opt(cmd, 'F', d, [:F :force])
	cmd = add_opt(cmd, 'G', d, [:G :labels])
	cmd = add_opt(cmd, 'L', d, [:L :range])
	cmd = add_opt(cmd, 'Q', d, [:Q :cut])
	cmd = add_opt(cmd, 'S', d, [:S :smooth])
	cmd = add_opt(cmd, 'T', d, [:T :ticks])
	cmd = add_opt(cmd, 'W', d, [:W :pen])
	cmd = add_opt(cmd, 'Z', d, [:Z :scale])

	cmd = finish_PS(d, cmd, output, K, O)
    return finish_PS_module(d, cmd, "-D", output, fname_ext, opt_T, K, "grdcontour", arg1, arg2)
end

# ---------------------------------------------------------------------------------------------------
grdcontour!(cmd0::String="", arg1=[], arg2=[]; K=true, O=true, first=false, kw...) =
	grdcontour(cmd0, arg1, arg2; K=true, O=true, first=false, kw...)

grdcontour(arg1, cmd0::String="", arg2=[]; K=false, O=false, first=true, kw...) =
	grdcontour(cmd0, arg1, arg2; K=K, O=O, first=first, kw...)

grdcontour!(arg1, cmd0::String="", arg2=[]; K=true, O=true, first=false, kw...) =
	grdcontour(cmd0, arg1, arg2; K=true, O=true, first=false, kw...)


text:
@testsection "constructors" begin
for T in (Float32, Float64, F64), dim in (1,2,3), order in (1,2,4)
    for op in (rand, zero, ones, randn)
        # Tensor, SymmetricTensor
        for TensorType in (Tensor, SymmetricTensor)
            TensorType == SymmetricTensor && order == 1 && continue
            N = Tensors.n_components(TensorType{order, dim})
            t = (@inferred (op)(TensorType{order, dim}))::TensorType{order, dim, Float64}
            t = (@inferred (op)(TensorType{order, dim, T}))::TensorType{order, dim, T}
            t = (@inferred (op)(TensorType{order, dim, T, N}))::TensorType{order, dim, T}
            t = (@inferred (op)(t))::TensorType{order, dim, T}

            op == zero && @test zero(TensorType{order, dim, T}) == zeros(T, size(t))
            op == ones && @test ones(TensorType{order, dim, T}) == ones(T, size(t))
        end
        # Vec
        if order == 1
            (@inferred (op)(Vec{dim}))::Tensor{order, dim, Float64}
            (@inferred (op)(Vec{dim, T}))::Tensor{order, dim, T}
        end
    end
    # Special Vec constructor
    if order == 1
        t = ntuple(i -> T(i), dim)
        @test (@inferred Vec(t))::Tensor{1,dim,T,dim} == Vec{dim}(t)
        @test (@inferred Vec(t...))::Tensor{1,dim,T,dim} == Vec{dim}(t)
    end
    for TensorType in (Tensor, SymmetricTensor), (func, el) in ((zeros, zero), (ones, one))
        TensorType == SymmetricTensor && order == 1 && continue
        order == 1 && func == ones && continue # one not supported for Vec's
        N = Tensors.n_components(TensorType{order, dim})
        tens_arr1 = func(TensorType{order, dim}, 1)
        tens_arr2 = func(TensorType{order, dim, T}, 2, 2)
        tens_arr3 = func(TensorType{order, dim, T, N}, 3, 3, 3)
        @test tens_arr1[1] == tens_arr2[1, 1] == tens_arr3[1, 1, 1] == el(TensorType{order, dim, T})
        @test eltype(tens_arr1) == TensorType{order, dim, Float64, N}
        @test eltype(tens_arr2) == eltype(tens_arr3) == TensorType{order, dim, T, N}
    end
end
for dim in (1, 2, 3)
    # Heterogeneous tuple/type unstable function
    z(i, jkl...) = i % 2 == 0 ? 0 : float(0)
    @test Vec{dim}(ntuple(z, dim))::Vec{dim,Float64} ==
          Vec{dim}(z)::Vec{dim,Float64}
    # @test Vec{dim,Float32}(ntuple(z, dim))::Vec{dim,Float32} ==
    #       Vec{dim,Float32}(z)::Vec{dim,Float32}
    for order in (1, 2, 4)
        N = Tensors.n_components(Tensor{order,dim})
        @test Tensor{order,dim}(ntuple(z, N))::Tensor{order,dim,Float64} ==
              Tensor{order,dim}(z)::Tensor{order,dim,Float64}
        @test Tensor{order,dim,Float32}(ntuple(z, N))::Tensor{order,dim,Float32} ==
              Tensor{order,dim,Float32}(z)::Tensor{order,dim,Float32}
        order == 1 && continue
        N = Tensors.n_components(SymmetricTensor{order,dim})
        @test SymmetricTensor{order,dim}(ntuple(z, N))::SymmetricTensor{order,dim,Float64} ==
              SymmetricTensor{order,dim}(z)::SymmetricTensor{order,dim,Float64}
        @test SymmetricTensor{order,dim,Float32}(ntuple(z, N))::SymmetricTensor{order,dim,Float32} ==
              SymmetricTensor{order,dim,Float32}(z)::SymmetricTensor{order,dim,Float32}
    end
end
# Number type which is not <: Real but <: Number (Tensors#154)
@test Vec{3, NotReal}((1, 2, 3)) isa Vec{3, NotReal}
end # of testset

@testsection "diagm, one" begin
for T in (Float32, Float64), dim in (1,2,3)
    # diagm
    v = rand(T, dim)
    vt = (v...,)

    @test (@inferred diagm(Tensor{2, dim}, v))::Tensor{2, dim, T} == diagm(0 => v)
    @test (@inferred diagm(Tensor{2, dim}, vt))::Tensor{2, dim, T} == diagm(0 => v)
    @test (@inferred diagm(SymmetricTensor{2, dim}, v))::SymmetricTensor{2, dim, T} == diagm(0 => v)
    @test (@inferred diagm(SymmetricTensor{2, dim}, vt))::SymmetricTensor{2, dim, T} == diagm(0 => v)

    v = rand(T); vv = v * ones(T, dim)
    @test (@inferred diagm(Tensor{2, dim}, v))::Tensor{2, dim, T} == diagm(0 => vv)
    @test (@inferred diagm(SymmetricTensor{2, dim}, v))::SymmetricTensor{2, dim, T} == diagm(0 => vv)

    # one
    @test one(Tensor{2, dim, T}) == diagm(Tensor{2, dim}, one(T)) == Matrix(I, dim, dim)
    @test one(SymmetricTensor{2, dim, T}) == diagm(SymmetricTensor{2, dim}, one(T)) == Matrix(I, dim, dim)

    M = 1 # dummy
    @test one(Tensor{2, dim, T, M}) == one(Tensor{2, dim, T})
    @test one(SymmetricTensor{2, dim, T, M}) == one(SymmetricTensor{2, dim, T})

    _I =  (@inferred one(Tensor{2, dim, T}))::Tensor{2, dim, T}
    II = (@inferred one(Tensor{4, dim, T}))::Tensor{4, dim, T}
    I_sym =  (@inferred one(SymmetricTensor{2, dim, T}))::SymmetricTensor{2, dim, T}
    II_sym = (@inferred one(SymmetricTensor{4, dim, T}))::SymmetricTensor{4, dim, T}
    for i in 1:dim, j in 1:dim
        if i == j
            @test _I[i,j] == T(1)
            @test I_sym[i,j] == T(1)
        else
            @test _I[i,j] == T(0)
            @test I_sym[i,j] == T(0)
        end
        for k in 1:dim, l in 1:dim
            if i == k && j == l
                @test II[i,j,k,l] == T(1)
                if i == l && j == k
                    @test II_sym[i,j,k,l] == T(1)
                else
                    @test II_sym[i,j,k,l] == T(1) / 2
                end
            else
                @test II[i,j,k,l] == T(0)
                if i == l && j == k
                    @test II_sym[i,j,k,l] == T(1) / 2
                else
                    @test II_sym[i,j,k,l] == T(0)
                end
            end
        end
    end
end
end # of testset

@testsection "base vectors" begin
for T in (Float32, Float64, F64), dim in (1,2,3)
    eᵢ_func(i) = Tensor{1, dim, T}(j->j==i ? one(T) : zero(T))

    a = rand(Vec{dim, T})
    for i in 1:dim
        @test eᵢ(a, i) == eᵢ(typeof(a), i) == eᵢ(a)[i] == eᵢ(typeof(a))[i] == eᵢ_func(i)
    end

    b = zero(a)
    for i in 1:dim
        @test a[i] == eᵢ(a, i) ⋅ a
        b += eᵢ(a, i) * a[i]
    end
    @test a ≈ b
end
end # of testset

@testsection "simple math" begin
for T in (Float32, Float64), dim in (1,2,3), order in (1,2,4), TensorType in (Tensor, SymmetricTensor)
    TensorType == SymmetricTensor && order == 1 && continue
    t = rand(TensorType{order, dim, T})

    # Binary tensor tensor: +, -
    @test (@inferred t + t)::TensorType{order, dim} == Array(t) + Array(t)
    @test (@inferred 2*t)::TensorType{order, dim} == 2 * Array(t)
    @test (@inferred t - t)::TensorType{order, dim} == Array(t) - Array(t)
    @test (@inferred 0*t)::TensorType{order, dim} == 0 * Array(t)

    # Binary tensor number: *, /
    @test (@inferred 0.5 * t)::TensorType{order, dim} ≈ 0.5 * Array(t)
    @test (@inferred t * 0.5)::TensorType{order, dim} ≈ Array(t) * 0.5
    @test (@inferred t / 2.0)::TensorType{order, dim} ≈ Array(t) / 2.0

    # Unary: +, -
    @test (@inferred +t)::TensorType{order, dim} == zero(t) + t
    @test (@inferred -t)::TensorType{order, dim} == zero(t) - t

    if order == 2
        # Power by literal integer
        fm3, fm2, fm1, f0, fp1, fp2, fp3 = t -> t^-3, t -> t^-2, t -> t^-1, t -> t^0, t -> t^1, t -> t^2, t -> t^3
        @test (@inferred fm3(t))::typeof(t) ≈ inv(t) ⋅ inv(t) ⋅ inv(t)
        @test (@inferred fm2(t))::typeof(t) ≈ inv(t) ⋅ inv(t)
        @test (@inferred fm1(t))::typeof(t) ≈ inv(t)
        @test (@inferred f0(t))::typeof(t)  ≈ one(t)
        @test (@inferred fp1(t))::typeof(t) ≈ t
        @test (@inferred fp2(t))::typeof(t) ≈ t ⋅ t
        @test (@inferred fp3(t))::typeof(t) ≈ t ⋅ t ⋅ t
    end
end
end # of testset

@testsection "construct from function" begin
for T in (Float32, Float64, F64)
    for dim in (1,2,3)
        fi = (i) -> cos(i)
        fij = (i,j) -> cos(i) + sin(j)
        fijkl = (i, j, k ,l) -> cos(i) + sin(j) + tan(k) + exp(l)

        vf = (@inferred Vec{dim, T}(fi))::Tensor{1, dim, T}
        af = (@inferred Tensor{1, dim, T}(fi))::Tensor{1, dim, T}
        Af = (@inferred Tensor{2, dim, T}(fij))::Tensor{2, dim, T}
        AAf = (@inferred Tensor{4, dim, T}(fijkl))::Tensor{4, dim, T}
        Af_sym = (@inferred SymmetricTensor{2, dim, T}(fij))::SymmetricTensor{2, dim, T}
        AAf_sym = (@inferred SymmetricTensor{4, dim, T}(fijkl))::SymmetricTensor{4, dim, T}

        for i in 1:dim
            @test vf[i] == T(fi(i))
            @test af[i] == T(fi(i))
            for j in 1:dim
                @test Af[i,j] == T(fij(i, j))
                for k in 1:dim, l in 1:dim
                    @test AAf[i,j,k,l] == T(fijkl(i,j,k,l))
                end
            end
        end

        for j in 1:dim, i in j:dim
            @test Af_sym[i,j] == T(fij(i, j))
            for l in 1:dim, k in l:dim
                 @test AAf_sym[i,j,k,l] == T(fijkl(i,j,k,l))
            end
        end

        # Test fully specified constructors
        for vec in (vf,af)
            @test vec == typeof(vec)(i->vec[i])
        end
        for t2 in (Af, Af_sym)
            @test t2 == typeof(t2)((i,j)->t2[i,j])
        end
        for t4 in (AAf, AAf_sym)
            @test t4 == typeof(t4)((i,j,k,l) -> t4[i,j,k,l])
        end
    end
end
end # of testset

@testsection "construct from Array" begin
for (T1, T2) in ((Float32, Float64), (Float64, Float32)), order in (1,2,4), dim in (1,2,3)
    At = rand(Tensor{order, dim})
    gen_data = rand(T1, size(At))

    @test (@inferred Tensor{order, dim}(gen_data))::Tensor{order, dim, T1}     ≈ gen_data
    @test (@inferred Tensor{order, dim, T2}(gen_data))::Tensor{order, dim, T2} ≈ gen_data

    if order != 1
        As = rand(SymmetricTensor{order, dim})
        gen_data = rand(T1, size(As))
        At = Tensor{order, dim, T1}(gen_data)
        Ats = symmetric(At)
        gen_sym_data_full = Array(Ats)
        gen_sym_data = T1[Ats.data[i] for i in 1:length(Ats.data)]

        @test (@inferred SymmetricTensor{order, dim}(gen_sym_data_full))::SymmetricTensor{order, dim, T1}     ≈ Ats
        @test (@inferred SymmetricTensor{order, dim}(gen_sym_data))::SymmetricTensor{order, dim, T1}          ≈ Ats
        @test (@inferred SymmetricTensor{order, dim, T2}(gen_sym_data_full))::SymmetricTensor{order, dim, T2} ≈ Ats
        @test (@inferred SymmetricTensor{order, dim, T2}(gen_sym_data))::SymmetricTensor{order, dim, T2}      ≈ Ats
    end
end
end # of testset

@testsection "indexing" begin
for T in (Float32, Float64, F64), dim in (1,2,3), order in (1,2,4)
    if order == 1
        data = rand(T, dim)
        vect = Tensor{order, dim, T}(data)
        for i in 1:dim+1
            if i > dim
                @test_throws BoundsError vect[i]
            else
                @test vect[i] ≈ data[i]
            end
        end
        @test (@inferred vect[:])::Vec{dim, T} == vect
    elseif order == 2
        data = rand(T, dim, dim)
        symdata = data + data'
        S = Tensor{order,dim, T}(data)
        Ssym = SymmetricTensor{order,dim, T}(symdata)
        @test_throws ArgumentError S[:]
        @test_throws ArgumentError Ssym[:]
        for i in 1:dim+1, j in 1:dim+1
            if i > dim || j > dim
                @test_throws BoundsError S[i, j]
                @test_throws BoundsError Ssym[i, j]
            else
                @test S[i, j] ≈ data[i, j]
                @test Ssym[i, j] ≈ symdata[i, j]
                # Slice
                @test (@inferred S[i,:])::Tensor{1, dim, T} ≈ data[i,:]
                @test (@inferred S[:,j])::Tensor{1, dim, T} ≈ data[:,j]
                @test (@inferred Ssym[i,:])::Tensor{1, dim, T} ≈ symdata[i,:]
                @test (@inferred Ssym[:,j])::Tensor{1, dim, T} ≈ symdata[:,j]
            end
        end
    elseif order == 4
        data = rand(T,dim,dim,dim,dim)
        S = Tensor{order,dim, T}(data)
        Ssym = symmetric(S)
        symdata = Array(Ssym)
        @test_throws ArgumentError S[:]
        @test_throws ArgumentError Ssym[:]
        for i in 1:dim+1, j in 1:dim+1, k in 1:dim+1, l in 1:dim+1
            if i > dim || j > dim || k > dim || l > dim
                @test_throws BoundsError S[i, j, k, l]
                @test_throws BoundsError Ssym[i, j, k, l]
            else
                @test S[i, j, k, l] ≈ data[i, j, k, l]
                @test Ssym[i, j, k, l] ≈ symdata[i, j, k, l]
            end
        end
    end
end
end # of testset

@testsection "norm, trace, det, inv, eig" begin
for T in (Float32, Float64, F64), dim in (1,2,3)
    # norm
    for order in (1,2,4)
        t = rand(Tensor{order, dim, T})
        @test (@inferred norm(t))::T ≈ sqrt(sum(abs2, Array(t)[:]))
        if order != 1
            t_sym = rand(SymmetricTensor{order, dim, T})
            @test (@inferred norm(t_sym))::T ≈ sqrt(sum(abs2, Array(t_sym)[:]))
        end
    end

    # trace, vol, dev, det, inv (only for second order tensors)
    t = rand(Tensor{2, dim, T})
    t_sym = rand(SymmetricTensor{2, dim, T})

    @test (@inferred tr(t))::T == sum([t[i,i] for i in 1:dim])
    @test (@inferred tr(t_sym))::T == sum([t_sym[i,i] for i in 1:dim])

    @test tr(t) ≈ mean(t)*3.0
    @test tr(t_sym) ≈ mean(t_sym)*3.0

    @test (@inferred vol(t))::Tensor{2, dim, T} ≈ mean(t) * Matrix(I, dim, dim)
    @test (@inferred vol(t_sym))::SymmetricTensor{2, dim, T} ≈ mean(t_sym) * Matrix(I, dim, dim)

    @test (@inferred dev(t))::Tensor{2, dim, T} ≈ Array(t) - 1/3*tr(t)* Matrix(I, dim, dim)
    @test (@inferred dev(t_sym))::SymmetricTensor{2, dim, T} ≈ Array(t_sym) - 1/3*tr(t_sym)* Matrix(I, dim, dim)

    @test (@inferred det(t))::T ≈ det(Array(t))
    @test (@inferred det(t_sym))::T ≈ det(Array(t_sym))

    @test (@inferred inv(t))::Tensor{2, dim, T} ≈ inv(Array(t))
    @test (@inferred inv(t_sym))::SymmetricTensor{2, dim, T} ≈ inv(Array(t_sym))

    # inv for fourth order tensors
    Random.seed!(1234)
    AA = rand(Tensor{4, dim, T})
    AA_sym = rand(SymmetricTensor{4, dim, T})
    @test AA ⊡ (@inferred inv(AA))::Tensor{4, dim, T} ≈ one(Tensor{4, dim, T})
    @test AA_sym ⊡ (@inferred inv(AA_sym))::SymmetricTensor{4, dim, T} ≈ one(SymmetricTensor{4, dim, T})

    E = @inferred eigen(t_sym)
    Λ, Φ = E
    Λa, Φa = eigen(Array(t_sym))

    @test Λ ≈ (@inferred eigvals(t_sym)) ≈ eigvals(E) ≈ Λa
    @test Φ ≈ (@inferred eigvecs(t_sym)) ≈ eigvecs(E)
    for i in 1:dim
        # scale with first element of eigenvector to account for possible directions
        @test Φ[:, i]*Φ[1, i] ≈ Φa[:, i]*Φa[1, i]
    end

    # test eigenfactorizations for a diagonal tensor
    v = rand(T, dim)
    d_sym = diagm(SymmetricTensor{2, dim, T}, v)
    E = @inferred eigen(d_sym)
    Λ, Φ = E
    Λa, Φa = eigen(Symmetric(Array(d_sym)))

    @test Λ ≈ (@inferred eigvals(d_sym)) ≈ eigvals(E) ≈ Λa
    @test Φ ≈ (@inferred eigvecs(d_sym)) ≈ eigvecs(E)

    # sqrt
    Apd = tdot(t_sym)
    @test sqrt(Apd) ⋅ sqrt(Apd) ≈ Apd
end
end # of testset

@testsection "eigen(::FourthOrderTensor)" begin
for T in (Float32, Float64), dim in (1, 2, 3)
    Random.seed!(123)
    # construct positive definite Voigt-tensor
    n = dim*dim - div((dim-1)*dim, 2)
    A = rand(T, n, n); A = A'A + I
    Aval, Avec = eigen(A)
    perm = sortperm(Aval)
    Aval = Aval[perm]
    Avec = [Avec[:, i] for i in perm]

    # construct tensor with "inverse scaling"
    S = fromvoigt(SymmetricTensor{4,dim,T}, A)
    S′ = SymmetricTensor{4,dim}((i,j,k,l) -> k == l ? S[i,j,k,l] : S[i,j,k,l]/2)

    E = eigen(S′)
    @test eigvals(E) ≈ Aval
    @test tovoigt.(eigvecs(E)) ≈ Avec
    for i in 1:n
        @test S′ ⊡ E.vectors[i] ≈ E.values[i] * E.vectors[i]
    end
    a, b = E # iteration
    @test a == eigvals(E) == E.values
    @test b == eigvecs(E) == E.vectors
end
end

# https://en.wikiversity.org/wiki/Continuum_mechanics/Tensor_algebra_identities
@testsection "tensor identities" begin
for T in (Float32, Float64, F64)
    for dim in (1,2,3)
        # Identities with second order and first order
        A = rand(Tensor{2, dim, T})
        B = rand(Tensor{2, dim, T})
        C = rand(Tensor{2, dim, T})
        I = one(Tensor{2, dim, T})
        a = rand(Tensor{1, dim, T})
        b = rand(Tensor{1, dim, T})

        @test A ⊡ B ≈ (A' ⋅ B) ⊡ one(A)
        @test A ⊡ (a ⊗ b) ≈ (A ⋅ b) ⋅ a
        @test (A ⋅ a) ⋅ (B ⋅ b) ≈ (A' ⋅ B) ⊡ (a ⊗ b)
        @test (A ⋅ a) ⊗ b ≈ A ⋅ (a ⊗ b)
        @test a ⊗ (A ⋅ b) ≈ (A ⋅ (b ⊗ a))'
        @test a ⊗ (A ⋅ b) ≈ (a ⊗ b) ⋅ A'

        @test A ⊡ I ≈ tr(A)
        @test det(A) ≈ det(A')
        @test tr(inv(A) ⋅ A) ≈ dim
        @test inv(A) ⋅ A ≈ A \ A ≈ I
        @test inv(A) ⋅ a ≈ A \ a

        @test (I ⊗ I) ⊡ A ≈ tr(A) * I
        @test (I ⊗ I) ⊡ A ⊡ A ≈ tr(A)^2

        @test A ⋅ a ≈ a ⋅ A'

        A_sym = rand(SymmetricTensor{2, dim})
        B_sym = rand(SymmetricTensor{2, dim})
        C_sym = rand(SymmetricTensor{2, dim})
        I_sym = one(SymmetricTensor{2, dim})

        @test A_sym ⊡ I_sym ≈ tr(A_sym)
        @test det(A_sym) ≈ det(A_sym')

        @test (I_sym ⊗ I_sym) ⊡ A_sym ≈ tr(A_sym) * I_sym
        @test ((I_sym ⊗ I_sym) ⊡ A_sym) ⊡ A_sym ≈ tr(A_sym)^2
    end
end

for T in (Float32, Float64, F64)
    for dim in (1,2,3)
        # Identities with identity tensor
        II = one(Tensor{4, dim, T})
        I = one(Tensor{2, dim, T})
        AA = rand(Tensor{4, dim, T})
        A = rand(Tensor{2, dim, T})
        II_sym = one(SymmetricTensor{4, dim, T})
        I_sym = one(SymmetricTensor{2, dim, T})
        AA_sym = rand(SymmetricTensor{4, dim, T})
        A_sym = rand(SymmetricTensor{2, dim, T})

        @test II ⊡ AA ≈ AA
        @test AA ⊡ II ≈ AA
        @test II ⊡ A ≈ A
        @test A ⊡ II ≈ A
        @test II ⊡ A ⊡ A ≈ (tr(A' ⋅ A))

        @test II ⊡ AA_sym ≈ AA_sym
        @test AA_sym ⊡ II ≈ AA_sym
        @test II ⊡ A_sym ≈ A_sym
        @test A_sym ⊡ II ≈ A_sym
        @test II ⊡ A_sym ⊡ A_sym ≈ (tr(A_sym' ⋅ A_sym))

        @test II_sym ⊡ AA_sym ≈ AA_sym
        @test AA_sym ⊡ II_sym ≈ AA_sym
        @test II_sym ⊡ A_sym ≈ A_sym
        @test A_sym ⊡ II_sym ≈ A_sym
        @test II_sym ⊡ A_sym ⊡ A_sym ≈ (tr(A_sym' ⋅ A_sym))
    end
end
end # of testset

@testsection "promotion/conversion" begin
T = Float32
WIDE_T = widen(T)
for dim in (1,2,3), order in (1,2,4)

    tens = Tensor{order, dim, T, dim^order}
    tens_wide = Tensor{order, dim, WIDE_T, dim^order}

    @test promote_type(tens, tens) == tens
    @test promote_type(tens_wide, tens) == tens_wide
    @test promote_type(tens, tens_wide) == tens_wide
    @test promote_type(tens_wide, tens_wide) == tens_wide

    At = rand(tens)
    Bt = rand(tens_wide)
    @test isa((@inferred At + Bt), tens_wide)
    @test isa((@inferred convert(Tensor, At)), tens)
    @test isa((@inferred convert(Tensor{order, dim, T}, At)), tens)
    @test isa((@inferred convert(tens, At)), tens)
    @test isa((@inferred convert(Tensor{order, dim, WIDE_T}, At)), tens_wide)
    @test isa((@inferred convert(tens_wide, At)), tens_wide)

    if order != 1
        M = Tensors.n_components(SymmetricTensor{order, dim})
        sym = SymmetricTensor{order, dim, T, M}
        sym_wide = SymmetricTensor{order, dim, WIDE_T, M}

        @test promote_type(sym, sym) == sym
        @test promote_type(sym_wide, sym) == sym_wide
        @test promote_type(sym, sym_wide) == sym_wide
        @test promote_type(sym_wide, sym_wide) == sym_wide

        @test promote_type(sym, tens) == tens
        @test promote_type(tens, sym) == tens
        @test promote_type(sym_wide, tens) == tens_wide
        @test promote_type(tens, sym_wide) == tens_wide
        @test promote_type(sym, tens_wide) == tens_wide
        @test promote_type(tens_wide, sym) == tens_wide
        @test promote_type(sym_wide, tens_wide) == tens_wide

        As = rand(sym)
        Bs = rand(sym_wide)
        @test isa((@inferred As + Bs), sym_wide)

        @test isa((@inferred convert(SymmetricTensor, As)), sym)
        @test isa((@inferred convert(SymmetricTensor{order, dim, T}, As)), sym)
        @test isa((@inferred convert(sym, As)), sym)
        @test isa((@inferred convert(SymmetricTensor{order, dim, WIDE_T}, As)), sym_wide)
        @test isa((@inferred convert(sym_wide, As)), sym_wide)

        # SymmetricTensor -> Tensor
        @test (@inferred convert(Tensor, As))::tens ≈ As ≈ Array(As)
        @test (@inferred convert(Tensor{order, dim}, As))::tens ≈ As ≈ Array(As)
        @test (@inferred convert(Tensor{order, dim, WIDE_T}, As))::tens_wide ≈ As ≈ Array(As)
        @test (@inferred convert(tens_wide, As))::tens_wide ≈ As ≈ Array(As)

        # Tensor -> SymmetricTensor
        if dim != 1
            @test_throws InexactError convert(SymmetricTensor, At)
            @test_throws InexactError convert(SymmetricTensor{order, dim}, At)
            @test_throws InexactError convert(SymmetricTensor{order, dim, WIDE_T}, At)
            @test_throws InexactError convert(typeof(Bs), At)
        end
    end
end
end  # of testset

@testset "issues #100, #101: orders of eigenvectors" begin
    for (i, j) in ((1, 2), (2, 1))
        S = SymmetricTensor{2,2,Float64}((i, 0, j))
        Λ = diagm(SymmetricTensor{2,2}, eigvals(S))
        Φ = eigvecs(S)
        @test Φ ⋅ Λ ⋅ Φ' ≈ S
    end
    for (i, j, k) in ((1, 2, 3), (1, 3, 2), (2, 1, 3), (2, 3, 1), (3, 1, 2), (3, 2, 1))
        S = SymmetricTensor{2,3,Float64}((i, 0, 0, j, 0, k))
        Λ = diagm(SymmetricTensor{2,3}, eigvals(S))
        Φ = eigvecs(S)
        @test Φ ⋅ Λ ⋅ Φ' ≈ S
    end
end

@testsection "exceptions" begin
    # normal multiplication
    A = rand(Tensor{2, 3})
    B = rand(Tensor{2, 3})
    @test_throws Exception A*B
    @test_throws Exception A'*B

    AA = rand(Tensor{4, 2})
    A2 = rand(Tensor{2, 2})
    @test_throws DimensionMismatch A + A2
    @test_throws DimensionMismatch AA - A

    # transpose/adjoint of Vec
    x = rand(Vec{3})
    @test_throws ArgumentError x'
    @test_throws ArgumentError transpose(x)
    @test_throws ArgumentError adjoint(x)
end # of testset


text:
# Juno.progress
function progress(f; name = "")
    _id = gensym()
    @info name progress=0.0 _id=_id
    try
        f(_id)
    finally
        @info name progress="done" _id=_id
    end
end

function demo1()
    progress(name="outer") do id
        for i in 1:10
            sleep(1e-3)
            @info "outer" progress=i/10
        end
    end
end

function demo2()
    progress(name="outer") do id
        for i in 1:10
            progress(name="inner") do id
                for j in 1:10
                    sleep(1e-3)
                    @info "inner" progress=i/10
                end
            end
            @info "outer" progress=i/10
        end
    end
end


text:
"""
	sphdistance(cmd0::String="", arg1=nothing, kwargs...)

Create Voronoi distance, node, or natural nearest-neighbor grid on a sphere

Full option list at [`sphdistance`]($(GMTdoc)sphdistance.html)

Parameters
----------

- **C** | **save_mem** :: [Type => Bool]

    For large data sets you can save some memory (at the expense of more processing).
    ($(GMTdoc)sphdistance.html#a)
- **E** | **what_quantity** :: [Type => Str]   ``Arg = d|n|z[dist]``

    Specify the quantity that should be assigned to the grid nodes.
    ($(GMTdoc)sphdistance.html#e)
- **G** | **grid** | **outgrid** :: [Type => Str]

    Output grid file name. Note that this is optional and to be used only when saving
    the result directly on disk. Otherwise, just use the G = sphdistance(....) form.
    ($(GMTdoc)sphdistance.html#g)
- **I** | **inc** :: [Type => Str or Number]

    *x_inc* [and optionally *y_inc*] is the grid spacing.
    ($(GMTdoc)sphdistance.html#i)
- **L** | **dist_unit** :: [Type => Str]      ``Arg = d|e|f|k|M|n|u``

    Specify the unit used for distance calculations.
    ($(GMTdoc)sphdistance.html#l)
- **N** | **nodes** :: [Type => Str]      ``Arg = nodes``

    Read the information pertaining to each Voronoi polygon (the unique node lon, lat and polygon area)
    from a separate file.
    ($(GMTdoc)sphdistance.html#n)
- **Q** | **voronoi** :: [Type => Str]     ``Arg = voronoifile``

    Append the name of a file with pre-calculated Voronoi polygons.
    ($(GMTdoc)sphdistance.html#q)
- $(GMT.opt_R)
- $(GMT.opt_V)
- $(GMT.opt_b)
- $(GMT.opt_d)
- $(GMT.opt_e)
- $(GMT.opt_h)
- $(GMT.opt_i)
- $(GMT.opt_r)
- $(GMT.opt_swap_xy)
"""
function sphdistance(cmd0::String="", arg1=nothing, arg2=nothing; kwargs...)

	length(kwargs) == 0 && occursin(" -", cmd0) && return monolitic("sphdistance ", cmd0, arg1, arg2)

	d = KW(kwargs)
	cmd = parse_common_opts(d, "", [:I :R :V_params :b :d :e :h :i :r :yx])
	cmd = parse_these_opts(cmd, d, [[:C :save_mem], [:E :what_quantity], [:G :grid :outgrid], [:L :dist_unit]])

	(arg1 === nothing) ? N_used = 0 : N_used = 1
	symbs = [[:Q :voronoi], [:N :nodes]];	flags = "QN"	# Process option -Q & -N
	for k = 1:2
		if ((val = find_in_dict(d, symbs[k])[1]) !== nothing)
			cmd *= " -" * flags[k]
			if (isa(val, GMTdataset) || isa(val, Array{GMTdataset}) || (isa(val, Array{<:Number}) && k == 2) )
				(N_used == 0) ? arg1 = val : arg2 = val
			else
				cmd *= arg2str(val)
			end
		end
	end

	common_grd(d, "sphdistance " * cmd, arg1, arg2)
end

# ---------------------------------------------------------------------------------------------------
sphdistance(arg1; kw...) = sphdistance("", arg1; kw...)


text:
@testset "Solvers" begin
  data1D = readgeotable(joinpath(datadir,"data1D.tsv"), delim='\t', coordnames=[:x])
  data2D = readgeotable(joinpath(datadir,"data2D.tsv"), delim='\t', coordnames=[:x,:y])
  grid1D = RegularGrid{Float64}(100)
  grid2D = RegularGrid{Float64}(100,100)

  @testset "Kriging" begin
    problem1D = EstimationProblem(data1D, grid1D, :value)
    problem2D = EstimationProblem(data2D, grid2D, :value)

    global_kriging = Kriging(
      :value => (variogram=GaussianVariogram(range=35.,nugget=0.),)
    )
    nearest_kriging = Kriging(
      :value => (variogram=GaussianVariogram(range=35.,nugget=0.), maxneighbors=3)
    )
    local_kriging = Kriging(
      :value => (variogram=GaussianVariogram(range=35.,nugget=0.),
                 maxneighbors=3, neighborhood=BallNeighborhood(100.))
    )

    solvers = [global_kriging, nearest_kriging, local_kriging]
    snames  = ["GlobalKriging", "NearestKriging", "LocalKriging"]

    solutions1D = [solve(problem1D, solver) for solver in solvers]
    solutions2D = [solve(problem2D, solver) for solver in solvers]

    # basic checks
    for solution in solutions2D
      M, V = solution[:value]
      @test M[26,26] == 1.
      @test M[51,76] == 0.
      @test M[76,51] == 1.
    end

    if visualtests
      gr(size=(800,400))
      for i in 1:2
        solution, sname = solutions1D[i], snames[i]
        @plottest plot(solution) joinpath(datadir,sname*"1D.png") !istravis
      end
      for i in 1:2
        solution, sname = solutions2D[i], snames[i]
        @plottest contourf(solution) joinpath(datadir,sname*"2D.png") !istravis
      end
      # TODO: test local_kriging
    end
  end

  @testset "SeqGaussSim" begin
    nreals = 3

    @testset "Conditional" begin
      problem = SimulationProblem(data2D, grid2D, :value, nreals)

      solver = SeqGaussSim(
        :value => (variogram=GaussianVariogram(range=35.),
                   neighborhood=BallNeighborhood(10.))
      )

      Random.seed!(2017)
      solution = solve(problem, solver)

      # basic checks
      reals = solution[:value]
      @test all(reals[i][26,26] == 1. for i in 1:nreals)
      @test all(reals[i][51,76] == 0. for i in 1:nreals)
      @test all(reals[i][76,51] == 1. for i in 1:nreals)

      if visualtests
        gr(size=(800,400))
        # @plottest plot(solution) joinpath(datadir,"SGSCond2D.png") !istravis
      end
    end

    @testset "Unconditional" begin
      problem = SimulationProblem(grid2D, :value => Float64, nreals)

      solver = SeqGaussSim(
        :value => (variogram=GaussianVariogram(range=35.),
                   neighborhood=BallNeighborhood(10.))
      )

      Random.seed!(2017)
      solution = solve(problem, solver)

      if visualtests
        gr(size=(800,400))
        # @plottest plot(solution) joinpath(datadir,"SGSUncond2D.png") !istravis
      end
    end
  end
end


text:
using .ProximalOperators

"""
    ls_sparse_spectral_lpv(Y,X,V,w,Nv::Int; λ = 1, coulomb = false, normalize=true)

Perform LPV spectral estimation using the method presented in
Bagge Carlson et al. "Linear Parameter-Varying Spectral Decomposition."
modified to include a sparsity-promoting L1 group-lasso penalty on the coefficients.
The groups are based on frequency, meaning a solution in which either all parameters
For a particular frequency are zero, or all are non-zero.
This is useful in the identification of frequency components among a large set of possible frequencies.
See the paper or README For additional details.

- `Y` output
- `X` sample locations
- `V` scheduling signal
- `w` frequency vector
- `Nv` number of basis functions
- `λ` Regularization parameter
- `coulomb` Assume discontinuity at `v=0` (useful For signals where, e.g., Coulomb friction might cause issues.)\n
- `normalize` Use normalized basis functions (See paper For details).

See `?ADMM` for keyword arguments to control the solver.

See also `psd`, `ls_spectral_lpv` and `ls_windowpsd_lpv`
"""
function ls_sparse_spectral_lpv(y::AbstractVector{S}, X::AbstractVector{S}, V::AbstractVector{S},
    w, Nv::Integer;
    λ         = 1,
    coulomb   = false,
    normalize = true,
    kwargs...) where S

    w        = w[:]
    T        = length(y)
    Nf       = length(w)
    K        = basis_activation_func(V,Nv,normalize,coulomb)
    M(w,X,V) = vec(vec(exp.(im.*w.*X))*K(V)')'
    As       = zeros(complex(S),T,ifelse(coulomb,2,1)*Nf*Nv)

    for n = 1:T
        As[n,:] = M(w,X[n],V[n])
    end

    x      = zeros(S, 2size(As,2))                 # Initialize with standard least squares
    inds   = reshape(1:length(x), Nf, :)'[:] # Permute parameters so that groups are adjacent
    inds   = vcat(inds...)
    # x    = [real.(params); imag.(params)][inds]
    Φ      = [real.(As) imag.(As)][:,inds]
    proxf  = ProximalOperators.LeastSquares(Φ,y,iterative=true)

    gs     = ntuple(f->NormL2(S(λ)), Nf)
    indsg  = ntuple(f->((f-1)*2Nv+1:f*2Nv, ) ,Nf)
    proxg  = SlicedSeparableSum(gs, indsg)
    local x, z
    try
        x,z      = ADMM(x, proxf, proxg; kwargs...)
    catch e
        if e isa InterruptException
            @info "Aborting"
            z = copy(x)
        else
            rethrow(e)
        end
    end
    z      = z[sortperm(inds)]            # Sortperm is inverse of inds
    params = complex.(z[1:end÷2], z[end÷2+1:end])
    SpectralExt(y, X, V, w, Nv, λ, coulomb, normalize, params, nothing)
end

"""`x,f = ls_sparse_spectral(y,t,f=default_freqs(t), [window::AbstractVector]; λ=1,
proxg      = ProximalOperators.NormL1(λ),
kwargs...)`

perform spectral estimation using the least-squares method with (default) a L1-norm penalty on the
Fourier coefficients, change kwarg `proxg` to e.g. `NormL0(λ)` for a different behavior or ` proxg = IndBallL0(4)` if the number of frequencies is known in advance. Promotes a sparse spectrum. See `?ADMM` for keyword arguments to control the solver.

`y` is the signal to be analyzed
`t` is the sampling points
`f` is a vector of frequencies
"""
function ls_sparse_spectral(y::AbstractArray{T},t,f=default_freqs(t);
    init       = false,
    λ          = T(1),
    proxg      = NormL1(λ),
    kwargs...) where T

    A,zerofreq  = get_fourier_regressor(T.(t),T.(f))
    params = init ? fourier_solve(A,y,zerofreq,λ) : fill(zero(T), length(f)) # Initialize with standard least squares
    if zerofreq === nothing
        x = [real.(params); imag.(params)]
    else
        x = [real.(params); imag.(params[2:end])]
    end
    proxf  = ProximalOperators.LeastSquares(A,y, iterative=true)
    x,z    = ADMM(x, proxf, proxg; kwargs...)
    params = fourier2complex(z, zerofreq)
    params, f
end

function ls_sparse_spectral(y::AbstractArray{T},t,f, W;
    init       = false,
    λ          = T(1),
    proxg      = NormL1(T(λ)),
    kwargs...) where T

    A,zerofreq  = get_fourier_regressor(t,f)
    params = init ? fourier_solve(A,y,zerofreq,λ) : fill(zero(T), length(f)) # Initialize with standard least squares
    if zerofreq === nothing
        x = [real.(params); imag.(params)]
    else
        x = [real.(params); imag.(params[2:end])]
    end
    Wd     = Diagonal(W)
    Q      = A'Wd*A
    q      = A'Wd*y
    proxf  = ProximalOperators.Quadratic(Q, q, iterative=true)
    # proxf  = ProximalOperators.LeastSquares(A,Wd*y, iterative=true)
    x,z = ADMM(x, proxf, proxg; kwargs...)
    params = fourier2complex(z, zerofreq)
    params, f
end

"""
    ADMM(x,proxf,proxg;
    iters      = 10000,   # ADMM maximum number of iterations
    tol        = 1e-5,    # ADMM tolerance
    printerval = 100,     # Print this often
    cb(x,z)    = nothing, # Callback function
    μ          = 0.05`)   # ADMM tuning parameter. If results oscillate, lower this value.
"""
function ADMM(x::AbstractArray{T},proxf,proxg;
    iters      = 10000,
    tol        = 1e-5,
    printerval = 100,
    cb         = nothing,
    μ          = T(0.05)) where T

    @assert 0 ≤ μ ≤ 1 "μ should be ≤ 1"
    μ = T(μ)

    z     = copy(x)
    u     = zeros(T,size(x))
    tmp   = similar(u)
    for i = 1:iters
        tmp .= z.-u
        prox!(x, proxf, tmp, μ)
        tmp .= x .+ u
        prox!(z, proxg, tmp, μ)
        tmp .= x .- z
        u  .+= tmp

        nxz = norm(tmp)
        if i % printerval == 0
            @printf("%d ||x-z||₂ %.10f\n", i,  nxz)
            if cb != nothing
                cb(x,z)
            end
        end
        if nxz < tol
            @printf("%d ||x-z||₂ %.10f\n", i,  nxz)
            @info("||x-z||₂ ≤ tol")
            break
        end
    end
    x,z
end


text:
using LinearAlgebra
import Base.@kwdef

export steering, beamform
export Bartlett, Capon, Music

abstract type Beamformer end

"""
Frequency-domain Bartlett beamformer.
"""
struct Bartlett <: Beamformer end

@doc """
Frequency-domain Capon beamformer with diagonal loading factor `γ`.
""" Capon
@kwdef struct Capon <: Beamformer
  γ = 0.0
end

@doc """
Frequency-domain MUSIC beamformer with `nsignals` signals.
""" Music
@kwdef struct Music <: Beamformer
  nsignals = 1
end

"""
$(SIGNATURES)
Computes steering delays for specified receiver positions `rxpos`, signal
propagation speed `c`, and angles `θ`. The delays are computed with a far-field
assumption, i.e., for plane incoming waves.

# Examples:
```julia-repl
julia> steering(0.0:1.0:5.0, 1500.0, range(0.0, π; length=181))
6×181 Array{Float64,2}:
  0.00166667    0.00166641   …  -0.00166641   -0.00166667
  0.001         0.000999848     -0.000999848  -0.001
  0.000333333   0.000333283     -0.000333283  -0.000333333
 -0.000333333  -0.000333283      0.000333283   0.000333333
 -0.001        -0.000999848      0.000999848   0.001
 -0.00166667   -0.00166641   …   0.00166641    0.00166667

julia> rxpos = [  # can be 2D or 3D coordinates
  0.0  1.0  2.0  3.0  4.0  5.0
  0.0  0.0  0.0  0.0  0.0  0.0
];
julia> steering(rxpos, 1500.0, range(0.0, π; length=181))
6×181 Array{Float64,2}:
  0.00166667    0.00166641   …  -0.00166641   -0.00166667
  0.001         0.000999848     -0.000999848  -0.001
  0.000333333   0.000333283     -0.000333283  -0.000333333
 -0.000333333  -0.000333283      0.000333283   0.000333333
 -0.001        -0.000999848      0.000999848   0.001
 -0.00166667   -0.00166641   …   0.00166641    0.00166667
```
"""
function steering(rxpos::AbstractMatrix, c, θ)
  # TODO: add support for 2D steering
  nsensors = size(rxpos, 2)
  ndir = length(θ)
  pos0 = sum(rxpos; dims=2) / nsensors
  rxpos = rxpos .- pos0
  size(rxpos, 1) == 1 && (rxpos = vcat(rxpos, zeros(1, nsensors)))
  size(rxpos, 1) > 2 && (rxpos = rxpos[1:2,:])
  sd = zeros(nsensors, ndir)
  for k ∈ 1:ndir
    for j ∈ 1:nsensors
      sd[j,k] = -rxpos[:,j]' * [cos(θ[k]), sin(θ[k])] / c
    end
  end
  sd
end

steering(rxpos::AbstractVector, c, θ) = steering(collect(rxpos'), c, θ)

"""
    beamform(s, sd; fs=framerate(s))

Broadband time-domain delay-and-sum beamformer. Takes in passband or baseband
signals `s` and produces beamformer output for all directions specified by the
steering delays `sd`. The beamformer output is a timeseries signal for each
steering direction.

# Example:
```julia-repl
julia> x = cw(100.0, 1.0, 44100.0);
julia> sd = steering(0.0:1.0:3.0, 1500.0, range(0.0, π; length=181));
julia> bfo = beamform([x x x x], sd)
44100-frame, 181-channel SampleBuf{Complex{Float64}, 2}
1.0s sampled at 44100.0Hz
  ⋮
```
"""
function beamform(s, sd; fs=framerate(s))
  ndir = size(sd, 2)
  bfo = zeros(eltype(s), (nframes(s), ndir))
  for k ∈ 1:ndir
    for j ∈ 1:nchannels(s)
      bfo[:,k] .+= padded(s[:,j], 0; delay=round(Int, -sd[j,k]*fs))
    end
  end
  signal(bfo, fs)
end

"""
    beamform(s, f, n, sd; fs=framerate(s), method=Bartlett())
    beamform(s, f, sd; fs=framerate(s), method=Bartlett())

Narrowband frequency-domain beamformer. Takes in passband signals `s` and
produces beamformer output for all directions specified by the steering
delays `sd`. The beamformer output is an energy estimate (or equivalent)
for each steering direction. The beamforming only uses a narrowband signal
cenetered at frequency `f` with a bandwidth of about `fs/n`.

If `n` is not specified, or is 1, then the input signal is assumed to be
narrowband, and centered at frequency `f`.

The narrowband assumption requires that the bandwidth be no greater than
about 5/T, where T is the maximum time taken for a signal to propagate through
the array.

Several beamforming methods are available:
- [`Bartlett`](@ref)`()`
- [`Capon`](@ref)`(γ)`
- [`Music`](@ref)`(nsignals)`

Custom beamformers can be implemented by creating a subtype of
`SignalAnalysis.Beamformer` and implementing the `SignalAnalysis.beamformer()`
method dispatched on that type.

# Example:
```julia-repl
julia> x = cw(100.0, 1.0, 44100.0);
julia> sd = steering(0.0:1.0:3.0, 1500.0, range(0.0, π; length=181));
julia> bfo = beamform([x x x x], 100.0, 4096, sd; method=Capon(0.1))
181-element Array{Float64,1}:
 0.12406290296318974
 0.1240975045516605
 ⋮
 0.12406290296318974
```
"""
function beamform(s, f, n, sd; fs=framerate(s), method=Bartlett())
  s̃ = n > 1 ? goertzel(s, f, n; fs=fs) : signal(s, fs)
  diff([extrema(abs.(sd))...])[1] > 5/framerate(s̃) && @warn "Narrowband assumption not met, try increasing n"
  R = samples(s̃)' * samples(s̃)
  sv = cis.(2π * f * sd)
  beamformer(method, R, sv)
end

beamform(s, f, sd; fs=framerate(s), method=Bartlett()) = beamform(s, f, 1, sd; fs=fs, method=method)

# TODO: broadband beamformer

### various frequency-domain beamformer implementations
#
# beamformer(opt::type, R, sv)
#   R: narrowband covariance matrix
#   sv: narrowband steering vector
#   opt: beamformer-specific options

function beamformer(::Bartlett, R, sv)
  # TODO: add support for shading
  [abs.(sv[:,j]' * R * sv[:,j]) for j ∈ 1:size(sv,2)]
end

function beamformer(opt::Capon, R, sv)
  opt.γ != 0 && (R .+= opt.γ * I(size(R, 1)))
  [1 ./ abs.(sv[:,j]' * inv(R) * sv[:,j]) for j ∈ 1:size(sv,2)]
end

function beamformer(opt::Music, R, sv)
  F = eigen(R)
  ndx = sortperm(F.values; rev=true)
  Q = F.vectors[:, ndx[opt.nsignals+1:end]]
  Q = Q * Q'
  [1 ./ abs.(sv[:,j]' * Q * sv[:,j]) for j ∈ 1:size(sv,2)]
end


text:
include("LinearSolvers.jl")
@reexport using Gridap.LinearSolvers

include("NonLinearSolvers.jl")
@reexport using Gridap.NonLinearSolvers

include("NLSolvers.jl")
@reexport using Gridap.NLSolvers


text:
using  LinearAlgebra
using Plots

function newtons_method(∇f, H, x, ϵ, k_max)
    k, Δ = 1, Inf
    while abs(Δ) > ϵ && k ≤ k_max
        Δ =  ∇f(x)/abs(H(x))
        x -= Δ
        k += 1
    end
    return x
end

#definicao da funcao
f = x-> -(0.4/((x^2+1))^(1/2))+((x^2+1)^(1/2))*(1-0.4/(x^2+1))-x    #multiplicada por -1 para encontrar o minimo
g = x-> -(-(0.4/((x^2+1))^(1/2))+((x^2+1)^(1/2))*(1-0.4/(x^2+1))-x) #funcao original

#definicao da primeira derivada
df(x) = 1.2*x/((x^2+1)^(3/2))+(1-0.4/(x^2+1))*x/((x^2+1)^(1/2))-1

#definicao da segunda derivada
ddf(x) = (-0.6*x^2+1.8)/((x^2+1)^(5/2))

#converge para o ponto de minimo
initial_point = 1       #ponto inicial
convergence = 1e-8      #convergencia
maximum_iterations = 10 #numero maximo de iteracoes
xs = newtons_method(df,ddf,initial_point,convergence,maximum_iterations)
#xs,f(xs),ddf(xs),df(xs)

println("initial_point = ", initial_point)
println("x* = ", xs)
println("f(x*) = ", g(xs))

plot(g,min(initial_point,xs)-1.,max(initial_point,xs)+1.,lw=2,draw_arrow="true",label="Function")
scatter!([initial_point],[g(initial_point)])
scatter!([xs],[g(xs)])


text:
abstract KernelFilter
type AverageFilter <: KernelFilter
  xband :: Int64
  yband :: Int64
end

function AverageFilter()
  return AverageFilter(1,1)
end

type FourierFilter <: KernelFilter
  xband :: Float64
  yband :: Float64
end
function FourierFilter()
  return FourierFilter(2.0,2.0)
end

function xysmooth( signal :: Array{Float64,3}, filter :: KernelFilter)
  res = zeros(size(signal))
  for iz = 1:size(signal,3)
    res[:,:,iz] = xysmooth(signal[:,:,iz], filter)
  end
  return res
end

function xysmooth( x:: Array{Float64,2}, filter :: AverageFilter)
  y = zeros(size(x))
  xband = filter.xband
  yband = filter.yband
  for ix = 1:size(x,1)
    for iy = 1:size(x,2)
      px :: UnitRange{Int64} = max(1,ix-xband):min(size(x,1),ix + xband)
      py :: UnitRange{Int64} = max(1,iy-yband):min(size(x,2),iy + yband)
      y[ix,iy] = mean(x[px,py])
    end
  end
  return y
end

function xysmooth( ss :: Array{Float64,2}, filter :: FourierFilter)

  xband = filter.xband
  yband = filter.yband
  signal =ss

  kx = collect(1:(size(signal,1)/xband)) / size(signal,1)
  ky = collect(1:(size(signal,2)/yband)) / size(signal,2)

  fresx = zeros(size(signal))
  for ix = kx
    psix = ix*pi*collect(1:size(signal,1))
    fsinx  = mean(sin(psix) .* signal,1).* sin(psix)/(0.5)
    fcosx  = mean(cos(psix) .* signal,1).* cos(psix)/(0.5)
    fresx += fsinx + fcosx
  end
  signal = fresx

  fresy = zeros(size(signal))
  for iy = ky
    psiy = iy*pi*collect(1:size(signal,2))
    psiy = psiy'
    fsiny  = mean(sin(psiy) .* signal,2).* sin(psiy)/(0.5)
    fcosy  = mean(cos(psiy) .* signal,2).* cos(psiy)/(0.5)
    fresy += fsiny + fcosy
  end

  return fresy
end

export xysmooth
export AverageFilter, FourierFilter
export KernelFilter


text:
@testset "Node composition plots" begin

nodeCompositionPlot( mpSim, 12, "Master", "Junior", "Career", "Branch A",
    "Branch B" )
nodeCompositionPlot( mpSim, 12, "Junior", "Career", "Branch A", "Branch B",
    plotType=:stacked )
#=
nodeCompositionPlot( mpSim, 12, "Master", "Junior", "Career", "Branch A",
    "Branch B", showPlot=false, filename="plot/node composition test 1" )
nodeCompositionPlot( mpSim, 12, "Junior", "Career", "Branch A", "Branch B",
    showPlot=false, plotType=:stacked, filename="plot/node composition test 2" )
=#
end  # @testset "Node composition plots"


text:
function PrettyPrint.pprint_impl(io, v::LocalVar, indent, newline)
    print(io, "LocalVar($(v.sym), is_mutable=$(v.is_mutable.x), is_shared=$(v.is_shared.x))")
end

function PrettyPrint.pprint_impl(io, p::Pair, indent, newline)
    pprint(io, p.first, indent, newline)
    print(io, "=>")
    pprint(io, p.second, indent + 2, false)
end

function PrettyPrint.pprint_impl(io, d::Dict, indent, newline)
   PrettyPrint.pprint_for_seq(io, "{", "}", collect(d), indent, newline)
end

Base.show(io::IO, scope :: Scope) = pprint(io, scope)


text:
export AbstractTrajectory

"""
    AbstractTrajectory

A trace is used to record some useful information
during the interactions between agents and environments.

Required Methods:

- `Base.haskey(t::AbstractTrajectory, s::Symbol)`
- `Base.getindex(t::AbstractTrajectory, s::Symbol)`
- `Base.keys(t::AbstractTrajectory)`
- `Base.push!(t::AbstractTrajectory, kv::Pair{Symbol})`
- `Base.pop!(t::AbstractTrajectory, s::Symbol)`
- `Base.empty!(t::AbstractTrajectory)`

Optional Methods:

- `isfull`

"""
abstract type AbstractTrajectory end

function Base.push!(t::AbstractTrajectory; kwargs...)
    for kv in kwargs
        push!(t, kv)
    end
end

"""
    Base.pop!(t::AbstractTrajectory, s::Symbol...)

`pop!` out one element of the traces specified in `s`
"""
function Base.pop!(t::AbstractTrajectory, s::Tuple{Vararg{Symbol}})
    NamedTuple{s}(pop!(t, x) for x in s)
end

Base.pop!(t::AbstractTrajectory) = pop!(t, keys(t))

function Base.empty!(t::AbstractTrajectory)
    for s in keys(t)
        empty!(t[s])
    end
end

#####
# patch code
#####

# avoid showing the inner structure
function AbstractTrees.children(t::StructTree{<:AbstractTrajectory})
    Tuple(k => StructTree(t.x[k]) for k in keys(t.x))
end

@deprecate get_trace(t::AbstractTrajectory, s::Symbol) t[s]


text:
using Genie, Genie.Renderer, Genie.Renderer.Html

@test html(path"formview.jl.html").body |> String === """<form method="POST" enctype="multipart/form-data" action="/new">
Example select<select class="form-control" id="exampleFormControlSelect1"><option>1</option><option>2</option><option>3</option><option>4</option><option>5</option></select>
<button class="btn btn-primary" type="submit">Submit</button></form>"""


text:
module GLM_

# -------------------------------------------------------------------
# TODO
# - return feature names in the report
# - return feature importance curve to report using `features`
# - handle binomial case properly, needs MLJ API change for weighted
# samples (y/N ~ Be(p) with weights N)
# - handle levels properly (see GLM.jl/issues/240); if feed something
# with levels, the fit will fail.
# - revisit and test Poisson and Negbin regression once there's a clear
# example we can test on (requires handling levels which deps upon GLM)
# - test Logit, Probit etc on Binomial once binomial case is handled
# -------------------------------------------------------------------

import MLJBase
import MLJBase: metadata_pkg, metadata_model
import Distributions
using Parameters
using Tables

import ..GLM

export LinearRegressor, LinearBinaryClassifier, LinearCountRegressor

##
## DESCRIPTIONS
##

const LR_DESCR = "Linear regressor (OLS) with a Normal model."
const LBC_DESCR = "Linear binary classifier with specified link (e.g. logistic)."
const LCR_DESCR = "Linear count regressor with specified link and distribution (e.g. log link and poisson)."

###
## Helper functions
###

"""
augment_X(X, b)

Augment the matrix `X` with a column of ones if the intercept is to be fitted (`b=true`), return
`X` otherwise.
"""
function augment_X(X::Matrix, b::Bool)::Matrix
    b && return hcat(X, ones(eltype(X), size(X, 1), 1))
    return X
end

"""
glm_report(fitresult)

Report based on the `fitresult` of a GLM model.
"""
glm_report(fitresult) = ( deviance     = GLM.deviance(fitresult),
                          dof_residual = GLM.dof_residual(fitresult),
                          stderror     = GLM.stderror(fitresult),
                          vcov         = GLM.vcov(fitresult) )

####
#### REGRESSION TYPES
####

# LinearRegressor        --> Probabilistic w Continuous Target
# LinearCountRegressor   --> Probabilistic w Count Target
# LinearBinaryClassifier --> Probabilistic w Binary target // logit,cauchit,..
# MulticlassClassifier   --> Probabilistic w Multiclass target

@with_kw_noshow mutable struct LinearRegressor <: MLJBase.Probabilistic
    fit_intercept::Bool      = true
    allowrankdeficient::Bool = false
end

@with_kw_noshow mutable struct LinearBinaryClassifier{L<:GLM.Link01} <: MLJBase.Probabilistic
    fit_intercept::Bool = true
    link::L             = GLM.LogitLink()
end

@with_kw_noshow mutable struct LinearCountRegressor{D<:Distributions.Distribution,L<:GLM.Link} <: MLJBase.Probabilistic
    fit_intercept::Bool = true
    distribution::D     = Distributions.Poisson()
    link::L 			= GLM.LogLink()
end

# Short names for convenience here

const GLM_MODELS = Union{<:LinearRegressor, <:LinearBinaryClassifier, <:LinearCountRegressor}

####
#### FIT FUNCTIONS
####

function MLJBase.fit(model::LinearRegressor, verbosity::Int, X, y)
    # apply the model
    features  = Tables.schema(X).names
    Xmatrix   = augment_X(MLJBase.matrix(X), model.fit_intercept)
    fitresult = GLM.glm(Xmatrix, y, Distributions.Normal(), GLM.IdentityLink())
    # form the report
    report    = glm_report(fitresult)
    cache     = nothing
    # return
    return fitresult, cache, report
end

function MLJBase.fit(model::LinearCountRegressor, verbosity::Int, X, y)
    # apply the model
    features  = Tables.schema(X).names
    Xmatrix   = augment_X(MLJBase.matrix(X), model.fit_intercept)
    fitresult = GLM.glm(Xmatrix, y, model.distribution, model.link)
    # form the report
    report    = glm_report(fitresult)
    cache     = nothing
    # return
    return fitresult, cache, report
end

function MLJBase.fit(model::LinearBinaryClassifier, verbosity::Int, X, y)
    # apply the model
    features  = Tables.schema(X).names
    Xmatrix   = augment_X(MLJBase.matrix(X), model.fit_intercept)
    decode    = y[1]
    y_plain   = MLJBase.int(y) .- 1 # 0, 1 of type Int
    fitresult = GLM.glm(Xmatrix, y_plain, Distributions.Bernoulli(), model.link)
    # form the report
    report    = glm_report(fitresult)
    cache     = nothing
    # return
    return (fitresult, decode), cache, report
end

function MLJBase.fitted_params(model::GLM_MODELS, fitresult)
    coefs = GLM.coef(fitresult)
    return (coef      = coefs[1:end-Int(model.fit_intercept)],
            intercept = ifelse(model.fit_intercept, coefs[end], nothing))
end

####
#### PREDICT FUNCTIONS
####

# more efficient than MLJBase fallback
function MLJBase.predict_mean(model::Union{LinearRegressor,<:LinearCountRegressor}, fitresult, Xnew)
    Xmatrix = augment_X(MLJBase.matrix(Xnew), model.fit_intercept)
    return GLM.predict(fitresult, Xmatrix)
end

function MLJBase.predict_mean(model::LinearBinaryClassifier, (fitresult, _), Xnew)
    Xmatrix = augment_X(MLJBase.matrix(Xnew), model.fit_intercept)
    return GLM.predict(fitresult, Xmatrix)
end

function MLJBase.predict(model::LinearRegressor, fitresult, Xnew)
    μ = MLJBase.predict_mean(model, fitresult, Xnew)
    σ̂ = GLM.dispersion(fitresult)
    return [GLM.Normal(μᵢ, σ̂) for μᵢ ∈ μ]
end

function MLJBase.predict(model::LinearCountRegressor, fitresult, Xnew)
    λ = MLJBase.predict_mean(model, fitresult, Xnew)
    return [GLM.Poisson(λᵢ) for λᵢ ∈ λ]
end

function MLJBase.predict(model::LinearBinaryClassifier, (fitresult, decode), Xnew)
    π = MLJBase.predict_mean(model, (fitresult, decode), Xnew)
    return [MLJBase.UnivariateFinite(MLJBase.classes(decode), [1-πᵢ, πᵢ]) for πᵢ in π]
end

# NOTE: predict_mode uses MLJBase's fallback

####
#### METADATA
####

# shared metadata
const GLM_REGS = Union{Type{<:LinearRegressor}, Type{<:LinearBinaryClassifier}, Type{<:LinearCountRegressor}}

metadata_pkg.((LinearRegressor, LinearBinaryClassifier, LinearCountRegressor),
    name="GLM",
    uuid="38e38edf-8417-5370-95a0-9cbb8c7f171a",
    url="https://github.com/JuliaStats/GLM.jl",
    julia=true,
    license="MIT",
    is_wrapper=false
    )

metadata_model(LinearRegressor,
    input=MLJBase.Table(MLJBase.Continuous),
    output=AbstractVector{MLJBase.Continuous},
    weights=false,
    descr=LR_DESCR
    )

metadata_model(LinearBinaryClassifier,
    input=MLJBase.Table(MLJBase.Continuous),
    target=AbstractVector{MLJBase.UnivariateFinite},
    weights=false,
    descr=LBC_DESCR
    )

metadata_model(LinearCountRegressor,
    input=MLJBase.Table(MLJBase.Continuous),
    target=AbstractVector{MLJBase.Count},
    weights=false,
    descr=LCR_DESCR
    )

end # module


text:
using BenchmarkTools

function benchmark(f, args...)
    for tag in (:Direct, :DataStructures, :BinaryHeaps)
        if applicable(f, Val(tag), args...)
            title = "$tag implementation"
            println(repeat("-", length(title)+1))
            println("$title:")
            println(repeat("-", length(title)+1))
            println()
            display(@benchmark $f($(Val(tag)),$args...))
            println()
            println()
        end
    end
end

module WalkDown
    function run(::Val{:Direct}, n)
        i = 1
        s = 0
        while heapleft(i) <= n
            s += i
            i = heapleft(i)
        end
        s += i
        return s
    end

    using BinaryHeaps: heapleft, walk
    function run(::Val{:BinaryHeaps}, n)
        s = Ref(0)
        walk(n, 1, Val(:down)) do i,np,nc
            s[] += i
            return Val(:left)
        end
        return s[]
    end
end

module Heapify
    using Base.Order

    using DataStructures
    run(::Val{:DataStructures}, v) = DataStructures.heapify!(v, Forward)

    using BinaryHeaps
    run(::Val{:BinaryHeaps}, v) = BinaryHeaps.heapify!(v, Forward)
end


text:
msg = "Rosetta code"

using Nettle
digest = hexdigest("sha256", msg)

# native
using SHA
digest1 = join(num2hex.(sha256(msg)))

@assert digest == digest1


text:
module AffineMapsTests

using Gridap.TensorValues
using Gridap.Arrays
using Gridap.Fields

origin = Point(1,1)
jacobian = TensorValue(2,0,0,2)

h = AffineMap(jacobian,origin)

x1 = Point(0,0)
x2 = Point(1,1)
x3 = Point(2,1)
x = [x1,x2,x3]

r = Point{2,Int}[(1, 1), (3, 3), (5, 3)]
∇r = TensorValue{2,2,Int,4}[(2, 0, 0, 2), (2, 0, 0, 2), (2, 0, 0, 2)]
test_field(h,x,r,grad=∇r)

end # module


text:
using IndirectArrays, MappedArrays
using Test, FixedPointNumbers, Colors

colors = [RGB(1,0,0) RGB(0,1,0);
          RGB(0,0,1) RGB(1,0,0)]
index0 = [1 3;
          2 1]
for indexT in (Int8, Int16, UInt8, UInt16)
    A = IndirectArray{indexT}(colors)
    @test eltype(A) == RGB{N0f8}
    @test size(A) == (2,2)
    @test ndims(A) == 2
    @test A[1,1] === A[1] === RGB(1,0,0)
    @test A[2,1] === A[2] === RGB(0,0,1)
    @test A[1,2] === A[3] === RGB(0,1,0)
    @test A[2,2] === A[4] === RGB(1,0,0)
    @test isa(eachindex(A), AbstractUnitRange)
    @test A.index == index0
end
x = IndirectArray(colors[:])
@test x == IndirectArray{UInt8}(colors[:])
xc = copy(x)
@test x == xc
xc[2], xc[3] = RGB(0,1,0), RGB(0,1,1)
@test xc == IndirectArray([RGB(1,0,0), RGB(0,1,0), RGB(0,1,1), RGB(1,0,0)])
@test append!(x, x) == IndirectArray([colors[:]; colors[:]])
@test append!(x, IndirectArray([RGB(1,0,0), RGB(1,1,0)])) ==
    IndirectArray([colors[:]; colors[:]; [RGB(1,0,0), RGB(1,1,0)]])
# Append with non-IndirectArray
@test append!(IndirectArray(colors[:]), IndirectArray(colors[:])) ==
      append!(IndirectArray(colors[:]), colors[:])

# Bounds checking upon construction
index_ob = copy(index0)
index_ob[1] = 5   # out-of-bounds
unsafe_ia(idx, vals) = (@inbounds ret = IndirectArray(idx, vals); ret)
  safe_ia(idx, vals) = (ret = IndirectArray(idx, vals); ret)
@test_throws BoundsError safe_ia(index_ob, colors[1:3])
# This requires inlining, which means it fails on Travis since we turn
# off inlining for better coverage stats
# B = unsafe_ia(index_ob, colors)
# @test_throws BoundsError B[1]
# @test B[2] == RGB(0,0,1)

# Non-Arrays
a = [0.1  0.4;
     0.33 1.0]
f(x) = round(Int, 99*x) + 1   # maps 0-1 to 1-100
m = mappedarray(f, a)
cmap = colormap("RdBu", 100)
img = IndirectArray(m, cmap)
@test img == [cmap[11] cmap[41];
              cmap[34] cmap[100]]


text:
# This file includes the theoretical functions of the ber performance cureves.

export berask, berpam, berpsk, berfsk, berqam, plotber

"""
    $SIGNATURES

Returns the probability of symbol error for the snr per bit `γb` and constallation size `M` for ASK and PAM signalling
"""
berask(γb, M) = 2 * (1 - 1 / M) * Q(sqrt(6log2(M) / (M^2 - 1) * dbtoval(γb)))
berpam = berask

"""
    $SIGNATURES

Returns the probability of symbol error for the snr per bit `γb` and constallation size `M` for PSK signalling
"""
berpsk(γs, M) = 1 - quadgk(θ -> pθ(θ, dbtoval(γs)), -π/M, π/M)[1]
pθ(θ, γs) = 1 / (2π) * exp(-γs * sin(θ) * sin(θ)) * quadgk(υ -> integrand(υ, γs, θ), 0, Inf)[1]
integrand(υ, γs, θ) = υ * exp(-1 / 2 * (υ - sqrt(2 * γs) * cos(θ))^2)

"""
    $SIGNATURES

Returns the probability of symbol error for the snr per bit `γb` and constallation size `M` for FSK signalling
"""
function berfsk(γs, M)
    k = log2(M)
    pe = 1 / sqrt(2π) * quadgk(x -> (1 - (1  - Q(x))^(M - 1)) * exp(-1/2 * (x - sqrt(2 * dbtoval(γs)))^2), -Inf, Inf)[1]
    # # Return probability of symbol error, not probability of bits error
    # pb = 2^(k -1) / (2^k - 1) * pe
end

"""
    $SIGNATURES

Returns the probability of symbol error for the snr per bit `γb` and constallation size `M` for QAM signalling
"""
function berqam(γb, M)
    qval = Q(√(3log2(M) / (M - 1) * dbtoval(γb)))
    4 * (1 - 1 / √(M)) * qval * (1 - (1 - 1 / √(M)) * qval)
end

#------------------------ Plots ---------------------------------------
"""
    $SIGNATURES

Plots the probability of symbol error versus snr per bit for the signaling `scheme`. `snr_per_bit_range` is the snr per bit range and `krange` is the symbol size of the M-ary signalling where ``M=2^k``.
"""
function plotber(;scheme="ASK", snr_per_bit_range=collect(-4 : 1 : 20), krange = 1 : 4,
    pltkwargs...)
    # Extract the ber function
    schemes = ["ASK", "PSK", "QAM", "FSK"]
    scheme in schemes || throw(ArgumentError("Expected `$schemes`, got `$scheme` instead"))
    @eval pe_vs_snr = $(Symbol("ber" * lowercase(scheme)))

    # Plot
    plt = plot(xlabel="γb", ylabel="Pe", title="M-$(scheme)", pltkwargs...)
    for k in krange
        M = 2^k
        γs = k * dbtoval.(snr_per_bit_range)
        plot!(snr_per_bit_range,  pe_vs_snr.(γs, M), yscale=:log10, label="M=$M")
    end
    plt
end


text:
"""
    Flowpipe{RSN<:AbstractReachSet} <: AbstractSolution

Wrapper of a sequence of sets (the solution of a continuous-post operator).

### Fields

- `reachsets` -- the list of [`AbstractReachSet`](@ref)s
"""
struct Flowpipe{RSN<:AbstractReachSet}
    reachsets::Vector{RSN}
end

# projection
function project(fp::Flowpipe, M::AbstractMatrix)
    return Flowpipe([project(X, M) for X in fp.reachsets])
end