text:
function uwdf!( S::GphysData,
                fname::String,
                full::Bool,
                memmap::Bool,
                strict::Bool,
                v::Integer)
  D = Dict{String,Any}()

  # Open data file
  fid = memmap ? IOBuffer(Mmap.mmap(fname)) : open(fname, "r")

  # Process main header
  N             = Int64(bswap(fastread(fid, Int16)))
  mast_fs       = bswap(fastread(fid, Int32))
  mast_lmin     = bswap(fastread(fid, Int32))
  mast_lsec     = bswap(fastread(fid, Int32))
  mast_nx       = bswap(fastread(fid, Int32))
  fastskip(fid, 24)
  extra         = fastread(fid, 10)
  fastskip(fid, 80)

  if v > 2
    println("mast_header:")
    println("N = ", N)
    println("mast_fs = ", mast_fs)
    println("mast_lmin = ", mast_lmin)
    println("mast_lsec = ", mast_lsec)
    println("mast_nx = ", mast_nx)
    println("extra = ", Char.(extra))
  end

  # Seek EOF to get number of structures
  fastseekend(fid)
  fastskip(fid, -4)
  nstructs = bswap(fastread(fid, Int32))
  v>0 && println(stdout, "nstructs = ", nstructs)
  structs_os = (-12*nstructs)-4
  tc_os = 0
  v>1 && println(stdout, "structs_os = ", structs_os)

  # Set version of UW seismic data file (char may be empty, leave code as-is!)
  uw2::Bool = extra[3] == 0x32 ? true : false
  chno = Array{Int32, 1}(undef, N)
  corr = Array{Int32, 1}(undef, N)

  # Read in UW2 data structures from record end
  if uw2
    fastseekend(fid)
    fastskip(fid, structs_os)
    for j = 1:nstructs
      structtag     = fastread(fid)
      fastskip(fid, 3)
      M             = bswap(fastread(fid, Int32))
      byteoffset    = bswap(fastread(fid, Int32))
      if structtag == 0x43 # 'C'
        N = Int64(M)
      elseif structtag == 0x54 # 'T'
        fpos        = fastpos(fid)
        fastseek(fid, byteoffset)
        chno = Array{Int32, 1}(undef, M)
        corr = Array{Int32, 1}(undef, M)
        n = 0
        @inbounds while n < M
          n += 1
          chno[n]   = fastread(fid, Int32)
          corr[n]   = fastread(fid, Int32)
        end
        chno .= (bswap.(chno) .+ 1)
        corr .= bswap.(corr)
        tc_os = -8*M
        fastseek(fid, fpos)
      end
    end
  end
  v>0 && println(stdout, "Processing ", N , " channels.")

  # Write time corrections
  timecorr = zeros(Int64, N)
  if length(chno) > 0
    for n = 1:N
      # corr is in μs
      timecorr[chno[n]] = Int64(corr[n])
    end
  end

  # Read UW2 channel headers ========================================
  if uw2
    fastseekend(fid)
    fastskip(fid, -56*N + structs_os + tc_os)
    I32 = Array{Int32, 2}(undef, 5, N)    # chlen, offset, lmin, lsec (μs), fs,   unused: expan1
    I16 = Array{Int16, 2}(undef, 3, N)    # lta, trig, bias                       unused: fill
    U8  = Array{UInt8, 2}(undef, 24, N)   # name(8), tmp(4), compflg(4), chid(4), expan2(4)
    i = 0
    @inbounds while i < N
      i = i + 1
      I32[1,i] = bswap(fastread(fid, Int32))
      I32[2,i] = bswap(fastread(fid, Int32))
      I32[3,i] = bswap(fastread(fid, Int32))
      I32[4,i] = bswap(fastread(fid, Int32))
      I32[5,i] = bswap(fastread(fid, Int32))
      if full == true
        fastskip(fid, 4)
        I16[1,i] = bswap(fastread(fid, Int16))
        I16[2,i] = bswap(fastread(fid, Int16))
        I16[3,i] = bswap(fastread(fid, Int16))
        fastskip(fid, 2)
      else
        fastskip(fid, 12)
      end
      j = 0
      while j < 24
        j = j + 1
        U8[j,i] = fastread(fid)
      end
    end

    # Parse U8 --------------------------------------------
    # rows 01:08    channel name
    # rows 09:12    format code
    # rows 13:16    compflg(4)
    # rows 17:20    chid
    # rows 21:24    expan2
    if full == true
      j = 0
      while j < N
        j = j + 1
        i = 16
        while i < 24
          i = i + 1
          if getindex(U8, i, j) == 0x00
            setindex!(U8, 0x20, i, j)
          end
        end
      end
    end

    fastseek(fid, getindex(I32, 2, 1))
    buf = getfield(BUF, :buf)
    checkbuf_8!(buf, 4*maximum(I32[1,:]))
    id = BUF.id
    id[1] = 0x55
    id[2] = 0x57
    id[3] = 0x2e

    i = 0
    os = 0
    @inbounds while i < N
      i += 1
      fastskip(fid, os)
      nx = getindex(I32, 1, i)

      # Generate ID
      j = 3
      k = 1
      while j < 8 && k < 9
        c = getindex(U8, k, i)
        if c != 0x00
          j += 1
          id[j] = c
        end
        k += 1
      end
      id[j+1] = 0x2e
      id[j+2] = 0x2e
      j += 2
      J = j+3
      k = 13
      while j < J && k < 17
        c = getindex(U8, k, i)
        if c != 0x00
          j += 1
          id[j] = c
        end
        k += 1
      end

      # Save to SeisChannel
      C = SeisChannel()
      setfield!(C, :id, unsafe_string(pointer(id), j))
      setfield!(C, :fs, Float64(getindex(I32, 5, i))*1.0e-3)
      setfield!(C, :units, "m/s")
      if full == true
        D = getfield(C, :misc)
        D["lta"]    = I16[1,i]
        D["trig"]   = I16[2,i]
        D["bias"]   = I16[3,i]
        D["chid"]   = String(getindex(U8, 17:20, i))
        D["expan2"] = String(getindex(U8, 21:24, i))
        if i == 1
          D["mast_fs"]    = mast_fs*1.0f-3
          D["mast_lmin"]  = mast_lmin
          D["mast_lsec"]  = mast_lsec
          D["mast_nx"]    = mast_nx
          D["extra"]      = String(extra)

          # Go back to main header; grab what we skipped
          p = fastpos(fid)
          fastseek(fid, 18)         # we have the first few fields already
          D["mast_tape_no"]   = bswap(fastread(fid, Int16))
          D["mast_event_no"]  = bswap(fastread(fid, Int16))
          D["flags"]          = bswap.(fastread!(fid, Array{Int16, 1}(undef, 10)))
          fastskip(fid, 10)         # we have "extra" already
          comment             = fastread(fid, 80)
          D["comment"]        = String(comment[comment.!=0x00])

          # Return to where we were
          fastseek(fid, p)
        end
      end

      # Generate T
      t = Array{Int64,2}(undef,2,2)
      ch_time = 60000000*Int64(getindex(I32, 3, i)) +
                Int64(getindex(I32, 4, i)) +
                getindex(timecorr, i) -
                11676096000000000
      mk_t!(C, nx, ch_time)

      # Generate X
      x = Array{Float32,1}(undef, nx)
      fmt = getindex(U8, 9, i)
      if fmt == 0x53
        fast_readbytes!(fid, buf, 2*nx)
        fillx_i16_be!(x, buf, nx, 0)
      elseif fmt == 0x4c
        fast_readbytes!(fid, buf, 4*nx)
        fillx_i32_be!(x, buf, nx, 0)
      else
        fast_readbytes!(fid, buf, 4*nx)
        x .= bswap.(reinterpret(Float32, buf))[1:nx]
      end
      setfield!(C, :x, x)

      # Push to SeisData
      add_chan!(S, C, strict)

      if i < N
        os = getindex(I32, 2, i+1) - fastpos(fid)
      end
    end
  end
  close(fid)
  return nothing
end

function uwdf(fname::String, full::Bool, memmap::Bool, strict::Bool, v::Integer)
  S = SeisData()
  uwdf!(S, fname, full, memmap, strict, v)
  return S
end


text:
function solve_002()
    fib_vec = [1, 2]
    i = 3
    # Generate fibonacci numbers
    while fib_vec[end] < 4000000
        push!(fib_vec, fib_vec[i-1] + fib_vec[i-2])
        i += 1
    end
    # Sum even terms
    sum(filter(x -> x % 2 == 0, fib_vec))
end

ext: jl


text:
import MKL # must load before anything else
using SafeTestsets: @safetestset

@time @safetestset "util" begin include("util.jl") end
@time @safetestset "conv" begin include("conv.jl") end
@time @safetestset "categorical" begin include("categorical.jl") end
@time @safetestset "convrbm" begin include("convrbm.jl") end
@time @safetestset "pooling" begin include("pooling.jl") end
@time @safetestset "denseconvrbm" begin include("denseconvrbm.jl") end
@time @safetestset "hankel" begin include("hankel.jl") end
@time @safetestset "free energy" begin include("free_energy.jl") end


text:
## This also tests some vanilla pricing engines ##
using Base.Test
using QuantLib

cal = QuantLib.Time.TargetCalendar()
todaysDate = Date(1998, 5, 15)
settlementDate = Date(1998, 5, 17)

set_eval_date!(settings, todaysDate)

optionType = Put()
underlying = 36.0
strike = 40.0
dividendYield = 0.00
riskFreeRate = 0.06
vol = 0.20
mat = Date(1999, 5, 17)
dc = QuantLib.Time.Actual365()

exerciseDates = Date[settlementDate + Dates.Month(3 * i) for i = 1:4]

europeanExercise = EuropeanExercise(mat)
bermudanExercise = BermudanExercise(exerciseDates)
americanExercise = AmericanExercise(settlementDate, mat)

underlyingH = Quote(underlying)

flatTermStructure = FlatForwardTermStructure(settlementDate, cal, Quote(riskFreeRate), dc)
flatDividendTS = FlatForwardTermStructure(settlementDate, cal, Quote(dividendYield), dc)
flatVolTS = BlackConstantVol(settlementDate, cal, vol, dc)

payoff = PlainVanillaPayoff(optionType, strike)

bsmProcess = BlackScholesMertonProcess(underlyingH, flatTermStructure, flatDividendTS, flatVolTS)

# test a European Option
bsPE = AnalyticEuropeanEngine(bsmProcess) # black scholes PE
europeanOption = VanillaOption(payoff, europeanExercise, bsPE)

# American Option
baroneAdesiPE = BaroneAdesiWhaleyApproximationEngine(bsmProcess)
americanOption = VanillaOption(payoff, americanExercise, baroneAdesiPE)

# Bermudan Option
timeSteps = 801
fdBermudanPE = FDBermudanEngine(bsmProcess, CrankNelson, timeSteps, timeSteps - 1)
bermudanOption = VanillaOption(payoff, bermudanExercise, fdBermudanPE)

@test npv(europeanOption) == 3.8443077915968398
@test npv(americanOption) == 4.459627613776478
@test npv(bermudanOption) == 4.3608068625088645


text:
export Graph, graph
export predict
export forward
export backward
export update!
export zerograds!
export clear

"""
    mutable struct Graph
# Field
`backward`: store backward operations
"""
mutable struct Graph
    # 存储反向传播操作
    backward::Vector
    function Graph()
        new(Vector(undef,0))
    end
end

import Base.size
import Base.length
import Base.getindex
import Base.deepcopy

Base.size(g::Graph)           = size(g.backward)
Base.size(g::Graph, dim::Int) = size(g.backward, dim)
Base.length(g::Graph)         = length(g.backward)
Base.getindex(g::Graph, k...) =  g.backward[k...]

# -- 全局图变量，存储所有反向运算的中间
# -- 变量，但是在参数更新后需要将其置空
global graph = Graph()

function clear()
    graph.backward = [];
end

function backward()
    for i = length(graph):-1:1
        graph[i]()
    end
    clear()
end

max_stars_repo_name: sonosole/Delta

max_issues_repo_name: sonosole/Delta

max_forks_repo_name: sonosole/Delta


text:
struct Pendulum{S, A, R} <: ContEnv{S, A, R}
    state::Vector{S}
    current::Vector{S}
    observationspace::ObservationSpace{Vector{S}}
    actionspace::ActionSpace{Tuple{A, A}}
    function Pendulum{T}() where T <: AbstractFloat
        observationspace = ObservationSpace{Vector{T}}(3)
        actionspace = ActionSpace([(-one(T), one(T))])
        state = [rand(Uniform(-π, π)), rand(Uniform(-1.0, 1.0))]
        current = [cos(state[1]), sin(state[1]), state[2]]
        new{T, T, T}(state, current, observationspace, actionspace)
    end
end

Pendulum() = Pendulum{Float32}()

function obs!(env::Pendulum)
    env.current[1] = cos(env.state[1])
    env.current[2] = sin(env.state[1])
    env.current[3] = env.state[2]
end

step!(env::Pendulum{S, A, R}, action::Vector{A}) where {S, A, R} = step!(env, action[1])
function step!(env::Pendulum{S, A, R}, action::A) where {S, A, R}
    θ, θ̇ = env.state
    u = clamp(2.0 * action, -2.0, 2.0)
    costs = (((θ + π) % (2.0 * π)) - π) ^ 2.0 + 0.1 * (θ̇ ^ 2.0) + 0.001 * (u ^ 2.0)
    θ̇ += (-15.0 * sin(θ + π) + 3.0 * u) * 0.05
    θ += θ̇ * 0.05
    θ̇ = clamp(θ̇, -8.0, 8.0)
    env.state .= θ, θ̇
    obs!(env)
    return env.current, -R(costs), false
end

function reset!(env::Pendulum)
    env.state .= rand(Uniform(-π, π)), rand(Uniform(-1.0, 1.0))
    obs!(env)
    return env.current
end


text:
struct UniformBinning <: AbstractBinningAlgorithm
    nbins::Int
end

function perform(binning::UniformBinning, predictions::AbstractMatrix{<:Real},
                 labels::AbstractVector{<:Integer})
    # create empty dictionary of bins
    T = eltype(predictions)
    bins = Dict{Int,Bin{T}}()

    # bin all samples
    nclasses = size(predictions, 1)
    @inbounds for i in 1:length(labels)
        prediction = view(predictions, :, i)
        bin = get!(bins, binid(binning, prediction), Bin(T, nclasses))
        adddata!(bin, prediction, labels[i])
    end

    bins
end

function binid(alg::UniformBinning, prediction::AbstractVector{<:Real})
    id = 1
    @inbounds for p in prediction
        id *= binid(alg, p)
    end
    id
end

function binid(alg::UniformBinning, p::Real)
    # check argument
    (p < zero(p) || p > one(p)) &&
        throw(ArgumentError("predictions must be between 0 and 1"))

    # handle special case p = 0
    iszero(p) && return 1

    ceil(alg.nbins * p)
end


text:
using Plots
using Statistics
using ProgressBars
L=101
function BoundaryCondition(i)
    if i==L+1
        return 1
    elseif i==0
        return L
    else
        return i
    end
end
function height_cal(i,arr) #calculating the height of the layer in a column
    i=BoundaryCondition(i)
    row=L #starting from the top layer
    height=L
    for j in 1:L
        if arr[row,i]==0.0 #when the entry is empty
            height-=1 #decreasing height by one
            row-=1
        else
            break #exit the loop for the next column
        end
    end
    return height
end
function neighbor_checking(arr, col)
    h=[]
    for i in (col-1):1:(col+1)
        height= height_cal(i,arr)
        if i==col
            push!(h, (height+1))
        else
            push!(h, height)
        end
    end
    return maximum(h)
end
function width_cal(arr, L)
    firstx=0
    firsty=0
    col=1 #starting from the beggining
    for i in 1:L
        for row in 1:L
            if arr[row, col]!=0.0
                firstx=col
                firsty=row
                break
            end
        end
        if firstx!=0 && firsty!=0
            break
        else
            col+=1
        end
    end
    lastx=0
    lasty=0
    col=L #starting from the end
    for i in 1:L
        for row in 1:L
            if arr[row, col]!=0.0
                lastx=col
                lasty=row
                break
            end
        end
        if lastx!=0 && lasty!=0
            break
        else
            col-=1
        end
    end
    d=sqrt((firstx-lastx)^2 + (firsty-lasty)^2)
    return d
end
function deposing(L,n)
    WidthList=[]
    color=1
    arr=zeros((L,L))
    meanList=[]
    stdList=[]
    finalP=n
    anim = Animation()
    for particle in ProgressBar(1:n)
        col=rand(1:L)
        height=neighbor_checking(arr, col)
        if height>=L
            finalP=particle
            break
        end
        if height!=1 || col==ceil(Int, L/2)
            arr[height, col]=color
        end
        if particle%(10*L*color)==0
            color+=1
        end
        WidthNum=width_cal(arr, L)
        push!(WidthList, WidthNum)
        p = heatmap(hcat(arr), c=cgrad(:roma, 10, categorical = true, scale = :exp), xlabel="L")
        frame(anim, p)
    end
    gif(anim,"C:\\Users\\Narges\\Documents\\GitHub\\computational_physics\\chapter3\\Fig\\3.4.gif",fps=50)

    return arr, WidthList, finalP
end
arr, WidthList, finalP=deposing(L,15000)
heatmap(hcat(arr), c=cgrad(:roma, 10, categorical = true, scale = :exp), xlabel="L")
savefig("C:\\Users\\Narges\\Documents\\GitHub\\computational_physics\\chapter3\\Fig\\3.4.png")
scatter(1:finalP, WidthList, xlabel="time",ylabel="width of the shrub" , markersize=0.001)
savefig("C:\\Users\\Narges\\Documents\\GitHub\\computational_physics\\chapter3\\Fig\\3.4_width.png")


text:
@testset "Degree" begin
    g5 = DiGraph(4)
    add_edge!(g5,1,2); add_edge!(g5,2,3); add_edge!(g5,1,3); add_edge!(g5,3,4)
    for g in testdigraphs(g5)
      @test @inferred(degree_centrality(g)) == [0.6666666666666666, 0.6666666666666666, 1.0, 0.3333333333333333]
      @test @inferred(indegree_centrality(g, normalize=false)) == [0.0, 1.0, 2.0, 1.0]
      @test @inferred(outdegree_centrality(g; normalize=false)) == [2.0, 1.0, 1.0, 0.0]
    end
end

max_stars_repo_name: AsileBcd/JohnsonCycles

max_issues_repo_name: AsileBcd/JohnsonCycles

max_forks_repo_name: AsileBcd/JohnsonCycles


text:
# support literal matches and regpexps

# p"..." creates a matcher for the given regular expression
macro p_str(s)
    Pattern(Regex(s))
end

macro P_str(s)
    Drop(Pattern(Regex(s)))
end

# e"..." creates a matcher for the given string
macro e_str(s)
    Equal(s)
end

macro E_str(s)
    Drop(Equal(s))
end

Opt(m::Matcher) = Alt(m, Epsilon())
Opt!(m::Matcher) = Alt!(m, Epsilon())

# sweet, sweet sugar

# match and discard
~(m::Matcher) = Drop(m)

# refuse
!(m::Matcher) = Not(Lookahead(m))

# match in sequence, result in array

# merge results
+(a::Seq, b::Seq) = Seq(vcat(a.matchers, b.matchers))
+(a::Seq, b::Matcher) = Seq(vcat(a.matchers, b))
+(a::Matcher, b::Seq) = Seq(vcat(a, b.matchers))
+(a::Matcher, b::Matcher) = Seq(a, b)

# separate results
# https://github.com/JuliaLang/julia/issues/11521
(&)(a::And, b::And) = And(vcat(a.matchers, b.matchers))
(&)(a::And, b::Matcher) = And(vcat(a.matchers, b))
(&)(a::Matcher, b::And) = And(vcat(a, b.matchers))
(&)(a::Matcher, b::Matcher) = And(a, b)

# alternates
|(a::Alt, b::Alt) = Alt(vcat(a.matchers, b.matchers))
|(a::Alt, b::Matcher) = Alt(vcat(a.matchers, b))
|(a::Matcher, b::Alt) = Alt(vcat(a, b.matchers))
|(a::Matcher, b::Matcher) = Alt(a, b)

# repeat via [lo:hi] or [n]
lastindex(m::Matcher) = typemax(Int)
size(m::Matcher, n) = endof(m)
axes(m::Matcher, n) = (n == 1) ? Base.OneTo(lastindex(m)) : 1
lastindex(m::Matcher, n) = last(axes(m, n))
getindex(m::Matcher, r::Int, s::Symbol...) = getindex(m, r:r; s...)
function getindex(m::Matcher, r::UnitRange, s::Symbol...)
    greedy = true
    flatten = true
    backtrack = true
    for x in s
        if x == :?
            greedy = false
        elseif x == :&
            flatten = false
        elseif x == :!
            backtrack = false
        else
            error("bad flag to []: $x")
        end
    end
    Repeat(m, r.start, r.stop; greedy=greedy, flatten=flatten, backtrack=backtrack)
end

# the raw Result instance
>=(m::Matcher, f::Applicable) = TransResult(m, f)
# interpolate multiple values (list or tuple)
>(m::Matcher, f::Applicable) = App(m, f)
# a single value
|>(m::Matcher, f::Applicable) = Appl(m, f)


text:
export Freeflyer

# Parameters from:
# asl_free_flyer/free_flyer_node/param/robots/enterprise.yaml
# asl_free_flyer/free_flyer_control/src/waypoint_py_controller/parameters.py
# tribal knowledge

mutable struct Freeflyer{T<:AbstractFloat} <: Robot
  mass_ff_min::T
  mass_ff_max::T
  mass_ff::T
  J_ff::T
  J_ff_inv::T
  J_w::T
  J_w_inv::T
  n_thrusters::Int
  r::T
  hard_limit_vel::T
  hard_limit_accel::T
  hard_limit_ω::T
  hard_limit_α::T
  btCollisionObject

  xb::Vector{T}
  Jcollision
end

function Freeflyer{T}() where T
  n_thrusters = 8
  r = 0.157

	mass_ff_min = 15.36
	mass_ff_max = 18.08
  mass_ff = 0.5*(mass_ff_min+mass_ff_max)

  J_ff = 0.184
  J_ff_inv = inv(J_ff)

  J_w = J_ff/6.43                   # Mark Mote calibration
  J_w_inv = inv(J_w)

  hard_limit_vel = 0.2              # set to match astrobee
  thrust_max = 2*0.185  # max thrust [N] from two thrusters
  hard_limit_accel = thrust_max/mass_ff

  hard_limit_ω = 20*π/180
  wheel_torque_limit = 0.593        # 84 oz-in (https://www.pololu.com/product/2822/specs)
  hard_limit_α = J_w_inv*wheel_torque_limit

  xb = [0.; 0.15; 0.]
  Jcollision = []
  # btCollisionObject = BulletCollision.convex_hull_cylinder(SVector{3}(zeros(3)), SVector{3}(0.,0.,0.5), r)
  btCollisionObjects = BulletCollision.BulletCollisionObjectPtr[]
  push!(btCollisionObjects, BulletCollision.convex_hull_cylinder(SVector{3}(zeros(3)), SVector{3}(0.,0.,1.), r))
  push!(btCollisionObjects, BulletCollision.convex_hull_cylinder(SVector{3}(xb...), SVector{3}(0.,0.,1.), r))

  btCollisionObject = BulletCollision.compound_collision_object(btCollisionObjects)

  # new Freeflyer instance
  return Freeflyer{T}(mass_ff_min,mass_ff_max,mass_ff,J_ff,J_ff_inv,J_w,J_w_inv,n_thrusters,r,
    hard_limit_vel,hard_limit_accel,hard_limit_ω,hard_limit_α,
    btCollisionObject,xb,Jcollision)
end
Freeflyer(::Type{T} = Float64; kwargs...) where {T} = Freeflyer{T}(; kwargs...)


text:
module SimpleFWF

using DataFrames, CSV

export readfwf

function fwfline(line, ranges; stripstr=true)
    data = String[]
    for r in ranges
        ldata = stripstr ? strip(SubString(line, r)) : SubString(line, r)
        push!(data, ldata)
    end
    return data
end

"""
    readfwf(source, ranges; stripstr=true, headers=true) -> DataFrame

Reads a fixed-width file `source` into a `DataFrame` using pre-defined column widths in `ranges`.

# Examples
```julia-repl
julia> df = readfwf("sample.dat", (10:43, 53:60))
```
"""
function readfwf(source, ranges; stripstr=true, headers=true)
    lines = eachline(source)

    # Construct DataFrame
    if headers
        firstline, lines = Base.Iterators.peel(lines)
        colnames = fwfline(firstline, ranges)
    else
        colnames = ["x_$i" for i in eachindex(ranges)]
    end
    columns = Dict(Symbol(c) => String[] for c in colnames)
    df = DataFrame(columns)

    # Read data
    for line in lines
        push!(df, fwfline(line, ranges; stripstr))
    end
    return df
end

end


text:
using Individual.Sampling
using Test

@testset "delay_geom_sample with scalar rate" begin

    @test length(delay_geom_sample(0, 1.0, 1.0)) == 0
    @test length(delay_geom_sample(10, 0.05, 1.0)) == 10

end

@testset "delay_geom_sample with vector rate" begin

    @test_throws ArgumentError delay_geom_sample(1, [1.0, 2.0], 1.0)

    @test length(delay_geom_sample(1, [1.0], 1.0)) == 1
    @test length(delay_geom_sample(10, fill(0.5, 10), 1.0)) == 10

end

@testset "bernoulli_sample with vector target, vector prob" begin

    @test_throws ArgumentError bernoulli_sample([1,2,3], [0.5,0.5,0.5,0.5])
    @test_throws ArgumentError bernoulli_sample([1,2,3], Float64[])
    @test_throws ErrorException bernoulli_sample([1,2,3], [1.0, Inf, 0.5])

    @test length(bernoulli_sample([1,2,3], zeros(3))) == 0
    @test length(bernoulli_sample([1,2,3], ones(3))) == 3
    @test length(bernoulli_sample([1,2,3], [1.0,0.5,1.0])) ∈ [2,3]

end

@testset "bernoulli_sample with vector target, vector rate" begin
    @test_throws ArgumentError bernoulli_sample([1,2,3], [0.5,0.5,0.5,0.5], 1.0)
    @test_throws ArgumentError bernoulli_sample([1,2,3], [0.5, -999.0, 0.5], 1.0)
    @test_throws ArgumentError bernoulli_sample([1,2,3], Float64[], 1.0)

    @test length(bernoulli_sample([1,2,3], zeros(3), 1.0)) == 0
    @test length(bernoulli_sample([1,2,3], fill(Inf, 3), 1.0)) == 3
    @test length(bernoulli_sample([1,2,3], [Inf,1.0,Inf], 1.0)) ∈ [2,3]

end

@testset "bernoulli_sample with vector target, scalar prob" begin

    @test_throws ArgumentError bernoulli_sample([1,2,3], -5.0)

    @test length(bernoulli_sample([1], 0.0)) == 0
    @test length(bernoulli_sample([1], 1.0)) == 1

    @test length(bernoulli_sample([1,2,3], 0.0)) == 0
    @test length(bernoulli_sample([1,2,3], 1.0)) == 3

end

@testset "bernoulli_sample with vector target, scalar rate" begin

    @test_throws ArgumentError bernoulli_sample([1,2,3], -999.0, 1.0)

    @test length(bernoulli_sample([1,2,3], 0.0, 1.0)) == 0
    @test length(bernoulli_sample([1,2,3], Inf, 1.0)) == 3

    @test length(bernoulli_sample([1], 0.0, 1.0)) == 0
    @test length(bernoulli_sample([1], Inf, 1.0)) == 1

    @test bernoulli_sample([1,2,3], Inf, 1.0) == [1,2,3]

end

@testset "bernoulli_sample with scalar target, scalar prob" begin

    @test_throws ArgumentError bernoulli_sample(1, -5.0)

    @test length(bernoulli_sample(1, 0.0)) == 0
    @test length(bernoulli_sample(1, 1.0)) == 1
    @test bernoulli_sample(1, 1.0) == [1]

end

@testset "bernoulli_sample with scalar target, scalar rate" begin

    @test_throws ArgumentError bernoulli_sample(1, -999.0, 1.0)

    @test length(bernoulli_sample(1, 0.0, 1.0)) == 0
    @test length(bernoulli_sample(1, Inf, 1.0)) == 1

    @test bernoulli_sample(1, Inf, 1.0) == [1]

end


text:
@testset "layer" begin
    T = Float32
    in_channel = 3
    in_channel_edge = 1
    out_channel = 5

    N = 4
    E = 4
    adj = T[0. 1. 0. 1.;
            1. 0. 1. 0.;
            0. 1. 0. 1.;
            1. 0. 1. 0.]
    fg = FeaturedGraph(adj)

    @testset "GCNConv" begin
        X = rand(T, in_channel, N)
        Xt = transpose(rand(T, N, in_channel))
        @testset "layer with graph" begin
            gc = GCNConv(fg, in_channel=>out_channel)
            @test size(gc.weight) == (out_channel, in_channel)
            @test size(gc.bias) == (out_channel,)
            @test GraphSignals.adjacency_matrix(gc.fg) == adj

            Y = gc(X)
            @test size(Y) == (out_channel, N)

            # Test with transposed features
            Y = gc(Xt)
            @test size(Y) == (out_channel, N)

            g = Zygote.gradient(() -> sum(gc(X)), Flux.params(gc))
            @test length(g.grads) == 2
        end

        @testset "layer without graph" begin
            gc = GCNConv(in_channel=>out_channel)
            @test size(gc.weight) == (out_channel, in_channel)
            @test size(gc.bias) == (out_channel,)
            @test !has_graph(gc.fg)

            fg = FeaturedGraph(adj, nf=X)
            fg_ = gc(fg)
            @test size(node_feature(fg_)) == (out_channel, N)
            @test_throws ArgumentError gc(X)

            # Test with transposed features
            fgt = FeaturedGraph(adj, nf=Xt)
            fgt_ = gc(fgt)
            @test size(node_feature(fgt_)) == (out_channel, N)

            g = Zygote.gradient(() -> sum(node_feature(gc(fg))), Flux.params(gc))
            @test length(g.grads) == 4
        end

        @testset "bias=false" begin
            @test length(Flux.params(GCNConv(2=>3))) == 2
            @test length(Flux.params(GCNConv(2=>3, bias=false))) == 1
        end
    end

    @testset "ChebConv" begin
        k = 6
        X = rand(T, in_channel, N)
        Xt = transpose(rand(T, N, in_channel))
        @testset "layer with graph" begin
            cc = ChebConv(fg, in_channel=>out_channel, k)
            @test size(cc.weight) == (out_channel, in_channel, k)
            @test size(cc.bias) == (out_channel,)
            @test GraphSignals.adjacency_matrix(cc.fg) == adj
            @test cc.k == k

            Y = cc(X)
            @test size(Y) == (out_channel, N)

            # Test with transposed features
            Y = cc(Xt)
            @test size(Y) == (out_channel, N)

            g = Zygote.gradient(() -> sum(cc(X)), Flux.params(cc))
            @test length(g.grads) == 2
        end

        @testset "layer without graph" begin
            cc = ChebConv(in_channel=>out_channel, k)
            @test size(cc.weight) == (out_channel, in_channel, k)
            @test size(cc.bias) == (out_channel,)
            @test !has_graph(cc.fg)
            @test cc.k == k

            fg = FeaturedGraph(adj, nf=X)
            fg_ = cc(fg)
            @test size(node_feature(fg_)) == (out_channel, N)
            @test_throws ArgumentError cc(X)

            # Test with transposed features
            fgt = FeaturedGraph(adj, nf=Xt)
            fgt_ = cc(fgt)
            @test size(node_feature(fgt_)) == (out_channel, N)

            g = Zygote.gradient(() -> sum(node_feature(cc(fg))), Flux.params(cc))
            @test length(g.grads) == 4
        end

        @testset "bias=false" begin
            @test length(Flux.params(ChebConv(2=>3, 3))) == 2
            @test length(Flux.params(ChebConv(2=>3, 3, bias=false))) == 1
        end
    end

    @testset "GraphConv" begin
        X = rand(T, in_channel, N)
        Xt = transpose(rand(T, N, in_channel))
        @testset "layer with graph" begin
            gc = GraphConv(fg, in_channel=>out_channel)
            @test adjacency_list(gc.fg) == [[2,4], [1,3], [2,4], [1,3]]
            @test size(gc.weight1) == (out_channel, in_channel)
            @test size(gc.weight2) == (out_channel, in_channel)
            @test size(gc.bias) == (out_channel,)

            Y = gc(X)
            @test size(Y) == (out_channel, N)

            # Test with transposed features
            Y = gc(Xt)
            @test size(Y) == (out_channel, N)

            g = Zygote.gradient(() -> sum(gc(X)), Flux.params(gc))
            @test length(g.grads) == 3
        end

        @testset "layer without graph" begin
            gc = GraphConv(in_channel=>out_channel)
            @test size(gc.weight1) == (out_channel, in_channel)
            @test size(gc.weight2) == (out_channel, in_channel)
            @test size(gc.bias) == (out_channel,)

            fg = FeaturedGraph(adj, nf=X)
            fg_ = gc(fg)
            @test size(node_feature(fg_)) == (out_channel, N)
            @test_throws ArgumentError gc(X)

            # Test with transposed features
            fgt = FeaturedGraph(adj, nf=Xt)
            fgt_ = gc(fgt)
            @test size(node_feature(fgt_)) == (out_channel, N)

            g = Zygote.gradient(() -> sum(node_feature(gc(fg))), Flux.params(gc))
            @test length(g.grads) == 5
        end

        @testset "bias=false" begin
            @test length(Flux.params(GraphConv(2=>3))) == 3
            @test length(Flux.params(GraphConv(2=>3, bias=false))) == 2
        end
    end

    @testset "GATConv" begin
        adj1 = [1 1 0 1;
                1 1 1 0;
                0 1 1 1;
                1 0 1 1]
        fg1 = FeaturedGraph(adj1)

        # isolated_vertex
        adj2 = [1 0 0 1;
                0 1 0 0;
                0 0 1 1;
                1 0 1 1]
        fg2 = FeaturedGraph(adj2)

        X = rand(T, in_channel, N)
        Xt = transpose(rand(T, N, in_channel))

        @testset "layer with graph" begin
            for heads = [1, 2], concat = [true, false], adj_gat in [adj1, adj2]
                fg_gat = FeaturedGraph(adj_gat)
                gat = GATConv(fg_gat, in_channel=>out_channel, heads=heads, concat=concat)

                @test size(gat.weight) == (out_channel * heads, in_channel)
                @test size(gat.bias) == (out_channel * heads,)
                @test size(gat.a) == (2*out_channel, heads)

                Y = gat(X)
                @test size(Y) == (concat ? (out_channel*heads, N) : (out_channel, N))

                # Test with transposed features
                Y = gat(Xt)
                @test size(Y) == (concat ? (out_channel*heads, N) : (out_channel, N))

                g = Zygote.gradient(() -> sum(gat(X)), Flux.params(gat))
                @test length(g.grads) == 3
            end
        end

        @testset "layer without graph" begin
            for heads = [1, 2], concat = [true, false], adj_gat in [adj1, adj2]
                fg_gat = FeaturedGraph(adj_gat, nf=X)
                gat = GATConv(in_channel=>out_channel, heads=heads, concat=concat)
                @test size(gat.weight) == (out_channel * heads, in_channel)
                @test size(gat.bias) == (out_channel * heads,)
                @test size(gat.a) == (2*out_channel, heads)

                fg_ = gat(fg_gat)
                Y = node_feature(fg_)
                @test size(Y) == (concat ? (out_channel*heads, N) : (out_channel, N))
                @test_throws ArgumentError gat(X)

                # Test with transposed features
                fgt = FeaturedGraph(adj_gat, nf=Xt)
                fgt_ = gat(fgt)
                @test size(node_feature(fgt_)) == (concat ? (out_channel*heads, N) : (out_channel, N))

                g = Zygote.gradient(() -> sum(node_feature(gat(fg_gat))), Flux.params(gat))
                @test length(g.grads) == 5
            end
        end

        @testset "bias=false" begin
            @test length(Flux.params(GATConv(2=>3))) == 3
            @test length(Flux.params(GATConv(2=>3, bias=false))) == 2
        end
    end

    @testset "GatedGraphConv" begin
        num_layers = 3
        X = rand(T, in_channel, N)
        Xt = transpose(rand(T, N, in_channel))
        @testset "layer with graph" begin
            ggc = GatedGraphConv(fg, out_channel, num_layers)
            @test adjacency_list(ggc.fg) == [[2,4], [1,3], [2,4], [1,3]]
            @test size(ggc.weight) == (out_channel, out_channel, num_layers)

            Y = ggc(X)
            @test size(Y) == (out_channel, N)

            # Test with transposed features
            Y = ggc(Xt)
            @test size(Y) == (out_channel, N)

            g = Zygote.gradient(() -> sum(ggc(X)), Flux.params(ggc))
            @test length(g.grads) == 6
        end

        @testset "layer without graph" begin
            ggc = GatedGraphConv(out_channel, num_layers)
            @test size(ggc.weight) == (out_channel, out_channel, num_layers)

            fg = FeaturedGraph(adj, nf=X)
            fg_ = ggc(fg)
            @test size(node_feature(fg_)) == (out_channel, N)
            @test_throws ArgumentError ggc(X)

            # Test with transposed features
            fgt = FeaturedGraph(adj, nf=Xt)
            fgt_ = ggc(fgt)
            @test size(node_feature(fgt_)) == (out_channel, N)

            g = Zygote.gradient(() -> sum(node_feature(ggc(fg))), Flux.params(ggc))
            @test length(g.grads) == 8
        end
    end

    @testset "EdgeConv" begin
        X = rand(T, in_channel, N)
        Xt = transpose(rand(T, N, in_channel))
        @testset "layer with graph" begin
            ec = EdgeConv(fg, Dense(2*in_channel, out_channel))
            @test adjacency_list(ec.fg) == [[2,4], [1,3], [2,4], [1,3]]

            Y = ec(X)
            @test size(Y) == (out_channel, N)

            # Test with transposed features
            Y = ec(Xt)
            @test size(Y) == (out_channel, N)

            g = Zygote.gradient(() -> sum(ec(X)), Flux.params(ec))
            @test length(g.grads) == 2
        end

        @testset "layer without graph" begin
            ec = EdgeConv(Dense(2*in_channel, out_channel))

            fg = FeaturedGraph(adj, nf=X)
            fg_ = ec(fg)
            @test size(node_feature(fg_)) == (out_channel, N)
            @test_throws ArgumentError ec(X)

            # Test with transposed features
            fgt = FeaturedGraph(adj, nf=Xt)
            fgt_ = ec(fgt)
            @test size(node_feature(fgt_)) == (out_channel, N)

            g = Zygote.gradient(() -> sum(node_feature(ec(fg))), Flux.params(ec))
            @test length(g.grads) == 4
        end
    end

    @testset "GINConv" begin
        X = rand(Float32, in_channel, N)
        Xt = transpose(rand(Float32, N, in_channel))
        nn = Flux.Chain(Dense(in_channel, out_channel))
        eps = 0.001

        @testset "layer with graph" begin
            gc = GINConv(FeaturedGraph(adj), nn, eps)
            @test size(gc.nn.layers[1].weight) == (out_channel, in_channel)
            @test size(gc.nn.layers[1].bias) == (out_channel, )
            @test GraphSignals.adjacency_matrix(gc.fg) == adj

            fg = FeaturedGraph(adj, nf=X)
            Y = gc(fg)
            @test size(node_feature(Y)) == (out_channel, N)

            # Test with transposed features
            Y = gc(FeaturedGraph(adj, nf=Xt))
            @test size(node_feature(Y)) == (out_channel, N)

            g = Zygote.gradient(() -> sum(node_feature(gc(fg))), Flux.params(gc))
            @test length(g.grads) == 8
        end
    end

    @testset "CGConv" begin
        fg = FeaturedGraph(adj)
        X = rand(Float32, in_channel, N)
        E = rand(Float32, in_channel_edge, ne(fg))
        Xt = transpose(rand(Float32, N, in_channel))
        @testset "layer with graph" begin
            cgc = CGConv(FeaturedGraph(adj),
                         (in_channel, in_channel_edge))
            @test size(cgc.Wf) == (in_channel, 2 * in_channel + in_channel_edge)
            @test size(cgc.Ws) == (in_channel, 2 * in_channel + in_channel_edge)
            @test size(cgc.bf) == (in_channel,)
            @test size(cgc.bs) == (in_channel,)

            Y = cgc(X, E)
            @test size(Y) == (in_channel, N)

            Yg = cgc(FeaturedGraph(adj, nf=X, ef=E))
            @test size(node_feature(Yg)) == (in_channel, N)
            @test edge_feature(Yg) == E
        end
    end
end


text:
# Packages required for testing
using Test
using Random
using Logging

# Package Under Test
using GridWorldPathFollowing
using Vec
using LinearAlgebra

# Set logging level
global_logger(SimpleLogger(stderr, Logging.Debug))

# Fix randomness during tests
Random.seed!(0)

# Check equality of two arrays
@inline function array_isapprox(x::AbstractArray{F},
                  y::AbstractArray{F};
                  rtol::F=sqrt(eps(F)),
                  atol::F=zero(F)) where {F<:AbstractFloat}

    # Easy check on matching size
    if length(x) != length(y)
        return false
    end

    for (a,b) in zip(x,y)
        if !isapprox(a,b, rtol=rtol, atol=atol)
            return false
        end
    end
    return true
end

# Check if array equals a single value
@inline function array_isapprox(x::AbstractArray{F},
                  y::F;
                  rtol::F=sqrt(eps(F)),
                  atol::F=zero(F)) where {F<:AbstractFloat}

    for a in x
        if !isapprox(a,y, rtol=rtol, atol=atol)
            return false
        end
    end
    return true
end

# Define package tests
@time @testset "GridWorldPathFollowing Package Tests" begin
    testdir = joinpath(dirname(@__DIR__), "test")
    @time @testset "GridWorldPathFollowing.Utils" begin
        include(joinpath(testdir, "test_utils.jl"))
    end
    @time @testset "GridWorldPathFollowing.GridPaths" begin
        include(joinpath(testdir, "test_grid_paths.jl"))
    end
    @time @testset "GridWorldPathFollowing.Trajectories" begin
        include(joinpath(testdir, "test_trajectories.jl"))
    end
    @time @testset "GridWorldPathFollowing.RobotModels" begin
        include(joinpath(testdir, "test_robot_models.jl"))
    end
    # @time @testset "GridWorldPathFollowing.Deploy" begin
    #     include(joinpath(testdir, "test_deploy.jl"))
    # end
end


text:
#=
P34 (**) Calculate Euler's totient function phi(m).
Euler's so-called totient function phi(m) is defined as the number of positive integers r (1 <= r < m) that are coprime to m.
Example: m = 10: r = 1,3,7,9; thus phi(m) = 4. Note the special case: phi(1) = 1.

* (totient-phi 10)
4

Find out what the value of phi(m) is if m is a prime number. Euler's totient function plays an important role in one of the most widely used public key cryptography methods (RSA). In this exercise you should use the most primitive method to calculate this function (there are smarter ways that we shall discuss later).
=#

# From P32 but more terse
function gcd(a, b)
    b == 0 ? a : gcd(b, mod(a, b))
end

# From P33
function coprime(a, b)
    gcd(a, b) == 1
end

function totient_phi(m)
    # Starting at 2, not 1 (special case)
    coprime_to_m = filter(r -> coprime(r, m), 2:m)
    length(coprime_to_m) + 1
end

@assert totient_phi(10) == 4

println("Tests passed: JL-34.jl")


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMBase.jl/blob/master/LICENSE

const Node = Vector{Float64}

abstract type AbstractPoint end

mutable struct Point{P<:AbstractPoint}
    id :: Int
    weight :: Float64
    coords :: Tuple{Vararg{Float64}}
    fields :: Dict{String, AbstractField}
    properties :: P
end

function setindex!(point::Point, val::Pair{Float64, T}, field_name) where T
    point.fields[field_name] = field(val)
end

function getindex(point::Point, field_name)
    return point.fields[field_name]
end

function getindex(point::Point, idx::Int)
    return point.coords[idx]
end

function haskey(point::Point, field_name)
    return haskey(point.fields, field_name)
end

function (point::Point)(field_name, time)
    interpolate(point.fields[field_name], time)
end

function start(point::Point)
    return start(point.coords)
end

function done(point::Point, i)
    return done(point.coords, i)
end

function next(point::Point, i)
    return next(point.coords, i)
end

function update!(point::Point, field_name, val::Pair{Float64, T}) where T
    if haskey(point, field_name)
        update!(point[field_name], val)
    else
        point[field_name] = val
    end
end

#= TODO: in future
type Node <: AbstractPoint
end

type MaterialPoint <: AbstractPoint
end
=#

mutable struct IntegrationPoint <: AbstractPoint
end

const IP = Point{IntegrationPoint}

function IP(id, weight, coords::Tuple)
    return IP(id, weight, coords, Dict(), IntegrationPoint())
end

function IP(id, weight, coords::Vector)
    warn("Consider giving coords as Tuple.")
    return IP(id, weight, tuple(coords...), Dict(), IntegrationPoint())
end


text:
function check_status(status)
  if status != 0
    error("libbbn-ni-r failed with error number $status")
  end
end

function open_load_run()
    check_status(ccall((:open_load_run, "libbbn-ni-r"), Cint, ()))
end

function set_numSamples(numSamples)
  check_status(ccall((:set_numSamples, "libbbn-ni-r"), Cint, (Cuint,), numSamples))
end

function get_numSamples()
  numSamples = Ref{Cuint}(0)
  check_status(ccall((:get_numSamples, "libbbn-ni-r"), Cint, (Ref{Cuint},), numSamples))
  return numSamples[]
end

function set_sampleInterval(sampleInterval)
  check_status(ccall((:set_sampleInterval, "libbbn-ni-r"), Cint, (Cdouble,), sampleInterval))
end

function get_sampleInterval()
  sampleInterval = Ref{Cdouble}(0)
  check_status(ccall((:get_sampleInterval, "libbbn-ni-r"), Cint, (Ref{Cdouble},), sampleInterval))
  return sampleInterval[]
end

function get_sampleCt()
  sampleCt = Ref{Cint}(0)
  check_status(ccall((:get_sampleCt, "libbbn-ni-r"), Cint, (Ref{Cint},), sampleCt))
  return sampleCt[]
end

function enable_acquisition()
  check_status(ccall((:enable_acquisition, "libbbn-ni-r"), Cint, ()))
end

function reset_fifo()
  check_status(ccall((:reset_fifo, "libbbn-ni-r"), Cint, ()))
end

function disable_acquisition()
  check_status(ccall((:disable_acquisition, "libbbn-ni-r"), Cint, ()))
end

function trigger()
  check_status(ccall((:trigger, "libbbn-ni-r"), Cint, ()))
end

function trigger_out()
  check_status(ccall((:trigger_out, "libbbn-ni-r"), Cint, ()))
end

function transfer_waveform(numPts, timeOut)
  data = Array(Float64, numPts)
  check_status(ccall((:transfer_waveform, "libbbn-ni-r"), Cint,
                      (Cuint, Ptr{Cdouble}, Cdouble),
                      numPts, data, timeOut))
  return data
end

function get_analogOut()
  val = Ref{Cdouble}(0)
  check_status(ccall((:get_analogOut, "libbbn-ni-r"), Cint, (Ref{Cdouble},), val))
  return val[]
end

function get_analogIn()
  val = Ref{Cdouble}(0)
  check_status(ccall((:get_analogIn, "libbbn-ni-r"), Cint, (Ref{Cdouble},), val))
  return val[]
end

max_stars_repo_name: BBN-Q/NI-RSeries

max_issues_repo_name: BBN-Q/NI-RSeries

max_forks_repo_name: BBN-Q/NI-RSeries


text:
abstract type CompoundingFrequency end
Base.Broadcast.broadcastable(x::T) where {T<:CompoundingFrequency} = Ref(x)

"""
    Continuous()

A type representing continuous interest compounding frequency.

# Examples

```julia-repl
julia> Rate(0.01,Continuous())
Rate(0.01, Continuous())
```

See also: [`Periodic`](@ref)
"""
struct Continuous <: CompoundingFrequency end

"""

```julia-repl
julia> Continuous(0.01)
Rate(0.01, Continuous())
```

See also: [`Periodic`](@ref)
"""
Continuous(rate) = Rate(rate, Continuous())

"""
    Periodic(frequency)

A type representing periodic interest compounding with the given frequency

# Examples

Creating a semi-annual bond equivalent yield:

```julia-repl
julia> Rate(0.01,Periodic(2))
Rate(0.01, Periodic(2))
```

See also: [`Continuous`](@ref)
"""
struct Periodic <: CompoundingFrequency
    frequency::Int
end

"""
    Periodic(rate,frequency)

A convenience constructor for Rate(rate,Periodic(frequency)).

# Examples

Creating a semi-annual bond equivalent yield:

```julia-repl
julia> Periodic(0.01,2)
Rate(0.01, Periodic(2))
```

See also: [`Continuous`](@ref)
"""
Periodic(x, frequency) = Rate(x, Periodic(frequency))

struct Rate{N<:Real,T<:CompoundingFrequency} <: AbstractYield
    value::N
    compounding::T
end

# make rate a broadcastable type
Base.Broadcast.broadcastable(ic::T) where {T<:Rate} = Ref(ic)

"""
    Rate(rate[,frequency=1])
    Rate(rate,frequency::CompoundingFrequency)

Rate is a type that encapsulates an interest `rate` along with its compounding `frequency`.

Periodic rates can be constructed via `Rate(rate,frequency)` or `Rate(rate,Periodic(frequency))`. If not given a second argument, `Rate(rate)` is equivalent to `Rate(rate,Periodic(1))`.

Continuous rates can be constructed via `Rate(rate, Inf)` or `Rate(rate,Continuous())`.

# Examples

```julia-repl
julia> Rate(0.01,Continuous())
Rate(0.01, Continuous())

julia> Rate(0.01,Periodic(2))
Rate(0.01, Periodic(2))

julia> Rate(0.01)
Rate(0.01, Periodic(1))

julia> Rate(0.01,2)
Rate(0.01, Periodic(2))

julia> Rate(0.01,Periodic(4))
Rate(0.01, Periodic(4))

julia> Rate(0.01,Inf)
Rate(0.01, Continuous())

julia> Rate(0.01,Continuous())
Rate(0.01, Continuous())
```
"""
Rate(rate) = Rate(rate, Periodic(1))
Rate(x, frequency::T) where {T<:Real} = isinf(frequency) ? Rate(x, Continuous()) : Rate(x, Periodic(frequency))

"""
    convert(T::CompoundingFrequency,r::Rate)

Returns a `Rate` with an equivalent discount but represented with a different compounding frequency.

# Examples

```julia-repl
julia> r = Rate(Periodic(12),0.01)
Rate(0.01, Periodic(12))

julia> convert(Periodic(1),r)
Rate(0.010045960887181016, Periodic(1))

julia> convert(Continuous(),r)
Rate(0.009995835646701251, Continuous())
```
"""
function Base.convert(T::CompoundingFrequency, r::Rate{<:Real,<:CompoundingFrequency})
    convert(T, r, r.compounding)
end
function Base.convert(to::Continuous, r, from::Continuous)
    return r
end

function Base.convert(to::Periodic, r, from::Continuous)
    return Rate(to.frequency * (exp(r.value / to.frequency) - 1), to)
end

function Base.convert(to::Continuous, r, from::Periodic)
    return Rate(from.frequency * log(1 + r.value / from.frequency), to)
end

function Base.convert(to::Periodic, r, from::Periodic)
    c = convert(Continuous(), r, from)
    return convert(to, c, Continuous())
end

function rate(r::Rate{<:Real,<:CompoundingFrequency})
    r.value
end

function Base.isapprox(a::Rate{N,T}, b::Rate{N,T}; atol::Real = 0, rtol::Real = atol > 0 ? 0 : √eps()) where {T<:Periodic,N<:Real}
    c = convert(a.compounding,b)
    return (a.compounding.frequency == c.compounding.frequency) && isapprox(rate(a), rate(c); atol, rtol)
end

function Base.isapprox(a::Rate{N,T}, b::Rate{N,T}; atol::Real = 0, rtol::Real = atol > 0 ? 0 : √eps()) where {T<:Continuous,N<:Real}
    return isapprox(rate(a), rate(b); atol, rtol)
end

# the fallback for rates not of the same type
function Base.isapprox(a::T, b::N; atol::Real = 0, rtol::Real = atol > 0 ? 0 : √eps()) where {T<:Rate,N<:Rate}
    return isapprox(convert(b.compounding, a), b; atol, rtol)
end

"""
    discount(rate, t)
    discount(rate, from, to)

Discount `rate` for a time `t` or for an interval `(from, to)`. If `rate` is not a `Rate`, it will be assumed to be a `Periodic` rate compounded once per period, i.e. `Periodic(rate,1)`.

# Examples

```julia-repl
julia> discount(0.03, 10)
0.7440939148967249

julia> discount(Periodic(0.03, 2), 10)
0.7424704182237725

julia> discount(Continuous(0.03), 10)
0.7408182206817179

julia> discount(0.03, 5, 10)
0.8626087843841639
```
"""

discount(rate::T, t) where {T<:Real} = discount(Rate(rate), t)
discount(rate::Rate{<:Real, <:Continuous}, t) = exp(-rate.value * t)
discount(rate::Rate{<:Real, <:Periodic}, t) = (1 + rate.value / rate.compounding.frequency)^(-rate.compounding.frequency * t)
discount(rate, from, to) = discount(rate, to - from)

"""
    accumulation(rate, t)
    accumulation(rate, from, to)

Accumulate `rate` for a time `t` or for an interval `(from, to)`. If `rate` is not a `Rate`, it will be assumed to be a `Periodic` rate compounded once per period, i.e. `Periodic(rate,1)`.

    # Examples

```julia-repl
julia> accumulation(0.03, 10)
1.3439163793441222

julia> accumulation(Periodic(0.03, 2), 10)
1.3468550065500535

julia> accumulation(Continuous(0.03), 10)
1.3498588075760032

julia> accumulation(0.03, 5, 10)
1.1592740743
```
"""

accumulation(rate, t) = accumulation(Rate(rate), t)
accumulation(rate::Rate{<:Real, <:Continuous}, t) = exp(rate.value * t)
accumulation(rate::Rate{<:Real, <:Periodic}, t) = (1 + rate.value / rate.compounding.frequency)^(rate.compounding.frequency * t)
accumulation(rate, from, to) = accumulation(rate, to - from)

Base.zero(rate::T,t) where {T<:Rate} = rate
forward(rate::T,to) where {T<:Rate} = rate
forward(rate::T,from,to) where {T<:Rate} = rate

"""
    +(Yields.Rate, T<:Real)
    +(T<:Real, Yields.Rate)
    +(Yields.Rate,Yields.Rate)

The addition of a rate with a number will inherit the type of the `Rate`, or the first argument's type if both are `Rate`s.

# Examples

```julia-repl
julia> Yields.Periodic(0.01,2) + Yields.Periodic(0.04,2)
Yields.Rate{Float64, Yields.Periodic}(0.05000000000000004, Yields.Periodic(2))

julia> Yields.Periodic(0.04,2) + 0.01
Yields.Rate{Float64, Yields.Periodic}(0.05, Yields.Periodic(2))
```
"""
function Base.:+(a::Rate{N,T}, b::Real) where {N<:Real,T<:Continuous}
    return Continuous(a.value + b)
end
function Base.:+(a::Real, b::Rate{N,T}) where {N<:Real,T<:Continuous}
    return Continuous(b.value + a)
end

function Base.:+(a::Rate{N,T}, b::Real) where {N<:Real, T<:Periodic}
    return Periodic(a.value + b,a.compounding.frequency)
end
function Base.:+(a::Real, b::Rate{N,T}) where {N<:Real, T<:Periodic}
    return Periodic(b.value + a,b.compounding.frequency)
end

function Base.:+(a::Rate{N,T},b::Rate{N,T}) where {N<:Real,T<:CompoundingFrequency}
    a_rate = rate(a)
    b_rate = rate(convert(a.compounding,b))
    r = Rate(a_rate + b_rate,a.compounding)
    return r
end

"""
    -(Yields.Rate, T<:Real)
    -(T<:Real, Yields.Rate)
    -(Yields.Rate, Yields.Rate)

The addition of a rate with a number will inherit the type of the `Rate`, or the first argument's type if both are `Rate`s.

# Examples

```julia-repl
julia> Yields.Periodic(0.04,2) - Yields.Periodic(0.01,2)
Yields.Rate{Float64, Yields.Periodic}(0.030000000000000214, Yields.Periodic(2))

julia> Yields.Periodic(0.04,2) - 0.01
Yields.Rate{Float64, Yields.Periodic}(0.03, Yields.Periodic(2))

```
"""
function Base.:-(a::Rate{N,T}, b::Real) where {N<:Real,T<:Continuous}
    return Continuous(a.value - b)
end
function Base.:-(a::Real, b::Rate{N,T}) where {N<:Real,T<:Continuous}
    return Continuous( a - b.value)
end

function Base.:-(a::Rate{N,T}, b::Real) where {N<:Real, T<:Periodic}
    return Periodic(a.value - b, a.compounding.frequency)
end
function Base.:-(a::Real, b::Rate{N,T}) where {N<:Real, T<:Periodic}
    return Periodic(a - b.value, b.compounding.frequency)
end
function Base.:-(a::Rate{N,T},b::Rate{N,T}) where {N<:Real,T<:CompoundingFrequency}
    a_rate = rate(a)
    b_rate = rate(convert(a.compounding,b))
    r = Rate(a_rate - b_rate,a.compounding)
    return r
end

"""
    *(Yields.Rate, T<:Real)
    *(T<:Real, Yields.Rate)

The multiplication of a Rate with a scalar will inherit the type of the `Rate`, or the first argument's type if both are `Rate`s.
"""
function Base.:*(a::Rate{N,T}, b::Real) where {N<:Real,T<:Continuous}
    return Continuous(a.value * b)
end
function Base.:*(a::Real, b::Rate{N,T}) where {N<:Real,T<:Continuous}
    return Continuous( a * b.value)
end

function Base.:*(a::Rate{N,T}, b::Real) where {N<:Real, T<:Periodic}
    return Periodic(a.value * b, a.compounding.frequency)
end
function Base.:*(a::Real, b::Rate{N,T}) where {N<:Real, T<:Periodic}
    return Periodic(a * b.value, b.compounding.frequency)
end

"""
    /(Yields.Rate, T<:Real)
    /(T<:Real, Yields.Rate)

The division of a Rate with a scalar will inherit the type of the `Rate`, or the first argument's type if both are `Rate`s.
"""
function Base.:/(a::Rate{N,T}, b::Real) where {N<:Real,T<:Continuous}
    return Continuous(a.value / b)
end

# unclear if dividing a scalar by a rate should be allowed
# function Base.:/(a::Real, b::Rate{N,T}) where {N<:Real,T<:Continuous}
#     return Continuous( a / b.value)
# end

function Base.:/(a::Rate{N,T}, b::Real) where {N<:Real, T<:Periodic}
    return Periodic(a.value / b, a.compounding.frequency)
end

# unclear if dividing a scalar by a rate should be allowed
# function Base.:/(a::Real, b::Rate{N,T}) where {N<:Real, T<:Periodic}
#     return Periodic(a / b.value, b.compounding.frequency)
# end

"""
    <(Rate,Rate)

Convert the second argument to the periodicity of the first and compare the scalar rate values to determine if the first argument has a lower force of interest than the second.

# Examples

```julia-repl
julia> Yields.Periodic(0.03,100) < Yields.Continuous(0.03)
true
```
"""
function Base.:<(a::Rate{N,T},b::Rate{N,U}) where {N<:Real,T<:CompoundingFrequency,U<:CompoundingFrequency}
    bc = convert(a.compounding,b)
    return rate(a) < rate(bc)
end

"""
    >(Rate,Rate)

Convert the second argument to the periodicity of the first and compare the scalar rate values to determine if the first argument has a greater force of interest than the second.

# Examples

```julia-repl
julia> Yields.Periodic(0.03,100) > Yields.Continuous(0.03)
false
```
"""
function Base.:>(a::Rate{N,T},b::Rate{N,U}) where {N<:Real,T<:CompoundingFrequency,U<:CompoundingFrequency}
    bc = convert(a.compounding,b)
    return rate(a) > rate(bc)
end


text:
using Whirl
using Plots

gr()

Δt = 1e-4
n_steps = 10000
n_skips = 50

builder = @show sod_shock_tube_builder(number_of_particles = 1000)
stepper = build(builder)

anim = @animate for _ = 1:n_steps
    sphplot(get_positions(stepper), get_mass_densities(stepper))
    step!(stepper, Δt)
end every n_skips
gif(anim, "sod_shock_tube.gif", fps = 15)

max_stars_repo_name: lars-frogner/Whirl

max_issues_repo_name: lars-frogner/Whirl

max_forks_repo_name: lars-frogner/Whirl


text:
push!(LOAD_PATH, string(pwd(),"/src"))

using Revise
using PyPlot
using ModStokesFastDirect
import AnalyticDomains
import PUX
include("../src/plottools.jl")

mss = ModifiedStokesSolver
unitcircle = AnalyticDomains.starfish(n_arms=3, amplitude=0.00)

function estimate_fourier_resolution(U)
    Uhat = fftshift(fft(U))
    Ures = max( norm(Uhat[1:2,:], Inf),
                norm(Uhat[end-1:end,:], Inf),
                norm(Uhat[:,1:2], Inf),
                norm(Uhat[:,end-1:end], Inf) ) / maximum(abs.(Uhat))
    return Ures
end

curve = AnalyticDomains.starfish()
Ngrid = 800
numpanels = 400
pux_R = 0.15

# Dev
# curve = unitcircle
# Ngrid = 100
# numpanels = 100

panelorder = 16
method = :dense_gmres
pux_ep = 2.0

plotext = false
ploterror=true

# Define ref problem (linear oscillatory)
k1 = 1.0
k2 = 1.0
ksq = k1^2 + k2^2
phi_k(x, y) = exp(1im*(k1*x + k2*y))
fhat1 = 1.0
fhat2 = -2.0
kdotfhat = k1*fhat1 + k2*fhat2
# rhs
ffunc1(x, y) = real(fhat1*phi_k(x,y))
ffunc2(x, y) = real(fhat2*phi_k(x,y))
# solution
ufunc1(x, y, alpha) = real( (ksq*fhat1 - k1*kdotfhat)/(ksq*(ksq+alpha^2))*phi_k(x,y) )
ufunc2(x, y, alpha) = real( (ksq*fhat2 - k2*kdotfhat)/(ksq*(ksq+alpha^2))*phi_k(x,y) )
# derivatives of solution
ufunc1x(x, y, alpha) = real( (ksq*fhat1 - k1*kdotfhat)/(ksq*(ksq+alpha^2))*phi_k(x,y)*1im*k1 )
ufunc1y(x, y, alpha) = real( (ksq*fhat1 - k1*kdotfhat)/(ksq*(ksq+alpha^2))*phi_k(x,y)*1im*k2 )
ufunc2x(x, y, alpha) = real( (ksq*fhat2 - k2*kdotfhat)/(ksq*(ksq+alpha^2))*phi_k(x,y)*1im*k1 )
ufunc2y(x, y, alpha) = real( (ksq*fhat2 - k2*kdotfhat)/(ksq*(ksq+alpha^2))*phi_k(x,y)*1im*k2 )

function solve(pux_R, alpha)
    info(">>>>>>>>>>>>>>>>>>>>>> alpha=$alpha")
    @show Ngrid

    # Discretize
    println("* Discretizing")
    dcurve = mss.discretize(curve, numpanels, panelorder, equal_arclength=true)

    if method==:fastdirect
        # Fast direct prep
        println("* Fast direct build")
        @time Ainv = compressed_system_matrix_inverse(dcurve, alpha, 64, fastdirect_TOL)
    elseif method==:fmm
        # GMRES prep
        println("* GMRES FMM prep")
        @time flhs = mss.system_matvec(dcurve, alpha)
    elseif method==:dense_lu
        println("* Dense matrix assembly")
        @time MAT = mss.system_matrix(dcurve, alpha)
        println("* Dense matrix LU")
        @time MATL, MATU, MATp = lu(MAT)
    elseif method==:dense_gmres
        println("* Dense matrix assembly")
        @time MAT = mss.system_matrix(dcurve, alpha)
    else
        error("unknown method")
    end

    println("* PUX precompute")
    @time PUXStor, X, Y, Lgrid, Ngrid_new, interior = PUX.pux_precompute(curve, dcurve, Ngrid, pux_ep, R=pux_R)
    xbdry = dcurve.points[1,:]
    ybdry = dcurve.points[2,:]
    interior_points = [X[interior]' ; Y[interior]']
    numeval = size(interior_points, 2)
    println("* FMM tree precompute")
    @time fmm_precomp = mss.doublelayer_precomp(dcurve, interior_points, alpha)
    println("* Near eval precompute")
    @time near_weights = mss.doublelayer_near_weights(dcurve, interior_points, alpha; ifgrad=true)
    hgrid = Lgrid/Ngrid_new

    # Setup test problem: BCs and RHS
    F1 = ffunc1.(X,Y)
    F2 = ffunc2.(X,Y)
    ubc1 = ufunc1.(xbdry, ybdry, alpha)
    ubc2 = ufunc2.(xbdry, ybdry, alpha)

    # Extend RHS
    #println("* PUX eval")
    tic()
    FEXT1, FEXT2 = PUX.pux_eval_vec(F1, F2, PUXStor)
    time_pux = toq()
    # Compute particular solution
    #println("* Periodic FFT Solve")
    tic()
    U1, U2, upbdry1, upbdry2, dU1dx, dU1dy, dU2dx, dU2dy =
        mss.per_modstokes(FEXT1, FEXT2, Lgrid, alpha, xbdry, ybdry, ifgrad=true, zerok0=true)
    time_fft = toq()
    # Estimate particular solution resolution
    UP_resolution =  max(estimate_fourier_resolution(U1),
                         estimate_fourier_resolution(U2))
    UPgrad_resolution = max( estimate_fourier_resolution(dU1dx),
                             estimate_fourier_resolution(dU1dy),
                             estimate_fourier_resolution(dU2dx),
                             estimate_fourier_resolution(dU2dy) )
    # Integral equation RHS
    rhs1 = ubc1 - upbdry1
    rhs2 = ubc2 - upbdry2
    rhs = vec([rhs1'; rhs2'])
    # Solve integral equation
    tic()
    if method==:fastdirect
        #println("* Fast direct solve")
        dsol = apply_compressed_matrix(Ainv, rhs)
        density = reshape(dsol, 2, :)
    elseif method==:fmm
        #println("* GMRES FMM solve")
        density = mss.gmres_solve(dcurve, rhs, alpha, flhs=flhs)
    elseif method==:dense_lu
        #println("* LU apply")
        dsol = MATU\(MATL\rhs[MATp])
        density = reshape(dsol, 2, :)
    elseif method==:dense_gmres
        println("* GMRES dense solve")
        @time dsol, gmlog = IterativeSolvers.gmres(MAT, rhs; tol=1e-15, restart=100, maxiter=100,
                                             verbose=false, log=true)
        if !gmlog.isconverged
            warn("GMRES did not converge, residual = $(gmlog.data[:resnorm][end])")
        end
        density = reshape(dsol, 2, :)
    end
    time_solve = toq()
    # Estimate density resolution
    L = mss.legendre_matrix(panelorder)
    density_resolution = 0.0;
    for i=1:numpanels
        xcoeff = L*density[1, (1:panelorder) + panelorder*(i-1)]
        ycoeff = L*density[2, (1:panelorder) + panelorder*(i-1)]
        xres = max(abs(xcoeff[end]), abs(xcoeff[end-1]))
        yres = max(abs(ycoeff[end]), abs(ycoeff[end-1]))
        res = max(xres, yres)
        density_resolution = max(density_resolution, res)
    end
    # Evaluate layer potential at interior points
    #println("* FMM eval")
    tic()
    u, u1grad, u2grad = mss.doublelayer_fast(dcurve, density, interior_points, alpha; specquad=true, weights=near_weights, ifgrad=true, precomp=fmm_precomp)
    time_fmm = toq()
    @. begin
        # Add hom sol to part sol on grid
        U1[interior]   += u[1,:]
        U2[interior]   += u[2,:]
        dU1dx[interior] += u1grad[1,:]
        dU1dy[interior] += u1grad[2,:]
        dU2dx[interior] += u2grad[1,:]
        dU2dy[interior] += u2grad[2,:]
        # Clear exterior
        U1[!interior] = 0.0
        U2[!interior] = 0.0
        dU1dx[!interior] = 0.0
        dU1dy[!interior] = 0.0
        dU2dx[!interior] = 0.0
        dU2dy[!interior] = 0.0
    end
    println(@sprintf("Timings: pux=%.3f, fft=%.3f, solve=%.3f, fmm=%.3f", time_pux, time_fft, time_solve, time_fmm))
    info(@sprintf("Resolution est.: den=%.1e, uP=%.1e, uP'=%.1e",
                  density_resolution, UP_resolution, UPgrad_resolution))
    if plotext
        figure(4)
        clf()
        pcolor(X, Y, FEXT1)
        plot(xbdry, ybdry,"r")
        axis("equal")
    end
    sol = mss.Solution(U1, U2, dU1dx, dU1dy, dU2dx, dU2dy,
                       copy(F1), copy(F2), copy(ubc1), copy(ubc2))
    return hgrid, errcheck(sol, X, Y, interior_points, interior, alpha)
end

function errcheck(sol, X, Y, interior_points, interior, alpha)
    ### END STUFF TO PUT IN TIMESTEP LOOP
    # Check errors in test problem
    println("* Computing error")
    # Reference
    xint = interior_points[1,:]
    yint = interior_points[2,:]
    uref1 = ufunc1.(xint, yint, alpha)
    uref2 = ufunc2.(xint, yint, alpha)
    numeval = size(interior_points, 2)

    uref = [uref1 uref2]'

    maxnorm(v) = maximum(abs.(v))
    l2norm(v) = sqrt(sum(vec(v.^2)))

    # Error
    @show size(sol.U1)
    @show size(interior)
    u = [sol.U1[interior] sol.U2[interior]]'
    err = u .- uref
    maxerr = Array{Float64}(numeval)
    for i=1:numeval
        maxerr[i] = norm(err[:, i], Inf)
    end
    max_rel_err = maxnorm(err) / maxnorm(uref)
    l2_rel_err = l2norm(err) / l2norm(uref)
    println("Max rel err on grid: ", max_rel_err)
    println("L2 rel err on grid: ",  l2_rel_err)

    println("* Derivative errors (max rel)")
    grad = max(     norm(vec(sol.dU1dx[interior]), Inf),
                    norm(vec(sol.dU1dy[interior]), Inf),
                    norm(vec(sol.dU2dx[interior]), Inf),
                    norm(vec(sol.dU2dy[interior]), Inf) )

    gradnorm = max( norm(vec(sol.dU1dx[interior]), Inf),
                    norm(vec(sol.dU1dy[interior]), Inf),
                    norm(vec(sol.dU2dx[interior]), Inf),
    norm(vec(sol.dU2dy[interior]), Inf) )

    grad = [sol.dU1dx[interior]; sol.dU1dy[interior]; sol.dU2dx[interior]; sol.dU2dy[interior]]
    gradref = [ufunc1x.(xint, yint, alpha)
               ufunc1y.(xint, yint, alpha)
               ufunc2x.(xint, yint, alpha)
               ufunc2y.(xint, yint, alpha)]

    Egrad_max = maxnorm(grad.-gradref) / maxnorm(gradref)
    Egrad_l2 = l2norm(grad.-gradref) / l2norm(gradref)

    @show Egrad_max
    @show Egrad_l2

    if ploterror
        println("* Displaying")
        figure(1, figsize=(4,3))
        clf()
        # Put on grid
        E = zeros(sol.U1)
        E[interior] = maxerr / maxnorm(uref)

        logE = log10.(E+1e-16)
        #interior_pcolor(X, Y, logE, interior, vmin=-16, vmax=0)
        interior_pcolor(X, Y, logE, interior, cmap=PyPlot.cm[:coolwarm])
        clim(log10(l2_rel_err)-1, log10(l2_rel_err))
        #plot(dcurve.points[1,:], dcurve.points[2,:], ".k")
        axis("off")
        axis("equal")
        xlim(0.6, 1.3)
        ylim(-0.35, 0.35)
        cbar = colorbar(pad=0.0, shrink=0.8)
        cbar[:set_label]("Relative error")
        tight_layout(0.2)
        savefig(@sprintf("figs/starfish_errfield_alpha%.3f.png", alpha))
    end

    return max_rel_err, l2_rel_err, Egrad_max, Egrad_l2
end

alpha_list = 10 .^ (-2.0:0.25:4.0)

Ntests = length(alpha_list)
max_rel_err = zeros(Ntests)
l2_rel_err = zeros(Ntests)
grad_max_rel_err = zeros(Ntests)
grad_l2_rel_err = zeros(Ntests)
hgrid = zeros(Ntests)

include("alpha_variation_plot.jl")

tic()
for i=1:Ntests
    hgrid[i], errors = solve(pux_R, alpha_list[i])
    max_rel_err[i], l2_rel_err[i], grad_max_rel_err[i], grad_l2_rel_err[i] = errors
    plot_results()
end
toq()

include("alpha_variation_plot.jl")
plot_results(paperplots=true)
savefig("../docs/paper/figs/alpha_variation.pdf")


text:
@testset "SBitstream" begin
    x = SBitstream(0.5)
    y = SBitstream(0.1)
    w = SBitstream(-0.1)
    X = SBitstream.(0.5 * rand(2, 2))
    Y = SBitstream.(-0.5 * rand(2, 2))
    v = SBitstream.([0.5 * rand(), -0.5 * rand()])

    @testset "FP Values" begin
        @test x.value == 0.5 && y.value == 0.1
        @test float(x) == x.value
        @testset for op in (+, -, *)
            @test float(op(x, y)) == op(float(x), float(y))
            @test float(op(x, w)) == op(float(x), float(w))
        end
        @test float(y / x) == float(y) / float(x)
        @test_logs (:warn, "SBitstream can only be ∈ [-1, 1] (saturation occurring).") float(x / y) == 1
        @test float(x ÷ 2) == float(x) / 2
        @test float(w ÷ 2) == float(w) / 2
        @test_throws ErrorException x ÷ 0.2
        @test float(sqrt(x)) == sqrt(float(x))
        @test float(decorrelate(x)) == float(x)
        @test float.(X * Y) == float.(X) * float.(Y)
        @test float.(X * v) == float.(X) * float.(v)
        @test float(norm(v)) == norm(float.(v))
        @test float(max(x, y)) == max(float(x), float(y))
    end

    @testset "Sample Generation" begin
        xbits = generate(x, 100000)
        @test mean(map(xbit -> pos(xbit) - neg(xbit), xbits)) ≈ x.value rtol = 0.1
    end

    @testset "Bit-Level Simulation" begin
        T = 50000
        @testset for op in (+, -, *, max)
            sim = simulatable(op, x, y)
            @test begin
                z = 0
                for t in 1:T
                    bit = pop!(sim(op, x, y))
                    z += pos(bit) - neg(bit)
                end

                z / T
            end ≈ float(op(x, y)) rtol = 0.1
            sim = simulatable(op, y, x)
            @test begin
                z = 0
                for t in 1:T
                    bit = pop!(sim(op, y, x))
                    z += pos(bit) - neg(bit)
                end

                z / T
            end ≈ float(op(y, x)) rtol = 0.1
            sim = simulatable(op, x, w)
            @test begin
                z = 0
                for t in 1:T
                    bit = pop!(sim(op, x, w))
                    z += pos(bit) - neg(bit)
                end

                z / T
            end ≈ float(op(x, w)) rtol = 0.1
        end
        @testset "op = ÷" begin
            sim = simulatable(÷, x, 2)
            @test begin
                z = 0
                for t in 1:T
                    bit = pop!(sim(÷, x, 2))
                    z += pos(bit) - neg(bit)
                end

                z / T
            end ≈ float(x ÷ 2) rtol = 0.1
            sim = simulatable(÷, w, 2)
            @test begin
                z = 0
                for t in 1:T
                    bit = pop!(sim(÷, w, 2))
                    z += pos(bit) - neg(bit)
                end

                z / T
            end ≈ float(w ÷ 2) rtol = 0.1
        end
        @testset "op = /" begin
            sim = simulatable(/, y, x)
            @test begin
                z = 0
                for t in 1:T
                    bit = pop!(sim(/, y, x))
                    z += pos(bit) - neg(bit)
                end

                z / T
            end ≈ float(y / x) rtol = 0.1
        end
        @testset for op in (sqrt, decorrelate)
            sim = simulatable(op, x)
            @test begin
                z = 0
                for t in 1:T
                    bit = pop!(sim(op, x))
                    z += pos(bit) - neg(bit)
                end

                z / T
            end ≈ float(op(x)) rtol = 0.1
        end
        @testset "op = max(x, y, z)" begin
            sim = simulatable(max, x, y, w)
            @test begin
                z = 0
                for t in 1:T
                    bit = pop!(sim(max, x, y, w))
                    z += pos(bit) - neg(bit)
                end

                z / T
            end ≈ float(max(x, y, w)) rtol = 0.1
        end
        @testset "op = * (matrix, matrix)" begin
            sim = simulatable(*, X, Y)
            @test begin
                Z = zeros(2, 2)
                for t in 1:T
                    bit = pop!.(sim(*, X, Y))
                    Z .+= pos.(bit) .- neg.(bit)
                end

                Z ./ T
            end ≈ float.(X * Y) rtol = 0.1
            sim = simulatable(*, Y, X)
            @test begin
                Z = zeros(2, 2)
                for t in 1:T
                    bit = pop!.(sim(*, Y, X))
                    Z .+= pos.(bit) .- neg.(bit)
                end

                Z ./ T
            end ≈ float.(Y * X) rtol = 0.1
        end
        @testset "op = * (matrix, vector)" begin
            sim = simulatable(*, X, v)
            @test begin
                T = 100000
                Z = zeros(2)
                for t in 1:T
                    bit = pop!.(sim(*, X, v))
                    Z .+= pos.(bit) .- neg.(bit)
                end

                Z ./ T
            end ≈ float.(X * v) rtol = 0.1
            sim = simulatable(*, transpose(v), X)
            @test begin
                T = 100000
                Z = zeros(1, 2)
                for t in 1:T
                    bit = pop!.(sim(*, transpose(v), X))
                    Z .+= pos.(bit) .- neg.(bit)
                end

                Z ./ T
            end ≈ float.(transpose(v) * X) rtol = 0.1
        end
        @testset "op = norm" begin
            sim = simulatable(norm, v)
            @test begin
                T = 100000
                z = 0
                for t in 1:T
                    bit = pop!(sim(norm, v))
                    z += pos(bit) - neg(bit)
                end

                z / T
            end ≈ float(norm(v)) rtol = 0.1
        end
    end
end


text:
module DTW

function dtw(x,y; distance::Function=Distance.default)
	D = zeros(Float64, length(x)+1, length(y)+1)
	D[1,2:end] = Inf
	D[2:end,1] = Inf
	for i=1:length(x), j=1:length(y)
		D[i+1, j+1] = distance(x[i], y[j])
	end
	for i=1:length(x), j=1:length(y)
		D[i+1, j+1] += min(D[i,j], D[i,j+1], D[i+1,j])
	end

	D = D[2:end, 2:end]
	dist = D[end,end]/sum(size(D))
	dist, D, trackback(D)
end

function trackback(D)
	i,j = size(D)
	p,q = [i],[j]

	while i > 1 && j > 1
		tb = indmin([D[i-1,j-1], D[i-1,j], D[i,j-1]])
		@show(i,j,tb)
		tb in [1,2] && (i-=1)
		tb in [1,3] && (j-=1)
		unshift!(p,i)
		unshift!(q,j)
	end
	unshift!(p,1)
	unshift!(q,1)
	p,q
end

module Distance # contains distance measurement functions
default(x,y) = norm(x-y,1)
end

export dtw

end # module


text:
using Flux, Printf, Statistics, IPMeasures, LinearAlgebra, MLDataPattern, FluxExtensions
using IPMeasures: mmd, IMQKernel
using FluxExtensions: log_normal
using Plots

data(n) = Float32.(cholesky([1 0.9; 0.9 1]).U * randn(2,n))

model = (f = Chain(Dense(2,10,tanh),Dense(10,1)), g = Chain(Dense(1,10,tanh), Dense(10,1)), μ = Dense(1,2), σ = Chain(Dense(1,10,tanh),Dense(10,2)))
imq(z, c = 0.01f0) = mmd(IMQKernel(0.01f0), randn(Float32,1,size(z,2)), z)

lkl_full(model, x, z) = FluxExtensions.log_normal(x, model.μ(model.g(z)), softplus.(model.σ(model.g(z))))
lkl_fixed(model, x, z) = FluxExtensions.log_normal(x, model.μ(model.g(z)))

function closestz(model, x, steps = 100)
z = param(Flux.data(model.f(x)))
ps = Flux.Tracker.Params([z])
opt = ADAM()
_lkl(model, x, z) = mean(lkl(model, x, z) .- FluxExtensions.log_normal(z))
_info() = println("likelihood = ", _lkl(model, x, z))
li = Flux.data(_lkl(model, x, z))
# Flux.train!((i) -> -_lkl(model, x, z), ps, 1:steps, opt, cb = () -> Flux.throttle(_info(),5))
Flux.train!((i) -> -_lkl(model, x, z), ps, 1:steps, opt, cb = () -> Flux.throttle(_info(),5))
le = Flux.data(_lkl(model, x, z))
println("initial = ",li, " final = ",le)
Flux.data(z)
end

function manifoldz(model, x, steps = 100)
z = param(Flux.data(model.f(x)))
ps = Flux.Tracker.Params([z])
opt = ADAM()
_info() = println("likelihood = ", lkl(model, x, z))
li = Flux.data(mean(lkl(model, x, z)))
# Flux.train!((i) -> -mean(lkl(model, x, z)), ps, 1:steps, opt, cb = () -> Flux.throttle(_info(),5))
Flux.train!((i) -> -mean(lkl(model, x, z)), ps, 1:steps, opt)
le = Flux.data(mean(lkl(model, x, z)))
println("initial = ",li, " final = ",le)
Flux.data(z)
end

lkl = lkl_full

function loss(model, x)
z = model.f(x)
imq(z) - mean(lkl(model, x, z))
end

function loss(model, x, ϵ)
z = model.f(x .+ ϵ .* randn(eltype(x),size(x)))
imq(z) - mean(lkl(model, x, z))
end

function info(model, x)
z = model.f(x)
println("imq = ", imq(z), " lkl = ", mean(lkl(model, x, z)))
end

xv = data(1000)
x = data(10000)

ps = Flux.Tracker.Params()
foreach(m -> push!(ps,params(m)...),model)
opt = ADAM()
Flux.train!(x -> loss(model, getobs(x), 0.05), ps, RandomBatches((x,), 100,1000), opt, cb = () -> Flux.throttle(info(model, xv),5))

#histogram of latent data
frmodel = FluxExtensions.freeze(model)
z = frmodel.f(xv)
t2m(x) = reshape([x[1],x[2]],:,1)
xr = range(minimum(xv[1,:]) - 1 , maximum(xv[1,:])+ 1 , length = 100)
yr = range(minimum(xv[2,:]) - 1 , maximum(xv[2,:])+ 1 , length = 100)
histogram(z[:],bins = round(Int, sqrt(length(z))))

#sample reconstructed data
xx = frmodel.μ(frmodel.g(z)) .+ randn(2,size(z,2)) .* softplus.(frmodel.σ(frmodel.g(z)))
# xx = frmodel.μ(frmodel.g(z))
scatter(xx[1,:],xx[2,:])

#generate new data
z = randn(1,1000)
xx = frmodel.μ(frmodel.g(z)) .+ randn(2,size(z,2)) .* softplus.(frmodel.σ(frmodel.g(z)))
scatter(xx[1,:],xx[2,:])

#likelihood of latent layer
z = frmodel.f(xv)
contour(xr, yr, (x,y) -> log_normal(frmodel.f([x,y]))[1], title = "p(z)");
scatter!(xv[1,:], xv[2,:], zcolor = log_normal(z)[:], title = "p(z)")

#likelihood of output layer
contour(xr, yr, (x...) -> lkl(frmodel, t2m(x), frmodel.f(t2m(x)))[1], title = "p(x|z)");
scatter!(xv[1,:], xv[2,:], zcolor = lkl(frmodel, xv, z)[:], title = "p(x|z)")

#joined likelihood
contour(xr, yr, (x...) -> log_normal(frmodel.f(t2m(x)))[1] + lkl(frmodel, t2m(x),  frmodel.f(t2m(x)))[1], title = "p(x|z) + p(z)");
scatter!(xv[1,:], xv[2,:], zcolor = log_normal(z)[:] + lkl(frmodel, xv, z)[:], title = "p(x|z) + p(z)")

#optimized likelihood
function opt_lkl(model, x)
z = closestz(model, x)
log_normal(z) + lkl(frmodel, x, z)
end
contour(xr, yr, (x...) -> opt_lkl(frmodel, t2m(x))[1], title = "optimized");
scatter!(xv[1,:], xv[2,:], zcolor = opt_lkl(frmodel, xv)[:], title = "optimized")
heatmap(xr, yr, (x...) -> opt_lkl(frmodel, t2m(x))[1], title = "optimized")

#optimized reconstruction error
function opt_lkl(model, x)
z = manifoldz(model, x)
log_normal(z) + lkl(frmodel, x, z)
end
contour(xr, yr, (x...) -> opt_lkl(frmodel, t2m(x))[1], title = "optimized");
scatter!(xv[1,:], xv[2,:], zcolor = opt_lkl(frmodel, xv)[:], title = "optimized")
heatmap(xr, yr, (x...) -> opt_lkl(frmodel, t2m(x))[1], title = "optimized")


text:
# Renders and Films
export render_dir,
       render_user_dir,
       render_backend_dir,
       render_kind_dir,
       render_color_dir,
       render_ext,
       render_width,
       render_height,
       render_quality,
       render_exposure,
       set_render_dir,
       render_size,
       prepare_for_saving_file,
       render_pathname,
       render_view,
       rendering_with

# There is a render directory
const render_dir = Parameter(homedir())
# with a user-specific subdirectory
const render_user_dir = Parameter(".")
# with a backend-specific subdirectory
const render_backend_dir = Parameter(".")
# and with subdirectories for static images, movies, etc
const render_kind_dir = Parameter("Render")
# and with subdirectories for white, black, and colored renders
const render_color_dir = Parameter(".")
# containing files with different extensions
const render_ext = Parameter(".png")

render_pathname(name::String) =
    normpath(
      joinpath(
        render_dir(),
        render_user_dir(),
        render_backend_dir(),
        render_kind_dir(),
        render_color_dir(),
        "$(name)$(render_ext())"))

const render_width = Parameter(1024)
const render_height = Parameter(768)
const render_quality = Parameter{Real}(0) # [-1, 1]
const render_exposure = Parameter{Real}(0)  # [-3, +3]
const render_floor_width = Parameter(1000)
const render_floor_height = Parameter(1000)

render_size(width::Integer, heigth::Integer) =
    (render_width(width), render_height(heigth))

prepare_for_saving_file(path::String) =
    let p = normpath(path)
        mkpath(dirname(p))
        rm(p, force=true)
        p
    end

export film_active,
       film_filename,
       film_frame,
       start_film,
       frame_filename,
       save_film_frame,
       saving_film_frames

const film_active = Parameter(false)
const film_filename = Parameter("")
const film_frame = Parameter(0)
const saving_film_frames = Parameter(true)

start_film(name::String) =
    begin
        film_active(true)
        film_filename(name)
        film_frame(0)

    end

frame_filename(filename::String, i::Integer) =
    "$(filename)-frame-$(lpad(i,3,'0'))"

save_film_frame(obj::Any=true) =
  begin
    if saving_film_frames()
      with(render_kind_dir, "Film") do
        render_view(prepare_for_saving_file(render_pathname(frame_filename(film_filename(), film_frame()))),
                    current_backend())
        film_frame(film_frame() + 1)
      end
    end
    obj
  end

rendering_with(f;
    dir=render_dir(),
    user_dir=render_user_dir(),
    backend_dir=render_backend_dir(),
    kind_dir=render_kind_dir(),
    color_dir=render_color_dir(),
    ext=render_ext(),
    width=render_width(),
    height=render_height(),
    quality=render_quality(),
    exposure=render_exposure(),
    floor_width=render_floor_width(),
    floor_height=render_floor_height()) =
    with(render_dir, dir) do
        with(render_user_dir, user_dir) do
            with(render_backend_dir, backend_dir) do
                with(render_kind_dir, kind_dir) do
                    with(render_color_dir, color_dir) do
                        with(render_ext, ext) do
                            with(render_width, width) do
                                with(render_height, height) do
                                    with(render_quality, quality) do
                                        with(render_exposure, exposure) do
                                           with(render_floor_width, floor_width) do
                                               with(render_floor_height, floor_height) do
                                                   f()
                                               end
                                            end
                                        end
                                    end
                                end
                            end
                        end
                    end
                end
            end
        end
    end

export default_lens,
       track_still_target,
       track_moving_target,
       walkthrough,
       panning,
       lens_zoom,
       dolly_effect_back,
       dolly_effect_forth,
       set_view_save_frame,
       set_view_save_frames,
       interpolate_view_save_frames,
       select_camera_target_lens_positions

default_lens = Parameter(50)

set_view_save_frame(camera, target, lens=default_lens()) =
  begin
    set_view(camera, target, lens)
    save_film_frame()
  end

# Tracking

# Rotation around target

# target object remains still and the camera moves around it
# parameters: camera path (list of locations) + fixed target (1 location)
# number of locations in the path defines the number of frames

track_still_target(camera_path, target, lens=default_lens()) =
  for camera in camera_path
    set_view_save_frame(camera, target, lens)
  end

# Moving camera and target

# target object moves and the camera makes a similar movement
# parameters: camera (initial location) + target path (list of locations)
# number of locations in the path defines the number of frames

set_view_save_frames(cameras::Locs, targets::Locs, lens::Real) =
  for (camera, target) in zip(cameras, targets)
    set_view_save_frame(camera, target, lens)
  end

set_view_save_frames(cameras::Locs, targets::Locs, lenses::Vector{<:Real}) =
  for (camera, target, lens) in zip(cameras, targets, lenses)
    set_view_save_frame(camera, target, lens)
  end

twin_path(path, v) = map(p -> p+v, path)

track_moving_target(camera, targets, lens=default_lens()) =
  let v = camera - targets[1]
      cameras = twin_path(targets, v)
    set_view_save_frames(cameras, targets, lens)
  end

# Walkthroughs

# targets object moves and the camera follows a few positions behind
# parameters: path (list of locations) + number of distance locations

walkthrough(path, camera_spread, lens=default_lens()) =
  let target_path = path[camera_spread:end]
      camera_path = path[1:end-camera_spread]
    set_view_save_frames(camera_path, target_path, lens)
  end

# Panning

# Moving target, still camera

# Camera rotation
# parameters: fixed camera + target path
# number of locations in the path defines the number of frames

panning(camera, targets, lens=default_lens()) =
  for target in targets
    set_view_save_frame(camera, target, lens)
  end

#;(define my-path (rotation-path (x 0) -50 0 2pi 0 100)) ;example
#;(panning (x 0) my-path) ;example

# Zooming

# Lens zoom

# Changing lens focal length while camera position remains the same
# paremeters: fixed camera + fixed target + zomming scale (lens degrees to change) + number of frames
# zomming scale: delta > 1 => zooms in, delta < 1 => zooms out - |delta| >> zoom mais violeto

lens_zoom(camera, target, delta, frames, lens=default_lens()) =
  let lenses = range(lens, step=delta, length=frames)
    for lens in lenses
      set_view_save_frame(camera, target, lens)
    end
  end

# Dolly zoom

# Moving camera backwards as the lens increases the focal length
# parameters: camera start + fixed target + zomming scale (lens degrees to change) + number of frames

#=
dolly_effect(cur_camera, target, cur_lens, new_camera) =
  let cur_dist = distance(cur_camera, target)
      new_dist = distance(new_camera, target)
      new_lens = (cur_lens*new_dist)/cur_dist
    set_view_save_frame(new_camera, target, new_lens)
  end

dolly_effect_pull_push(back_forth, delta) =
  camera, target, lens = view()
  let d = distance(camera, target)
    let (new_camera_back, new_camera_forth) = (target+(camera-target)*(d+delta)/d, target+(camera-target)*(d-delta)/d)
      back/forth == 0 ? dolly_effect(camera, target, lens, new_camera_back) : dolly_effect(camera, target, lens, new_camera_forth)
    end
  end

dolly_zoom(back/forth, camera_start, target, lens, delta, frames) =
  begin
    set_view_save_frame(camera_start, target, lens)
    for i in 0:frames-1
      dolly_effect_pull/push(back/forth, delta)
    end
  end

current_view_dolly_zoom(back/forth, delta, frames) =
  for i in 0:frames-1
    dolly_effect_pull/push(back/forth, delta)
  end

=#
#;(view (xyz 5 -0.8 0.5) (xyz 5 0 0.5) 10)
#;(dolly-zoom (xyz 5 -0.8 0.5) (xyz 5 0 0.5) 0.05 100) ;example

# When we don't have access to the parameters of the current view
# parameters: zomming scale + camera start + fixed target + initial lens value + number of frames

dolly_effect_back(delta, camera, target, lens, frames) =
  let d = distance(camera, target)
      new_camera = target+(camera-target)*(d+delta)/d
      new_d = distance(new_camera, target)
      new_lens = (lens*new_d)/d
    if frames <= 0
      nothing
    else
      set_view_save_frame(camera, target, lens)
      dolly_effect_back(delta, new_camera, target, new_lens, frames-1)
    end
  end

dolly_effect_forth(delta, camera, target, lens, frames) =
  let d = distance(camera, target)
      new_camera = target+(camera-target)*(d-delta)/d
      new_d = distance(new_camera, target)
      new_lens = (lens*new_d)/d
    if frames <= 0
      nothing
    else
      set_view_save_frame(camera, target, lens)
      dolly_effect_forth(delta, new_camera, target, new_lens, frames-1)
    end
  end

#=
To collect the camera/target/lens path
=#

select_camera_target_lens_positions() =
  select_position() == nothing ?
    [] :
    [get_view(), select_camera_target_lens_path()...]

interpolate_view_save_frames(ctl_positions, nframes) =
  set_view_save_frames(
    map_division(
      in_world,
      open_spline_path(map(p -> p[1], ctl_positions)), nframes),
    map_division(
      in_world,
      open_spline_path(map(p -> p[2], ctl_positions)), nframes),
    map(p -> p[3], ctl_positions))

#=
To generate the actual film from the set of frames, we need ffmpeg.
=#

create_mp4_from_frames(name=film_filename()) =
  with(render_kind_dir, "Film") do
    with(film_filename, name) do
      with(render_ext, "-film.mp4") do
        let film_path = prepare_for_saving_file(render_pathname(name))
          with(render_ext, "-frame*.png") do
            let frames_path = render_pathname(name)
              println(frames_path)
              println(film_path)
              println("`ffmpeg -pattern_type glob -i '$(frames_path)' -c:v libx264 -vf fps=25 -pix_fmt yuv420p $(film_path)`")
            end
          end
        end
      end
    end
  end


text:
println("Starting tests with $(Threads.nthreads()) threads out of `Sys.CPU_THREADS = $(Sys.CPU_THREADS)`...")
using CheapThreads, Aqua, ForwardDiff
using Test

@testset "Range Map" begin
    function rangemap!(f::F, allargs, start, stop) where {F}
        dest = first(allargs)
        args = Base.tail(allargs)
        @inbounds @simd for i ∈ start:stop
            dest[i] = f(Base.unsafe_getindex.(args, i)...)
        end
        nothing
    end

    function tmap!(f::F, args::Vararg{AbstractArray,K}) where {K,F}
        dest = first(args)
        N = length(dest)
        mapfun! = (allargs, start, stop) -> rangemap!(f, allargs, start, stop)
        batch(mapfun!, (N, num_threads()), args...)
        dest
    end

    x = rand(1024); y = rand(length(x)); z = similar(x);
    foo(x,y) = exp(-0.5abs2(x-y))
    println("Running `tmap!` test...")
    @test tmap!(foo, z, x, y) ≈ foo.(x, y)

    # TODO: why does this crash Julia?
    # function slow_task!((x, digits, n), j, k)
    #     start = 1 + (n * (j - 1)) ÷ num_threads()
    #     stop =  (n* k) ÷ num_threads()
    #     target = 0.0
    #     for i ∈ start:stop
    #         target += occursin(digits, string(i)) ? 0.0 : 1.0 / i
    #     end
    #     x[1,j] = target
    # end
    # function slow_cheap(n, digits)
    #     x = zeros(8, num_threads())
    #     batch(slow_task!, (num_threads(), num_threads()), x, digits, n)
    #     sum(@view(x[1,1:end]))
    # end
    # slow_cheap(1000, "9")
end

@testset "ForwardDiff" begin
    x = randn(80);
    dxref = similar(x);
    dx = similar(x);
    f(x) = -sum(sum ∘ sincos, x)
    println("Running threaded ForwardDiff test...")
    CheapThreads.threaded_gradient!(f, dx, x, ForwardDiff.Chunk(8))
    ForwardDiff.gradient!(dxref, f, x, ForwardDiff.GradientConfig(f, x, ForwardDiff.Chunk(8), nothing))
    @test dx == dxref
end

println("Package tests complete. Running `Aqua` checks.")
Aqua.test_all(CheapThreads)