text:
using StaticArrays

abstract type AbstractProbe end

abstract type AbstractProbe1d <: AbstractProbe end
abstract type AbstractProbe2d <: AbstractProbe end

abstract type AbstractProbe3d <: AbstractProbe end

struct HWSupport
    "Support model"
    model::String
    "Support tag (storage and control)"
    tag::String
    "Nominal resistance of the support"
    R::Float64
    "Diameter of the support"
    D::Float64
    "Length of the support"
    L::Float64
    "Length of support cables"
    Lc::Float64
end

"""
    `HWSupport(;model="55H21", tag="", R=0.44, D=4.0, L=235, Lc=765)`

Stores information about Probe support. Noty strictly necessary but
could be useful for future reference

"""
HWSupport(;model="", tag="", R=0.0, D=0.0, L=0.0, Lc=0.0) = HWSupport(model, tag, R, D, L, Lc)

struct HWCable
    "Cable model/type"
    model::String
    "Cable tag (storage and control)"
    tag::String
    "Nominal resistance of the cable"
    R::Float64
    "Length of the cable"
    L::Float64
    "Connector types (left and right)"
    conn::Tuple{String,String}
end

"""
    `HWCable(model="", tag="", R=0.2, L=4.0, ("BNC","BNC"))`

Stores information about the cable used. Not strictly necessary but
could be useful.

"""
HWCable(;model="A1863 - 50Ω impedance", tag="", R=0.2, L=4.0, conn=("BNC", "BNC")) = HWCable(model, tag, R, L, conn)

struct HWBridge
    "Anemometer system model"
    model::String
    "Anemometer tag"
    tag::String
    "Specific bridge id"
    id::Int
    "Output offset"
    offset::Float64
    "Output gain"
    gain::Float64
    "Output low pass filter"
    lpfilt::Float64
    "Bridge ration"
    bridgeratio::String
    "Bridge amplifier gain"
    ampgain::Float64
    "Bridge filter"
    filter::Float64
    "Cable compensation"
    cablecomp::Int
    "Cooling interval"
    coolinterv::Float64
end

"""
    `HWBridge(model="", tag="", id=1, offset=0, gain=1, lpfilt=1e3,
              bridgeration="20:1", ampgain=8, filter=3,
              cablecomp=0, coolinterv=2)`

Information about hotwire hardware system.

"""
HWBridge(;model="Streamline", tag="", id=1, offset=0.0, gain=1.0,
         lpfilt=1e3, bridgeratio="20:1", ampgain=8, filter=3,
         cablecomp=0, coolinterv=2.0) =
             HWBridge(model, tag, id, offset, gain, lpfilt,
                      bridgeratio, ampgain, filter, cablecomp,
                      coolinterv)

mutable struct Probe1d{Anem<:AbstractThermalAnemometer,Calibr} <: AbstractProbe1d
    "Anemometer sensor information"
    sensor::Anem
    "Calibration data"
    cal::Calibr
    "Probe model"
    model::String
    "Probe tag"
    tag::String
    "Probe leads resistance"
    RL::Float64
    "Hardware information"
    bridge::HWBridge
    "Probe support"
    support::HWSupport
    "Connecting cable"
    cable::HWCable
end

"""
    `Probe1d(sensor, cal; model="model", tag="HW 1d", RL=0.5,
             bridge=HWBridge(), support=HWSupport(), cable=HWCable())`

Creates a 1d probe. This completely characterizes a 1d probe. It includes
information on calibration, bridge configuration, output filter, etc.

# Arguments:

 - `sensor`: an `AbstractThermalAnemometer` object
 - `cal`: `CalibrCurve` object
 - `model`: Probe model
 - `tag`: Probe id tag
 - `RL`: Lead resistance
 - `bridge`: Bridge configuration
 - `support`: `HWSupport` object specifying the probe support
 - `cable`: `HWCable` object with information on cable

"""
function Probe1d(sensor, cal; model="", tag="", RL=0.5,
                 bridge=HWBridge(id=1),
                 support=HWSupport(), cable=HWCable())
    return Probe1d(sensor, cal, model, tag, RL, bridge, support, cable)
end

(anem::Probe1d)(E, temp) = anem.cal(anem.sensor, E, temp)
(anem::Probe1d)(E) = anem.cal(anem.sensor, E, anem.cal.T0)

"""
    `velocity!(anem, E, T, 1)`

Calculate, in place, the velocity from an anemometer. In this function,
the input is an abstract array where 1 column with the bridge output is given by argument `idx` represent the output of the anemometer. A unique temperature should be used.

"""
function velocity!(anem::Probe1d, E::AbstractMatrix, T, idx=1)

    npts = size(E,1)
    for i in 1:npts
        U = anem(E[i,idx], T)
        E[i,idx] = U
    end
end

"""
    `velocity(anem, E, T, 1)`

Calculate the velocity from an anemometer. In this function,
the input is an abstract array where 1 column with the bridge output is given by argument `idx` represent the output of anemometer. A unique temperature should be used

"""
function velocity(anem::Probe1d, E::AbstractMatrix, T, idx=1)

    npts = size(E,1)
    U = zeros(npts)
    for i in 1:npts
        U[i] = anem(E[i,idx], T)
    end
    return U
end

mutable struct Probe2d{Anem<:AbstractThermalAnemometer,Calibr} <: AbstractProbe2d
    "Sensor information"
    sensor::NTuple{2,Anem}
    "Sensor calibration"
    cal::NTuple{2,Calibr}
    "k² directional calibration for each sensor"
    k²::NTuple{2,Float64}
    "h² directional calibration for each sensor"
    h²::NTuple{2,Float64}
    "Matrix with the cosine of each wire with respect to x, y and z axes"
    cosϕ::Matrix{Float64}
    "Probe model"
    model::String
    "Probe tag"
    tag::String
    "Probe leads resistance"
    RL::NTuple{2,Float64}
    "Hardware information"
    bridge::NTuple{2,HWBridge}
    "Probe support"
    support::HWSupport
    "Connecting cable"
    cable::NTuple{2,HWCable}

end

mutable struct Probe3d{Anem<:AbstractThermalAnemometer,Calibr} <: AbstractProbe3d
    "Sensor information"
    sensor::NTuple{3,Anem}
    "Sensor calibration"
    cal::NTuple{3,Calibr}
    "k² directional calibration for each sensor"
    k²::NTuple{3,Float64}
    "h² directional calibration for each sensor"
    h²::NTuple{3,Float64}
    "Cosine of wire direction with respect to axes x, y and z"
    cosϕ::SMatrix{3,3,Float64,9}
    "Probe model"
    model::String
    "Probe tag"
    tag::String
    "Probe leads resistance"
    RL::NTuple{3,Float64}
    "Factor to calculate effective velocity from calibration velocity"
    c2e::SVector{3,Float64}
    "Matrix to calculate velocity in wire coordinates from effective vel."
    A::SMatrix{3,3,Float64,9}
    "Hardware information"
    bridge::NTuple{3,HWBridge}
    "Probe support"
    support::HWSupport
    "Connecting cable"
    cable::NTuple{3,HWCable}
end

"""
    `Probe3d(sensor, cal, k², h²; ϕ=ϕ, model="model", tag="HW 1d", RL=0.5,
             bridge=br, support=su, cable=cc)`

Creates a 2d probe. This completely characterizes a 1d probe. It includes
information on calibration, directional calibration, bridge configuration,
 output filter, etc.
The effective velocity acting on a wire is given by
``
(Uᵉ)² = U₁² + h²⋅U₂² + k²⋅U₃²
``
In this equation, Uᵉ is the effective velocity, U₁ is the velocity normal
to the wire and in the direction of the probe body, U₂ is the velocity
normal to the flow but also normal to the body and U₃ is the velocity
along the wire.

The function assumes that the calibration was carried out with flow along
the x axis.

# Arguments:

 - `sensor`: a tuple of `AbstractThermalAnemometer` objects (3)
 - `cal`: `CalibrCurve` object
 - `model`: Probe model
 - `tag`: Probe id tag
 - `RL`: Lead resistance
 - `bridge`: Bridge configuration
 - `support`: `HWSupport` object specifying the probe support
 - `cable`: `HWCable` object with information on cable

"""
function Probe3d(sensor, cal, k², h²; ϕ=[125.264 45 114.094;
                                         125.264 135 114.094;
                                         125.264 90 35.264],
                 model="55R91", tag="", RL=0.5,
                 bridge=(HWBridge(), HWBridge(), HWBridge()),
                 support=HWSupport(),
                 cable=(HWCable(), HWCable(), HWCable()))

    # Calculate the cosine of the angles
    cosϕ = SMatrix{3,3,Float64}(cosd.(ϕ))
    c2e = SVector{3,Float64}( (1.0 + k²[1] + h²[1])*cosϕ[1,1]^2,
                              (1.0 + k²[2] + h²[2])*cosϕ[2,1]^2,
                              (1.0 + k²[3] + h²[3])*cosϕ[3,1]^2 )
    A = SMatrix{3,3,Float64}([k²[1] 1.0      h²[1];
                              h²[2] k²[2]    1.0;
                              1.0   h²[3]  k²[3]])
    return Probe3d((sensor[1], sensor[2], sensor[3]),
                   (cal[1], cal[2], cal[3]), (k²[1], k²[2], k²[3]),
                   (h²[1], h²[2], h²[3]), cosϕ, model, tag, (RL,RL,RL),
                   c2e, inv(A), (bridge[1],bridge[2],bridge[3]),
                   support, (cable[1],cable[2],cable[3]))
end

function (anem::Probe3d)(E₁, E₂, E₃, T)

    # Calibration curve and temperature correction for
    # each sensor and calculation of effective velocity for each wire
    Uᵉ = SVector(anem.cal[1](anem.sensor[1], E₁, T)^2 * anem.c2e[1],
                 anem.cal[2](anem.sensor[2], E₂, T)^2 * anem.c2e[2],
                 anem.cal[3](anem.sensor[3], E₃, T)^2 * anem.c2e[3])

    # Calculate speed in wire coordinates
    Uʷ = anem.A * Uᵉ

    U1 = sqrt(max(Uʷ[1], 0.0))
    U2 = sqrt(max(Uʷ[2], 0.0))
    U3 = sqrt(max(Uʷ[3], 0.0))

    return (-U1*anem.cosϕ[1,1] - U2*anem.cosϕ[2,1] - U3*anem.cosϕ[3,1],
            -U1*anem.cosϕ[1,2] - U2*anem.cosϕ[2,2] - U3*anem.cosϕ[3,2],
            -U1*anem.cosϕ[1,3] - U2*anem.cosϕ[2,3] - U3*anem.cosϕ[3,3])
end

"""
    `velocity!(anem, E, T, 1:3)`

Calculate, in place, the velocity from an anemometer. In this function,
the input is an abstract array where 3 columns given by argument `idx`
represent the output of each anemometer. A unique temperature should be
used

"""
function velocity!(anem::Probe3d, E::AbstractMatrix, T, idx=1:3)

    npts = size(E,1)
    for i in 1:npts
        U, V, W = anem(E[i,idx[1]], E[i,idx[2]], E[i,idx[3]], T)
        E[i,idx[1]], E[i,idx[2]], E[i,idx[3]] = U, V, W
    end

end

"""
    `velocity(anem, E, T, 1:3)`

Velocity calculation for arrays. A new array where each column represents
a component of the velocity is created

"""
function velocity(anem::Probe3d, E::AbstractMatrix, T, idx=1:3)
    npts = size(E,1)

    Uout = zeros(npts, 3)

    for i in 1:npts
        U, V, W = anem(E[i,idx[1]], E[i,idx[2]], E[i,idx[3]])
        Uout[i,idx[1]], Uout[i,idx[2]], Uout[i,idx[3]] = U, V, W
    end

    return Uout
end


text:
## gcmfaces type definition + methods

## type definition

abstract type AbstractGcmfaces{T, N} <: AbstractArray{T, N} end

"""
    gcmfaces{T, N}

gcmfaces data structure.
"""
struct gcmfaces{T, N} <: AbstractGcmfaces{T, N}
   nFaces::Int
   grTopo::String
   f::Array{Array{T,N},1}
   fSize::Array{NTuple{N, Int}}
   aSize::NTuple{N, Int}
end

"""
    gcmsubset{T, N}

gcmsubset data structure.
"""
struct gcmsubset{T, N} <: AbstractGcmfaces{T, N}
   nFaces::Int
   grTopo::String
   f::Array{Array{T,N},1}
   fSize::Array{NTuple{N, Int}}
   aSize::NTuple{N, Int}
   i::Array{Array{T,N},1}
   iSize::Array{NTuple{N, Int}}
end

## additional constructors for gcmfaces

"""
    gcmfaces(nFaces::Int,grTopo::String,::Type{T},fSize::Array{NTuple{N, Int}},
      aSize::NTuple{N,Int}) where {T,N}

gcmfaces constructor.
"""
function gcmfaces(nFaces::Int,grTopo::String,::Type{T},
  fSize::Array{NTuple{N, Int}},
  aSize::NTuple{N,Int}) where {T,N}
  f=Array{Array{T,N},1}(undef,nFaces)
  for a=1:nFaces
    f[a]=Array{T}(undef,fSize[a])
  end
  gcmfaces{T,N}(nFaces,grTopo,f,fSize,aSize)
end

"""
    gcmfaces(nFaces::Int,grTopo::String,
      v1::Array{Array{T,N},1}) where {T,N}

gcmfaces constructor.
"""
function gcmfaces(nFaces::Int,grTopo::String,
  v1::Array{Array{T,N},1}) where {T,N}
  fSize=fsize(v1)
  aSize=fsize(v1,0)
  gcmfaces{T,N}(nFaces,grTopo,deepcopy(v1),fSize,aSize)
#  gcmfaces(nFaces,grTopo,T,fs,as)
end

"""
    gcmfaces(A::AbstractGcmfaces{T, N}) where {T,N}

gcmfaces constructor.
"""
function gcmfaces(A::AbstractGcmfaces{T, N}) where {T,N}
  #should this be called similar? deepcopy?
  fSize=fsize(A)
  aSize=size(A)
  gcmfaces{T,N}(nFaces,grTopo,deepcopy(A.f),fSize,aSize)
#  gcmfaces(nFaces,grTopo,T,fSize,aSize)
end

"""
    gcmfaces()

gcmfaces constructor.
"""
function gcmfaces()
  if isdefined(MeshArrays,:nFaces)
    nFaces=MeshArrays.nFaces
    grTopo=MeshArrays.grTopo
    T=MeshArrays.ioPrec
    fSize=MeshArrays.facesSize
    aSize=(prod(MeshArrays.ioSize),1)
  else
    nFaces=5
    grTopo="llc"
    T=Float64
    fSize=[(90, 270), (90, 270), (90, 90), (270, 90), (270, 90)]
    aSize=(105300, 1);
  end
  gcmfaces(nFaces,grTopo,T,fSize,aSize)
end

## additional constructors for gcmsubset

#maybe: replace this constructor with one that gets A and sets f to view(A.f)
function gcmsubset(nFaces::Int,grTopo::String,::Type{T},
  fSize::Array{NTuple{N, Int}},aSize::NTuple{N,Int},
  dims::NTuple{N,Int}) where {T,N}
  f=Array{Array{T,N},1}(undef,nFaces)
  i=Array{Array{T,N},1}(undef,nFaces)
  iSize=Array{NTuple{N, Int},1}(undef,nFaces)
  for a=1:nFaces
    f[a]=Array{T}(undef,fSize[a])
    #eventually I will distribute across faces; for now I just use face 1:
    a==1 ? nloc=dims[1] : nloc=0
    tmp1=Base.tail(dims)
    iSize[a]=(nloc,tmp1...)
    i[a]=Array{T}(undef,iSize[a])
  end
  gcmsubset{T,N}(nFaces,grTopo,f,fSize,aSize,i,iSize)
end

## Convenience functions

"""
    fijind(A::gcmfaces,ij::Int)

fijind convenience function.
"""
function fijind(A::gcmfaces,ij::Int)
  f=0
  j=0
  k=0
  tmp1=0
  for iFace=1:A.nFaces
    tmpsize=fsize(A,iFace)
    tmp11=tmpsize[1]*tmpsize[2]
    tmp2=tmp1+tmp11
    if tmp1<ij<=tmp2
      f=iFace;
      tmp3=(ij-tmp1);
      k=Int(ceil(tmp3/tmpsize[1]))
      j=Int(tmp3-tmpsize[1]*(k-1))
    end;
    tmp1=tmp1+tmp11
  end
  return (f,j,k)
end

"""
    fsize(A::AbstractGcmfaces{T, N}) where {T,N}

fsize convenience function.
"""
function fsize(A::AbstractGcmfaces{T, N}) where {T,N}
  fs=Array{NTuple{N, Int}}(undef,A.nFaces)
  for i=1:A.nFaces
    fs[i]=size(A.f[i]);
  end
  return fs
end

"""
    fsize(A::AbstractGcmfaces{T, N},i::Int) where {T,N}

fsize convenience function.
"""
function fsize(A::AbstractGcmfaces{T, N},i::Int) where {T,N}
  if i>0
    fs=size(A.f[i])
  else
    tmp1=0
    for i=1:A.nFaces
      tmp1=tmp1+size(A.f[i],1)*size(A.f[i],2)
    end
    tmp2=size(A.f[1])
    fs=(tmp1,1,tmp2[3:end]...)
  end
end

"""
    fsize(A::Array{Array{T,N},1}) where {T,N}

fsize convenience function.
"""
function fsize(A::Array{Array{T,N},1}) where {T,N}
  fs=Array{NTuple{N, Int}}(undef,length(A))
  for i=1:length(A)
    fs[i]=size(A[i]);
  end
  return fs
end

"""
    fsize(A::Array{Array{T,N},1},i::Int) where {T,N}

fsize convenience function.
"""
function fsize(A::Array{Array{T,N},1},i::Int) where {T,N}
  if i>0
    fs=size(A[i])
  else
    tmp1=0
    for i=1:length(A)
      tmp1=tmp1+size(A[i],1)*size(A[i],2)
    end
    tmp2=size(A[1])
    fs=(tmp1,1,tmp2[3:end]...)
  end
end

## Interface Methods

Base.size(A::gcmfaces) = fsize(A, 0)
Base.size(A::gcmfaces, dim::Integer) = fsize(A, 0)[dim]
Base.size(A::gcmsubset) = fsize(A.i, 0)
Base.size(A::gcmsubset, dim::Integer) = fsize(A.i, 0)[dim]

#

function Base.getindex(A::AbstractGcmfaces{T, N}, I::Vararg{Union{Int,AbstractUnitRange,Colon}, N}) where {T,N}
  if typeof(I[1])<:Int
    (f,i,j)=fijind(A,I[1])
    J=Base.tail(Base.tail(I))
    J=(i,j,J...)
    val=A.f[f][J...]
  elseif typeof(I[1])<:AbstractUnitRange
    val=similar(A,eltype(A),length.(I))
    for iFace=1:A.nFaces
      @views val.f[iFace]=A.f[iFace]
    end
    #eventually I will distribute across faces; for now I just use face 1:
    k=0
    J=Base.tail(Base.tail(I))
    for kk=I[1]
      k+=1
      (f,i,j)=fijind(A,kk)
      tmp1=(k,1,J...)
      tmp2=(i,j,J...)
      val.i[1][tmp1...]=A.f[f][tmp2...]
    end
  elseif typeof(I[1])<:Colon
    #should this rather be a copy as the above?
    val=view(A,I...)
  else
    er1=typeof(A)
    er2=typeof(I[1])
    error("getindex not yet implemented for "*"$er1"*" and "*"$er2"*" indices")
  end
  return val
end

function Base.getindex(A::gcmsubset{T, N}, I::Vararg{Union{Int,AbstractUnitRange,Colon}, N}) where {T,N}
  #eventually I will distribute across faces; for now I just use face 1:
  return getindex(A.i[1],I...)
end

#

function Base.setindex!(A::AbstractGcmfaces{T, N}, v, I::Vararg{Int, N}) where {T,N}
  (f,i,j)=fijind(A,I[1])
  J=Base.tail(Base.tail(I))
  J=(i,j,J...)
  return (A.f[f][J...] = v)
end

function Base.setindex!(A::gcmsubset{T, N}, v, I::Vararg{Int, N}) where {T,N}
  #eventually I will distribute across faces; for now I just use face 1:
  return (A.i[1][I...] = v)
end

## view

function Base.view(a::AbstractGcmfaces{T, N}, I::Vararg{Union{Int,AbstractUnitRange,Colon}, N}) where {T,N}
  nFaces=a.nFaces;
  grTopo=a.grTopo;
  if !isa(I[1],Colon)|!isa(I[2],Colon)
    J=Base.tail(Base.tail(I))
    J=(:,:,J...)
  else
    J=I
  end
  Nout=length(size(view(a.f[1],J...)));
  v1=Array{Array{T,Nout}}(undef,nFaces);
  for iFace=1:nFaces
    v1[iFace]=view(a.f[iFace],J...);
  end
  c=gcmfaces(nFaces,grTopo,v1);
  return c;
end

# Custom pretty-printing

function Base.show(io::IO, z::AbstractGcmfaces{T, N}) where {T,N}

#    @printf io " MeshArrays instance with \n"
    if isa(z,gcmfaces)
      printstyled(io, " gcmfaces array \n",color=:normal)
      nm="face"
      fs=fsize(z.f)
    elseif isa(z,gcmsubset)
      printstyled(io, " gcmsubset array \n",color=:normal)
      fs=fsize(z.i)
      nm="subset"
    else
      error("unknown type")
    end
    printstyled(io, "  grid type   = ",color=:normal)
    printstyled(io, "$(z.grTopo)\n",color=:blue)
    printstyled(io, "  # of faces  = ",color=:normal)
    printstyled(io, "$(z.nFaces)\n",color=:blue)
    if ~isassigned(z.f);
      printstyled(io, "  data type   = ",color=:normal)
      printstyled(io, "unassigned\n",color=:green)
      printstyled(io, "  face sizes  = ",color=:normal)
      printstyled(io, "unassigned\n",color=:green)
    else
      printstyled(io, "  data type   = ",color=:normal)
      printstyled(io, "$(typeof(z.f[1][1]))\n",color=:blue)
      printstyled(io, "  $(nm) sizes  = ",color=:normal)
      printstyled(io, "$(fs[1])\n",color=:blue)
      for iFace=2:z.nFaces
        printstyled(io, "                ",color=:normal)
        printstyled(io, "$(fs[iFace])\n",color=:blue)
      end
    end

    return
end

#

function Base.similar(A::gcmfaces, ::Type{T}, dims::Dims) where {T}
  if prod(dims)==length(A)
    B=gcmfaces(A.nFaces,A.grTopo,T,A.fSize,A.aSize)
  else
    B=gcmsubset(A.nFaces,A.grTopo,T,A.fSize,A.aSize,dims)
  end
end

Base.BroadcastStyle(::Type{<:gcmfaces}) = Broadcast.ArrayStyle{gcmfaces}()

function Base.similar(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{gcmfaces}}, ::Type{ElType}) where ElType
    # Scan the inputs for the gcmfaces:
    A = find_gcmfaces(bc)
    # Create the gcmfaces output:
    similar(A)
end

find_gcmfaces(bc::Base.Broadcast.Broadcasted) = find_gcmfaces(bc.args)
find_gcmfaces(args::Tuple) = find_gcmfaces(find_gcmfaces(args[1]), Base.tail(args))
find_gcmfaces(x) = x
find_gcmfaces(a::gcmfaces, rest) = a
find_gcmfaces(::Any, rest) = find_gcmfaces(rest)

#

function Base.similar(A::gcmsubset, ::Type{T}, dims::Dims) where {T}
    B=gcmsubset(A.nFaces,A.grTopo,T,A.fSize,A.aSize,dims[1])
end

Base.BroadcastStyle(::Type{<:gcmsubset}) = Broadcast.ArrayStyle{gcmsubset}()

function Base.similar(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{gcmsubset}}, ::Type{ElType}) where ElType
    # Scan the inputs for the gcmsubset:
    A = find_gcmsubset(bc)
    # Create the gcmsubset output:
    similar(A)
end

find_gcmsubset(bc::Base.Broadcast.Broadcasted) = find_gcmsubset(bc.args)
find_gcmsubset(args::Tuple) = find_gcmsubset(find_gcmsubset(args[1]), Base.tail(args))
find_gcmsubset(x) = x
find_gcmsubset(a::gcmsubset, rest) = a
find_gcmsubset(::Any, rest) = find_gcmsubset(rest)

# basic operations

import Base: +, -, *, /, getindex
import Base: isnan, isinf, isfinite
import Base: maximum, minimum, sum, fill

function +(a::gcmfaces)
  c=similar(a)
  for iFace=1:a.nFaces
    c.f[iFace]=+a.f[iFace];
  end
  return c
end

function +(a::gcmfaces,b::gcmfaces)
  c=similar(a)
  for iFace=1:a.nFaces
    c.f[iFace]=a.f[iFace]+b.f[iFace];
  end
  return c
  #the following modifies a:
  #  c=a
  #  c.f=a.f .+ b.f
  #the following fails immutability:
  #  c=gcmfaces(a.nFaces,a.grTopo)
  #  c.f=a.f .+ b.f
end

function +(a::Number,b::gcmfaces)
  c=similar(b)
  for iFace=1:b.nFaces
    c.f[iFace]=a.+b.f[iFace];
  end
  return c
  #the following is deprecated synthax as of v0.7:
  #  c=gcmfaces(b.nFaces,b.grTopo)
  #  c.f=a .+ b.f
end

function +(a::gcmfaces,b::Number)
  c=b+a
  return c
end

function -(a::gcmfaces)
  c=similar(a)
  for iFace=1:a.nFaces
    c.f[iFace]=-a.f[iFace];
  end
  return c
end

function -(a::gcmfaces,b::gcmfaces)
  c=similar(a)
  for iFace=1:a.nFaces
    c.f[iFace]=a.f[iFace]-b.f[iFace];
  end
  return c
end

function -(a::Number,b::gcmfaces)
  c=similar(b)
  for iFace=1:b.nFaces
    c.f[iFace]=a.-b.f[iFace];
  end
  return c
end

function -(a::gcmfaces,b::Number)
  c=similar(a)
  for iFace=1:a.nFaces
    c.f[iFace]=a.f[iFace].-b;
  end
  return c
end

function *(a::gcmfaces,b::gcmfaces)
  c=similar(a);
  for iFace=1:a.nFaces
    c.f[iFace]=a.f[iFace].*b.f[iFace];
  end
  return c
end

function *(a::Number,b::gcmfaces)
  c=similar(b);
  for iFace=1:b.nFaces
    c.f[iFace]=a*b.f[iFace];
  end
  return c
end

function *(a::gcmfaces,b::Number)
  c=b*a
  return c
end

function /(a::gcmfaces,b::gcmfaces)
  c=similar(a);
  for iFace=1:a.nFaces
    c.f[iFace]=a.f[iFace]./b.f[iFace];
  end
  return c
end

function /(a::Number,b::gcmfaces)
  c=similar(b);
  for iFace=1:b.nFaces
    c.f[iFace]=a./b.f[iFace];
  end
  return c
end

function /(a::gcmfaces,b::Number)
  c=similar(a);
  for iFace=1:a.nFaces
    c.f[iFace]=a.f[iFace]./b;
  end
  return c
end

##

function isnan(a::gcmfaces)
    c=similar(a);
    for iFace=1:a.nFaces
      c.f[iFace]=isnan.(a.f[iFace]);
    end
    return c
end

function isinf(a::gcmfaces)
    c=similar(a);
    for iFace=1:a.nFaces
      c.f[iFace]=isinf.(a.f[iFace]);
    end
    return c
end

function isfinite(a::gcmfaces)
    c=similar(a);
    for iFace=1:a.nFaces
      c.f[iFace]=isfinite.(a.f[iFace]);
    end
    return c
end

function sum(a::gcmfaces)
    c=0.0;
    for iFace=1:a.nFaces
      tmp1=a.f[iFace];
      c=c+sum(tmp1);
    end
    return c
end

function maximum(a::gcmfaces)
    c=-Inf;
    for iFace=1:a.nFaces
      tmp1=a.f[iFace];
      c=max(c,maximum(tmp1));
    end
    return c
end

function minimum(a::gcmfaces)
    c=Inf;
    for iFace=1:a.nFaces
      tmp1=a.f[iFace];
      c=min(c,minimum(tmp1));
    end
    return c
end

function fill(val::Any,a::gcmfaces)
    c=similar(a);
    for iFace=1:a.nFaces
      c.f[iFace]=fill(val,fsize(a,iFace));
    end
    return c
end

ext: jl


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

using LightGraphs
using SparseArrays

function readNetworkFiles(nodesFilename::String, arcsFilename::String)
	net = Network()
	fGraph = net.fGraph # shorthand
	fGraph.nodes = readNodesFile(nodesFilename)
	(fGraph.arcs, arcTravelTimes) = readArcsFile(arcsFilename)

	# create map
	mp = Map()
	mp.xMin = minimum(node -> node.location.x, fGraph.nodes) - 1
	mp.xMax = maximum(node -> node.location.x, fGraph.nodes) + 1
	mp.yMin = minimum(node -> node.location.y, fGraph.nodes) - 1
	mp.yMax = maximum(node -> node.location.y, fGraph.nodes) + 1

	JEMSS.initGraph!(fGraph)
	JEMSS.setArcDistances!(fGraph, mp)
	JEMSS.checkGraph(fGraph, mp)
	JEMSS.initFNetTravels!(net, arcTravelTimes)
	JEMSS.createRGraphFromFGraph!(net)
	JEMSS.checkGraph(net.rGraph, mp)
	JEMSS.createRNetTravelsFromFNetTravels!(net)
	JEMSS.setCommonFNodes!(net, Int[])

	return net
end

function fNetTravelShortestPathData(net::Network, travelModeIndex::Int)
	# Return spTimes and spSuccs for fNetTravel:
	# spTimes[i,j] gives the travel time from fNode i to fNode j along shortest path
	# spSuccs[i,j] gives index of successor node on shortest path from fNode i to fNode j

	# Most of this code is copied from the function JEMSS.calcRNetTravelShortestPaths!,
	# but changed to be for fNetTravel instead of rNetTravel.

	# shorthand
	fGraph = net.fGraph
	numNodes = length(fGraph.nodes)

	# calculate all-pairs shortest-paths data for the network with Dijkstra's algorithm
	revGraph = LightGraphs.reverse(fGraph.light) # reverse graph
	revArcTimes = spzeros(FloatSpTime, numNodes, numNodes)
	for arc in fGraph.arcs
		@assert(revArcTimes[arc.toNodeIndex, arc.fromNodeIndex] == 0)
		revArcTimes[arc.toNodeIndex, arc.fromNodeIndex] = net.fNetTravels[travelModeIndex].arcTimes[arc.index]
	end
	spTimes = Array{FloatSpTime,2}(undef, numNodes, numNodes) # spTimes[i,j] is the shortest-path travel time from node i to node j, to be calculated
	spSuccs = Array{IntRNode,2}(undef, numNodes, numNodes) # spSuccs[i,j] gives the successor node of i on shortest path from node i to node j, to be calculated
	for i = 1:numNodes
		# calculate shortest paths for origin node i
		spData = LightGraphs.dijkstra_shortest_paths(revGraph, i, revArcTimes) # type: LightGraphs.DijkstraState{Float}
		spTimes[:,i] = spData.dists
		spSuccs[:,i] = spData.parents
	end

	# change spSuccs where needed
	# where a path i -> l has spSuccs[i,l] = k, but spSuccs[i,k] = j,
	# need to change all cases where spSuccs[i,-] = k to instead have spSuccs[i,-] = j
	# this case can happen when an arc (i,k) has the same travel time as a path (i,j,k)
	for fArc in fGraph.arcs
		i = fArc.fromNodeIndex
		j = k = fArc.toNodeIndex
		while j != spSuccs[i,j]
			j = spSuccs[i,j]
		end
		if j != k
			# for any paths from i that have successor node k, change successor to be j
			spSuccs[i, spSuccs[i,:] .== k] .= j
		end
	end

	return spTimes, spSuccs
end

@testset "shortest path travel time" begin
	# Check shortestPathTravelTime function.
	# Compare against a simpler (but slower) function that calculates travel times for the whole full network,
	# unlike shortestPathTravelTime which takes advantage of the reduced network.

	testNetworksFolder = "data/network/small"
	for testFolderName in readdir(testNetworksFolder)
		folder = joinpath(testNetworksFolder, testFolderName)
		net = readNetworkFiles(joinpath(folder, "nodes.csv"), joinpath(folder, "arcs.csv"))
		for travelModeIndex = 1:length(net.fNetTravels)
			rtol = eps(FloatSpTime)
			atol = 2 * maximum(net.rNetTravels[travelModeIndex].arcTimes) * eps(FloatSpTime) # in NetTravel type, precision of fNodeToRNodeTime and fNodeFromRNodeTime (used in sp time calculation) are affected by rArc travel time. Multiplied by 2 since this can affect the travel time at both the start and end of the path.

			# compare travel time from shortestPathTravelTime function with spTimes for full network
			spTimes, spSuccs = fNetTravelShortestPathData(net, travelModeIndex)
			numNodes = length(net.fGraph.nodes) # shorthand
			allPass = true
			for i = 1:numNodes, j = 1:numNodes
				t1 = shortestPathTravelTime(net, travelModeIndex, i, j)
				t2 = spTimes[i,j]
				allPass &= abs(t1-t2) <= atol + rtol*max(t1,t2) # equivalent to isapprox(t1, t2, atol = atol, rtol = rtol) for julia v0.6
			end
			@test allPass
		end
	end
end

@testset "shortest path distance" begin
	# Check shortestPathDistance function.
	# Compare against a simpler (but slower) function that calculates distances based on arc traversed in full network,
	# unlike shortestPathDistance which takes advantage of the reduced network.

	# Calculate distance along shortest path between two nodes, using fGraph.
	# spSuccs[i,j] gives successor of node i in shortest path from node i to node j, see fNetTravelShortestPathData.
	function shortestPathDistance2(net::Network, spSuccs::Array{IntRNode,2}, startFNode::Int, endFNode::Int)::Float
		dist = 0.0
		fNode = startFNode
		while fNode != endFNode
			nextFNode = spSuccs[fNode, endFNode]
			arcIndex = net.fGraph.nodePairArcIndex[fNode, nextFNode]
			dist += net.fGraph.arcs[arcIndex].distance
			fNode = nextFNode
		end
		return dist
	end

	testNetworksFolder = "data/network/small"
	for testFolderName in readdir(testNetworksFolder)
		folder = joinpath(testNetworksFolder, testFolderName)
		net = readNetworkFiles(joinpath(folder, "nodes.csv"), joinpath(folder, "arcs.csv"))
		for travelModeIndex = 1:length(net.fNetTravels)
			# compare shortestPathDistance and shortestPathDistance2 for each pair of nodes in fGraph
			spTimes, spSuccs = fNetTravelShortestPathData(net, travelModeIndex)
			numNodes = length(net.fGraph.nodes) # shorthand
			allPass = true
			for i = 1:numNodes, j = 1:numNodes
				d1 = shortestPathDistance(net, travelModeIndex, i, j) # from JEMSS
				d2 = shortestPathDistance2(net, spSuccs, i, j)
				allPass &= isapprox(d1, d2, rtol = eps())
			end
			@test allPass
		end
	end
end


text:
@doc "m_1" ->
macro m_1(x)
end

@doc """m_2""" ->
macro m_2(x)
end


text:
using EllZeroTrendFiltering, Test

include(joinpath(dirname(@__FILE__),"auxilliary_test_fcns.jl"))

Random.seed!(1234)

# Brute force up to nbr segments:
M = 7
# Total number of points:
N = 40
# Number of internal gridpoints:
n = 12

#Function
g = cumsum(randn(N).^2)

t_grids = Vector{Int}[1:N, [[1;sort!(shuffle(2:N-1)[1:n]);N] for i = 1:10]...]
@testset "Testing subgrid, regularized, with endpoints, i=$i" for (i,t_grid) in enumerate(t_grids)
    local I_sols  # TODO, can remove in julia 1.0?
    I_sol, Y_sol, f_sol = brute_force_multi(g, M, t_grid, tol=1e-3)

    I, Y, f = fit_pwl_constrained(g, M, t=t_grid)
    I_t = getindex.([t_grid], I)
    @test I_t == I_sol
    @test Y ≈ Y_sol atol = 1e-10 # Crazy that this is possible
    @test f ≈ f_sol atol = 1e-10

    ζ_vec = 10 .^ range(log10(f_sol[2]), stop=log10(0.8*(max(f_sol[end-1] - f_sol[end], f_sol[end-1]))), length=3)

    for ζ in ζ_vec
        local I2, Y2  # TODO, can remove in julia 1.0?
        I2, Y2, f2 = fit_pwl_regularized(g, ζ, t=t_grid)
        I2_t = t_grid[I2]
        #Make sure regularized solution is in constrained solution
        idx = findall([I2_t] .== I_t)
        @test length(idx) == 1
        if length(idx) == 1
            @test Y2 ≈ Y_sol[idx[1]] atol = 1e-10
            @test f2 ≈ f_sol[idx[1]] atol = 1e-10
        end
    end
end

println("Testing sub-grid that is not covering the full range of input, expect warnings.")
# Test when endpoints are not included
t_grids = Vector{Int}[1:N, [sort!(shuffle(2:N-1)[1:n]) for i = 1:10]...]
@testset "Testing subgrid, constrained, without endpoints, i=$i" for (i, t_grid) in enumerate(t_grids)
    I_sol, Y_sol, f_sol = brute_force_multi(g, M-2, t_grid, tol=1e-3)

    I, Y, f = fit_pwl_constrained(g, M-2, t=t_grid)
    I_t = getindex.([t_grid], I)
    @test I_t == I_sol
    @test Y ≈ Y_sol atol = 1e-10 # Crazy that this is possible
    @test f ≈ f_sol atol = 1e-10

    ζ_vec = 10 .^ range(log10(f_sol[2]), stop=log10(0.8*(max(f_sol[end-1] - f_sol[end], f_sol[end-1]))), length=3)

    for ζ in ζ_vec
        local I2, Y2  # TODO, can remove in julia 1.0?
        I2, Y2, f2 = fit_pwl_regularized(g, ζ, t=t_grid)
        I2_t = t_grid[I2]
        #Make sure regularized solution is in constrained solution
        idx = findall([I2_t] .== I_t)
        @test length(idx) == 1
        if length(idx) == 1
            @test Y2 ≈ Y_sol[idx[1]] atol = 1e-10
            @test f2 ≈ f_sol[idx[1]] atol = 1e-10
        end
    end
end


text:
#==============================================================================#
#=                              Simplified LCMC                               =#
#==============================================================================#

const LCMC_CATEGORIES = Dict(
    'A' => "Press: reportage",
    'B' => "Press: editorials",
    'C' => "Press: reviews",
    'D' => "Religion",
    'E' => "Skills, trades and hobbies",
    'F' => "Popular lore",
    'G' => "Biographies and essays",
    'H' => "Miscellaneous: reports and official documents",
    'J' => "Science: academic prose",
    'K' => "General fiction",
    'L' => "Mystery and detective fiction",
    'M' => "Science fiction",
    'N' => "Adventure and martial arts fiction",
    'P' => "Romantic fiction",
    'R' => "Humour")

"""
    SimplifiedLCMC([categories])

A character frequency dataset: Lancaster Corpus for Mandarin Chinese, simplified terms only,
based on simplified text corpus. See their
[website](https://www.lancaster.ac.uk/fass/projects/corpus/LCMC/default.htm) for more details about the corpus.

The character frequency can be based only on selected categories (see `CJKFrequencies.LCMC_CATEGORIES` for valid
 category keys and corresponding category names). Any invalid categories will be ignored.

## Examples
Loading all the categories:
```julia-repl
julia> charfreq(SimplifiedLCMC())
DataStructures.Accumulator{String,Int64} with 45411 entries:
  "一路…   => 1
  "舍得"   => 9
  "５８"   => 1
  "神农…   => 1
  "十点"   => 8
  "随从"   => 9
  "荡心…   => 1
  "尺码"   => 1
  ⋮      => ⋮
```
Or loading just a subset (argument can be any iterable):
```julia-repl
julia> charfreq(SimplifiedLCMC("ABEGKLMNR"))
DataStructures.Accumulator{String,Int64} with 35488 entries:
  "废…  => 1
  "蜷"  => 1
  "哇"  => 13
  "丰…  => 1
  "弊…  => 3
  "议…  => 10
  "滴"  => 28
  "美…  => 1
  ⋮    => ⋮
```

## Licensing/Copyright
Note: This corpus has some conflicting licensing information, depending on who is supplying the
data.

The original corpus is provided primarily for non-profit-making research. Be sure to see the full
[end user license agreement](https://www.lancaster.ac.uk/fass/projects/corpus/LCMC/lcmc/lcmc_license.htm).

Via the
[Oxford Text Archive](https://ota.bodleian.ox.ac.uk/repository/xmlui/handle/20.500.12024/2474),
this corpus is distributed under the
[CC BY-NC-SA 3.0](http://creativecommons.org/licenses/by-nc-sa/3.0/) license.
"""
struct SimplifiedLCMC
    categories::Set{Char}
    function SimplifiedLCMC(cats)
        lcmc = new(Set{String}())
        for cat in cats
            haskey(LCMC_CATEGORIES, cat) && push!(lcmc.categories, cat)
        end
        lcmc
    end
    SimplifiedLCMC() = new(keys(LCMC_CATEGORIES))
end

function charfreq(lcmc::SimplifiedLCMC)
    cf = CJKFrequency()
    for cat in lcmc.categories
        filename = joinpath(artifact"lcmc", "LCMC_$(cat).XML")
        doc = parse_file(filename)
        _words_from_xml(root(doc), cf)
    end
    cf
end

function _words_from_xml(xml_elem, accum)
    for c in child_nodes(xml_elem)
        if name(c) == "w"
            inc!(accum, content(c))
        else
            _words_from_xml(c, accum)
        end
    end
end

#==============================================================================#
#=                             Simplified Jun Da                              =#
#==============================================================================#
"""
    SimplifiedJunDa()

A character frequency
[dataset](https://lingua.mtsu.edu/chinese-computing/)
 of modern Chinese compiled by Jun Da, simplified single-character
words only.

Currently, only the modern Chinese dataset is fetched; however, in the future, the other lists may
also be provided as an option.

## Examples
```julia-repl
julia> charfreq(SimplifiedJunDa())
DataStructures.Accumulator{String,Int64} with 9932 entries:
  "蜷… => 837
  "哇… => 4055
  "湓… => 62
  "滴… => 8104
  "堞… => 74
  "狭… => 6901
  "尚… => 38376
  "懈… => 2893
  ⋮   => ⋮
```

## Licensing/Copyright
The original author maintains full copyright to the character frequency lists, but provides the
lists for research and teaching/learning purposes only, no commercial use without permission from
 the author. See their full disclaimer and copyright notice [here](https://lingua.mtsu.edu/chinese-computing/copyright.html).
"""
struct SimplifiedJunDa end

function charfreq(::SimplifiedJunDa)
    cf = CJKFrequency()
    pattern = r"^\d+\s+(\w)\s+(\d+)\s+\d+(?:\.\d+)\s+.+$"
    for line in eachline(joinpath(artifact"junda", "freq.txt"))
        m = match(pattern, line)
        m !== nothing && inc!(cf, m.captures[1], Base.parse(Int, m.captures[2]))
    end
    cf
end


text: Dict{String,Any}("mixed" => Dict{String,Any}("value" => Any[Dict{String,Any}("value" => Any[Dict{String,Any}("value" => "1","type" => "integer"), Dict{String,Any}("value" => "2","type" => "integer")],"type" => "array"), Dict{String,Any}("value" => Any[Dict{String,Any}("value" => "a","type" => "string"), Dict{String,Any}("value" => "b","type" => "string")],"type" => "array"), Dict{String,Any}("value" => Any[Dict{String,Any}("value" => "1.1","type" => "float"), Dict{String,Any}("value" => "2.1","type" => "float")],"type" => "array")],"type" => "array"))

max_stars_repo_name: SmiteWindows/julia

max_issues_repo_name: SmiteWindows/julia

max_forks_repo_name: SmiteWindows/julia


text:
int_rules_6_7_7 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 6.7.7*F^(c*(a+b*x))*hyper(d+e*x)^n =#
    @apply_utils Antiderivative((~F) ^ (~(c') * (~(a') + ~(b') * ~x)) * sinh(~(d') + ~(e') * ~x), ~x) => (-(~b) * ~c * log(~F) * (~F) ^ (~c * (~a + ~b * ~x)) * sinh(~d + ~e * ~x)) / ((~e) ^ 2 - (~b) ^ 2 * (~c) ^ 2 * log(~F) ^ 2) + (~e * (~F) ^ (~c * (~a + ~b * ~x)) * cosh(~d + ~e * ~x)) / ((~e) ^ 2 - (~b) ^ 2 * (~c) ^ 2 * log(~F) ^ 2) <-- FreeQ([~F, ~a, ~b, ~c, ~d, ~e], ~x) && NeQ((~e) ^ 2 - (~b) ^ 2 * (~c) ^ 2 * log(~F) ^ 2, 0)
    @apply_utils Antiderivative((~F) ^ (~(c') * (~(a') + ~(b') * ~x)) * cosh(~(d') + ~(e') * ~x), ~x) => (-(~b) * ~c * log(~F) * (~F) ^ (~c * (~a + ~b * ~x)) * cosh(~d + ~e * ~x)) / ((~e) ^ 2 - (~b) ^ 2 * (~c) ^ 2 * log(~F) ^ 2) + (~e * (~F) ^ (~c * (~a + ~b * ~x)) * sinh(~d + ~e * ~x)) / ((~e) ^ 2 - (~b) ^ 2 * (~c) ^ 2 * log(~F) ^ 2) <-- FreeQ([~F, ~a, ~b, ~c, ~d, ~e], ~x) && NeQ((~e) ^ 2 - (~b) ^ 2 * (~c) ^ 2 * log(~F) ^ 2, 0)
    @apply_utils Antiderivative((~F) ^ (~(c') * (~(a') + ~(b') * ~x)) * sinh(~(d') + ~(e') * ~x) ^ ~n, ~x) => ((-(~b) * ~c * log(~F) * (~F) ^ (~c * (~a + ~b * ~x)) * sinh(~d + ~e * ~x) ^ ~n) / ((~e) ^ 2 * (~n) ^ 2 - (~b) ^ 2 * (~c) ^ 2 * log(~F) ^ 2) + (~e * ~n * (~F) ^ (~c * (~a + ~b * ~x)) * cosh(~d + ~e * ~x) * sinh(~d + ~e * ~x) ^ (~n - 1)) / ((~e) ^ 2 * (~n) ^ 2 - (~b) ^ 2 * (~c) ^ 2 * log(~F) ^ 2)) - ((~n * (~n - 1) * (~e) ^ 2) / ((~e) ^ 2 * (~n) ^ 2 - (~b) ^ 2 * (~c) ^ 2 * log(~F) ^ 2)) * Antiderivative((~F) ^ (~c * (~a + ~b * ~x)) * sinh(~d + ~e * ~x) ^ (~n - 2), ~x) <-- FreeQ([~F, ~a, ~b, ~c, ~d, ~e], ~x) && (NeQ((~e) ^ 2 * (~n) ^ 2 - (~b) ^ 2 * (~c) ^ 2 * log(~F) ^ 2, 0) && GtQ(~n, 1))
    @apply_utils Antiderivative((~F) ^ (~(c') * (~(a') + ~(b') * ~x)) * cosh(~(d') + ~(e') * ~x) ^ ~n, ~x) => (-(~b) * ~c * log(~F) * (~F) ^ (~c * (~a + ~b * ~x)) * cosh(~d + ~e * ~x) ^ ~n) / ((~e) ^ 2 * (~n) ^ 2 - (~b) ^ 2 * (~c) ^ 2 * log(~F) ^ 2) + (~e * ~n * (~F) ^ (~c * (~a + ~b * ~x)) * sinh(~d + ~e * ~x) * cosh(~d + ~e * ~x) ^ (~n - 1)) / ((~e) ^ 2 * (~n) ^ 2 - (~b) ^ 2 * (~c) ^ 2 * log(~F) ^ 2) + ((~n * (~n - 1) * (~e) ^ 2) / ((~e) ^ 2 * (~n) ^ 2 - (~b) ^ 2 * (~c) ^ 2 * log(~F) ^ 2)) * Antiderivative((~F) ^ (~c * (~a + ~b * ~x)) * cosh(~d + ~e * ~x) ^ (~n - 2), ~x) <-- FreeQ([~F, ~a, ~b, ~c, ~d, ~e], ~x) && (NeQ((~e) ^ 2 * (~n) ^ 2 - (~b) ^ 2 * (~c) ^ 2 * log(~F) ^ 2, 0) && GtQ(~n, 1))
    @apply_utils Antiderivative((~F) ^ (~(c') * (~(a') + ~(b') * ~x)) * sinh(~(d') + ~(e') * ~x) ^ ~n, ~x) => (-(~b) * ~c * log(~F) * (~F) ^ (~c * (~a + ~b * ~x)) * sinh(~d + ~e * ~x) ^ (~n + 2)) / ((~e) ^ 2 * (~n + 1) * (~n + 2)) + ((~F) ^ (~c * (~a + ~b * ~x)) * cosh(~d + ~e * ~x) * sinh(~d + ~e * ~x) ^ (~n + 1)) / (~e * (~n + 1)) <-- FreeQ([~F, ~a, ~b, ~c, ~d, ~e, ~n], ~x) && (EqQ((~e) ^ 2 * (~n + 2) ^ 2 - (~b) ^ 2 * (~c) ^ 2 * log(~F) ^ 2, 0) && (NeQ(~n, -1) && NeQ(~n, -2)))
    @apply_utils Antiderivative((~F) ^ (~(c') * (~(a') + ~(b') * ~x)) * cosh(~(d') + ~(e') * ~x) ^ ~n, ~x) => (~b * ~c * log(~F) * (~F) ^ (~c * (~a + ~b * ~x)) * cosh(~d + ~e * ~x) ^ (~n + 2)) / ((~e) ^ 2 * (~n + 1) * (~n + 2)) - ((~F) ^ (~c * (~a + ~b * ~x)) * sinh(~d + ~e * ~x) * cosh(~d + ~e * ~x) ^ (~n + 1)) / (~e * (~n + 1)) <-- FreeQ([~F, ~a, ~b, ~c, ~d, ~e, ~n], ~x) && (EqQ((~e) ^ 2 * (~n + 2) ^ 2 - (~b) ^ 2 * (~c) ^ 2 * log(~F) ^ 2, 0) && (NeQ(~n, -1) && NeQ(~n, -2)))
    @apply_utils Antiderivative((~F) ^ (~(c') * (~(a') + ~(b') * ~x)) * sinh(~(d') + ~(e') * ~x) ^ ~n, ~x) => ((-(~b) * ~c * log(~F) * (~F) ^ (~c * (~a + ~b * ~x)) * sinh(~d + ~e * ~x) ^ (~n + 2)) / ((~e) ^ 2 * (~n + 1) * (~n + 2)) + ((~F) ^ (~c * (~a + ~b * ~x)) * cosh(~d + ~e * ~x) * sinh(~d + ~e * ~x) ^ (~n + 1)) / (~e * (~n + 1))) - (((~e) ^ 2 * (~n + 2) ^ 2 - (~b) ^ 2 * (~c) ^ 2 * log(~F) ^ 2) / ((~e) ^ 2 * (~n + 1) * (~n + 2))) * Antiderivative((~F) ^ (~c * (~a + ~b * ~x)) * sinh(~d + ~e * ~x) ^ (~n + 2), ~x) <-- FreeQ([~F, ~a, ~b, ~c, ~d, ~e], ~x) && (NeQ((~e) ^ 2 * (~n + 2) ^ 2 - (~b) ^ 2 * (~c) ^ 2 * log(~F) ^ 2, 0) && (LtQ(~n, -1) && NeQ(~n, -2)))
    @apply_utils Antiderivative((~F) ^ (~(c') * (~(a') + ~(b') * ~x)) * cosh(~(d') + ~(e') * ~x) ^ ~n, ~x) => ((~b * ~c * log(~F) * (~F) ^ (~c * (~a + ~b * ~x)) * cosh(~d + ~e * ~x) ^ (~n + 2)) / ((~e) ^ 2 * (~n + 1) * (~n + 2)) - ((~F) ^ (~c * (~a + ~b * ~x)) * sinh(~d + ~e * ~x) * cosh(~d + ~e * ~x) ^ (~n + 1)) / (~e * (~n + 1))) + (((~e) ^ 2 * (~n + 2) ^ 2 - (~b) ^ 2 * (~c) ^ 2 * log(~F) ^ 2) / ((~e) ^ 2 * (~n + 1) * (~n + 2))) * Antiderivative((~F) ^ (~c * (~a + ~b * ~x)) * cosh(~d + ~e * ~x) ^ (~n + 2), ~x) <-- FreeQ([~F, ~a, ~b, ~c, ~d, ~e], ~x) && (NeQ((~e) ^ 2 * (~n + 2) ^ 2 - (~b) ^ 2 * (~c) ^ 2 * log(~F) ^ 2, 0) && (LtQ(~n, -1) && NeQ(~n, -2)))
    @apply_utils Antiderivative((~F) ^ (~(c') * (~(a') + ~(b') * ~x)) * sinh(~(d') + ~(e') * ~x) ^ ~n, ~x) => ((ℯ ^ (~n * (~d + ~e * ~x)) * sinh(~d + ~e * ~x) ^ ~n) / (-1 + ℯ ^ (2 * (~d + ~e * ~x))) ^ ~n) * Antiderivative(((~F) ^ (~c * (~a + ~b * ~x)) * (-1 + ℯ ^ (2 * (~d + ~e * ~x))) ^ ~n) / ℯ ^ (~n * (~d + ~e * ~x)), ~x) <-- FreeQ([~F, ~a, ~b, ~c, ~d, ~e, ~n], ~x) && Not(IntegerQ(~n))
    @apply_utils Antiderivative((~F) ^ (~(c') * (~(a') + ~(b') * ~x)) * cosh(~(d') + ~(e') * ~x) ^ ~n, ~x) => ((ℯ ^ (~n * (~d + ~e * ~x)) * cosh(~d + ~e * ~x) ^ ~n) / (1 + ℯ ^ (2 * (~d + ~e * ~x))) ^ ~n) * Antiderivative(((~F) ^ (~c * (~a + ~b * ~x)) * (1 + ℯ ^ (2 * (~d + ~e * ~x))) ^ ~n) / ℯ ^ (~n * (~d + ~e * ~x)), ~x) <-- FreeQ([~F, ~a, ~b, ~c, ~d, ~e, ~n], ~x) && Not(IntegerQ(~n))
    @apply_utils Antiderivative((~F) ^ (~(c') * (~(a') + ~(b') * ~x)) * tanh(~(d') + ~(e') * ~x) ^ ~(n'), ~x) => Antiderivative(ExpandIntegrand(((~F) ^ (~c * (~a + ~b * ~x)) * (-1 + ℯ ^ (2 * (~d + ~e * ~x))) ^ ~n) / (1 + ℯ ^ (2 * (~d + ~e * ~x))) ^ ~n, ~x), ~x) <-- FreeQ([~F, ~a, ~b, ~c, ~d, ~e], ~x) && IntegerQ(~n)
    @apply_utils Antiderivative((~F) ^ (~(c') * (~(a') + ~(b') * ~x)) * coth(~(d') + ~(e') * ~x) ^ ~(n'), ~x) => Antiderivative(ExpandIntegrand(((~F) ^ (~c * (~a + ~b * ~x)) * (1 + ℯ ^ (2 * (~d + ~e * ~x))) ^ ~n) / (-1 + ℯ ^ (2 * (~d + ~e * ~x))) ^ ~n, ~x), ~x) <-- FreeQ([~F, ~a, ~b, ~c, ~d, ~e], ~x) && IntegerQ(~n)
    @apply_utils Antiderivative((~F) ^ (~(c') * (~(a') + ~(b') * ~x)) * sech(~(d') + ~(e') * ~x) ^ ~n, ~x) => (-(~b) * ~c * log(~F) * (~F) ^ (~c * (~a + ~b * ~x)) * (sech(~d + ~e * ~x) ^ ~n / ((~e) ^ 2 * (~n) ^ 2 - (~b) ^ 2 * (~c) ^ 2 * log(~F) ^ 2)) - ~e * ~n * (~F) ^ (~c * (~a + ~b * ~x)) * sech(~d + ~e * ~x) ^ (~n + 1) * (sinh(~d + ~e * ~x) / ((~e) ^ 2 * (~n) ^ 2 - (~b) ^ 2 * (~c) ^ 2 * log(~F) ^ 2))) + (~e) ^ 2 * ~n * ((~n + 1) / ((~e) ^ 2 * (~n) ^ 2 - (~b) ^ 2 * (~c) ^ 2 * log(~F) ^ 2)) * Antiderivative((~F) ^ (~c * (~a + ~b * ~x)) * sech(~d + ~e * ~x) ^ (~n + 2), ~x) <-- FreeQ([~F, ~a, ~b, ~c, ~d, ~e], ~x) && (NeQ((~e) ^ 2 * (~n) ^ 2 + (~b) ^ 2 * (~c) ^ 2 * log(~F) ^ 2, 0) && LtQ(~n, -1))
    @apply_utils Antiderivative((~F) ^ (~(c') * (~(a') + ~(b') * ~x)) * csch(~(d') + ~(e') * ~x) ^ ~n, ~x) => (-(~b) * ~c * log(~F) * (~F) ^ (~c * (~a + ~b * ~x)) * (csch(~d + ~e * ~x) ^ ~n / ((~e) ^ 2 * (~n) ^ 2 - (~b) ^ 2 * (~c) ^ 2 * log(~F) ^ 2)) - ~e * ~n * (~F) ^ (~c * (~a + ~b * ~x)) * csch(~d + ~e * ~x) ^ (~n + 1) * (cosh(~d + ~e * ~x) / ((~e) ^ 2 * (~n) ^ 2 - (~b) ^ 2 * (~c) ^ 2 * log(~F) ^ 2))) - (~e) ^ 2 * ~n * ((~n + 1) / ((~e) ^ 2 * (~n) ^ 2 - (~b) ^ 2 * (~c) ^ 2 * log(~F) ^ 2)) * Antiderivative((~F) ^ (~c * (~a + ~b * ~x)) * csch(~d + ~e * ~x) ^ (~n + 2), ~x) <-- FreeQ([~F, ~a, ~b, ~c, ~d, ~e], ~x) && (NeQ((~e) ^ 2 * (~n) ^ 2 + (~b) ^ 2 * (~c) ^ 2 * log(~F) ^ 2, 0) && LtQ(~n, -1))
    @apply_utils Antiderivative((~F) ^ (~(c') * (~(a') + ~(b') * ~x)) * sech(~(d') + ~(e') * ~x) ^ ~n, ~x) => (~b * ~c * log(~F) * (~F) ^ (~c * (~a + ~b * ~x)) * sech(~d + ~e * ~x) ^ (~n - 2)) / ((~e) ^ 2 * (~n - 1) * (~n - 2)) + ((~F) ^ (~c * (~a + ~b * ~x)) * sech(~d + ~e * ~x) ^ (~n - 1) * sinh(~d + ~e * ~x)) / (~e * (~n - 1)) <-- FreeQ([~F, ~a, ~b, ~c, ~d, ~e, ~n], ~x) && (EqQ((~e) ^ 2 * (~n - 2) ^ 2 - (~b) ^ 2 * (~c) ^ 2 * log(~F) ^ 2, 0) && (NeQ(~n, 1) && NeQ(~n, 2)))
    @apply_utils Antiderivative((~F) ^ (~(c') * (~(a') + ~(b') * ~x)) * csch(~(d') + ~(e') * ~x) ^ ~n, ~x) => (-(~b) * ~c * log(~F) * (~F) ^ (~c * (~a + ~b * ~x)) * csch(~d + ~e * ~x) ^ (~n - 2)) / ((~e) ^ 2 * (~n - 1) * (~n - 2)) - ((~F) ^ (~c * (~a + ~b * ~x)) * csch(~d + ~e * ~x) ^ (~n - 1) * cosh(~d + ~e * ~x)) / (~e * (~n - 1)) <-- FreeQ([~F, ~a, ~b, ~c, ~d, ~e, ~n], ~x) && (EqQ((~e) ^ 2 * (~n - 2) ^ 2 - (~b) ^ 2 * (~c) ^ 2 * log(~F) ^ 2, 0) && (NeQ(~n, 1) && NeQ(~n, 2)))
    @apply_utils Antiderivative((~F) ^ (~(c') * (~(a') + ~(b') * ~x)) * sech(~(d') + ~(e') * ~x) ^ ~n, ~x) => (~b * ~c * log(~F) * (~F) ^ (~c * (~a + ~b * ~x)) * sech(~d + ~e * ~x) ^ (~n - 2)) / ((~e) ^ 2 * (~n - 1) * (~n - 2)) + ((~F) ^ (~c * (~a + ~b * ~x)) * sech(~d + ~e * ~x) ^ (~n - 1) * sinh(~d + ~e * ~x)) / (~e * (~n - 1)) + (((~e) ^ 2 * (~n - 2) ^ 2 - (~b) ^ 2 * (~c) ^ 2 * log(~F) ^ 2) / ((~e) ^ 2 * (~n - 1) * (~n - 2))) * Antiderivative((~F) ^ (~c * (~a + ~b * ~x)) * sech(~d + ~e * ~x) ^ (~n - 2), ~x) <-- FreeQ([~F, ~a, ~b, ~c, ~d, ~e], ~x) && (NeQ((~e) ^ 2 * (~n - 2) ^ 2 - (~b) ^ 2 * (~c) ^ 2 * log(~F) ^ 2, 0) && (GtQ(~n, 1) && NeQ(~n, 2)))
    @apply_utils Antiderivative((~F) ^ (~(c') * (~(a') + ~(b') * ~x)) * csch(~(d') + ~(e') * ~x) ^ ~n, ~x) => ((-(~b) * ~c * log(~F) * (~F) ^ (~c * (~a + ~b * ~x)) * csch(~d + ~e * ~x) ^ (~n - 2)) / ((~e) ^ 2 * (~n - 1) * (~n - 2)) - ((~F) ^ (~c * (~a + ~b * ~x)) * csch(~d + ~e * ~x) ^ (~n - 1) * cosh(~d + ~e * ~x)) / (~e * (~n - 1))) - (((~e) ^ 2 * (~n - 2) ^ 2 - (~b) ^ 2 * (~c) ^ 2 * log(~F) ^ 2) / ((~e) ^ 2 * (~n - 1) * (~n - 2))) * Antiderivative((~F) ^ (~c * (~a + ~b * ~x)) * csch(~d + ~e * ~x) ^ (~n - 2), ~x) <-- FreeQ([~F, ~a, ~b, ~c, ~d, ~e], ~x) && (NeQ((~e) ^ 2 * (~n - 2) ^ 2 - (~b) ^ 2 * (~c) ^ 2 * log(~F) ^ 2, 0) && (GtQ(~n, 1) && NeQ(~n, 2)))
    #= Antiderivative((~F)^((~c')*((~a')+(~b')*(~x)))*sech((~d')+(~e')*(~x))^(~n'),~x) :=  2^n*Antiderivative(SimplifyIntegrand(F^(c*(a+b*x))*ℯ^(n*(d+e*x))/(1+ℯ^(2*(d+e*x)) )^n,x),x) <-- FreeQ([F,a,b,c,d,e],x) && IntegerQ(n) =#
    #= Antiderivative((~F)^((~c')*((~a')+(~b')*(~x)))*csch((~d')+(~e')*(~x))^(~n'),~x) :=  2^n*Antiderivative(SimplifyIntegrand(F^(c*(a+b*x))*ℯ^(-n*(d+e*x))/(1-ℯ^(-2*(d+e* x)))^n,x),x) <-- FreeQ([F,a,b,c,d,e],x) && IntegerQ(n) =#
    # Missing due to parse error: Antiderivative((~F)^((~c')*((~a') + (~b')*(~x)))*sech((~d') + (~e')*(~x))^(~n'), ~x) := 2^n*ℯ^(n*(d + e*x))*F^(c*(a + b*x))/(e*n + b*c*log(F))* HypergeometricFunctions._₂F₁(n, n/2 + b*c*log(F)/(2*e), 1 + n/2 + b*c*log(F)/(2*e), -ℯ^(2*(d + e*x))) <-- FreeQ([F, a, b, c, d, e], x) && IntegerQ(n)
    # Missing due to parse error: Antiderivative((~F)^((~c')*((~a') + (~b')*(~x)))*csch((~d') + (~e')*(~x))^(~n'), ~x) := (-2)^n*ℯ^(n*(d + e*x))* F^(c*(a + b*x))/(e*n + b*c*log(F))* HypergeometricFunctions._₂F₁(n, n/2 + b*c*log(F)/(2*e), 1 + n/2 + b*c*log(F)/(2*e), ℯ^(2*(d + e*x))) <-- FreeQ([F, a, b, c, d, e], x) && IntegerQ(n)
    @apply_utils Antiderivative((~F) ^ (~(c') * (~(a') + ~(b') * ~x)) * sech(~(d') + ~(e') * ~x) ^ ~(n'), ~x) => (((1 + ℯ ^ (2 * (~d + ~e * ~x))) ^ ~n * sech(~d + ~e * ~x) ^ ~n) / ℯ ^ (~n * (~d + ~e * ~x))) * Antiderivative(SimplifyIntegrand(((~F) ^ (~c * (~a + ~b * ~x)) * ℯ ^ (~n * (~d + ~e * ~x))) / (1 + ℯ ^ (2 * (~d + ~e * ~x))) ^ ~n, ~x), ~x) <-- FreeQ([~F, ~a, ~b, ~c, ~d, ~e], ~x) && Not(IntegerQ(~n))
    @apply_utils Antiderivative((~F) ^ (~(c') * (~(a') + ~(b') * ~x)) * csch(~(d') + ~(e') * ~x) ^ ~(n'), ~x) => (((1 - ℯ ^ (-2 * (~d + ~e * ~x))) ^ ~n * csch(~d + ~e * ~x) ^ ~n) / ℯ ^ (-(~n) * (~d + ~e * ~x))) * Antiderivative(SimplifyIntegrand(((~F) ^ (~c * (~a + ~b * ~x)) * ℯ ^ (-(~n) * (~d + ~e * ~x))) / (1 - ℯ ^ (-2 * (~d + ~e * ~x))) ^ ~n, ~x), ~x) <-- FreeQ([~F, ~a, ~b, ~c, ~d, ~e], ~x) && Not(IntegerQ(~n))
    @apply_utils Antiderivative((~F) ^ (~(c') * (~(a') + ~(b') * ~x)) * (~f + ~(g') * sinh(~(d') + ~(e') * ~x)) ^ ~(n'), ~x) => 2 ^ ~n * (~f) ^ ~n * Antiderivative((~F) ^ (~c * (~a + ~b * ~x)) * cosh((~d / 2 - (~f * pi) / (4 * ~g)) + (~e * ~x) / 2) ^ (2 * ~n), ~x) <-- FreeQ([~F, ~a, ~b, ~c, ~d, ~e, ~f, ~g], ~x) && (EqQ((~f) ^ 2 + (~g) ^ 2, 0) && ILtQ(~n, 0))
    @apply_utils Antiderivative((~F) ^ (~(c') * (~(a') + ~(b') * ~x)) * (~f + ~(g') * cosh(~(d') + ~(e') * ~x)) ^ ~(n'), ~x) => 2 ^ ~n * (~g) ^ ~n * Antiderivative((~F) ^ (~c * (~a + ~b * ~x)) * cosh(~d / 2 + (~e * ~x) / 2) ^ (2 * ~n), ~x) <-- FreeQ([~F, ~a, ~b, ~c, ~d, ~e, ~f, ~g], ~x) && (EqQ(~f - ~g, 0) && ILtQ(~n, 0))
    @apply_utils Antiderivative((~F) ^ (~(c') * (~(a') + ~(b') * ~x)) * (~f + ~(g') * cosh(~(d') + ~(e') * ~x)) ^ ~(n'), ~x) => 2 ^ ~n * (~g) ^ ~n * Antiderivative((~F) ^ (~c * (~a + ~b * ~x)) * sinh(~d / 2 + (~e * ~x) / 2) ^ (2 * ~n), ~x) <-- FreeQ([~F, ~a, ~b, ~c, ~d, ~e, ~f, ~g], ~x) && (EqQ(~f + ~g, 0) && ILtQ(~n, 0))
    @apply_utils Antiderivative((~F) ^ (~(c') * (~(a') + ~(b') * ~x)) * (~f + ~(g') * sinh(~(d') + ~(e') * ~x)) ^ ~(n'), ~x) => ((~f + ~g * sinh(~d + ~e * ~x)) ^ ~n / cosh((~d / 2 - (~f * pi) / (4 * ~g)) + (~e * ~x) / 2) ^ (2 * ~n)) * Antiderivative((~F) ^ (~c * (~a + ~b * ~x)) * cosh((~d / 2 - (~f * pi) / (4 * ~g)) + (~e * ~x) / 2) ^ (2 * ~n), ~x) <-- FreeQ([~F, ~a, ~b, ~c, ~d, ~e, ~f, ~g, ~n], ~x) && (EqQ((~f) ^ 2 + (~g) ^ 2, 0) && Not(IntegerQ(~n)))
    @apply_utils Antiderivative((~F) ^ (~(c') * (~(a') + ~(b') * ~x)) * (~f + ~(g') * cosh(~(d') + ~(e') * ~x)) ^ ~(n'), ~x) => ((~f + ~g * cosh(~d + ~e * ~x)) ^ ~n / cosh(~d / 2 + (~e * ~x) / 2) ^ (2 * ~n)) * Antiderivative((~F) ^ (~c * (~a + ~b * ~x)) * cosh(~d / 2 + (~e * ~x) / 2) ^ (2 * ~n), ~x) <-- FreeQ([~F, ~a, ~b, ~c, ~d, ~e, ~f, ~g, ~n], ~x) && (EqQ(~f - ~g, 0) && Not(IntegerQ(~n)))
    @apply_utils Antiderivative((~F) ^ (~(c') * (~(a') + ~(b') * ~x)) * (~f + ~(g') * cosh(~(d') + ~(e') * ~x)) ^ ~(n'), ~x) => ((~f + ~g * cosh(~d + ~e * ~x)) ^ ~n / sinh(~d / 2 + (~e * ~x) / 2) ^ (2 * ~n)) * Antiderivative((~F) ^ (~c * (~a + ~b * ~x)) * sinh(~d / 2 + (~e * ~x) / 2) ^ (2 * ~n), ~x) <-- FreeQ([~F, ~a, ~b, ~c, ~d, ~e, ~f, ~g, ~n], ~x) && (EqQ(~f + ~g, 0) && Not(IntegerQ(~n)))
    @apply_utils Antiderivative((~F) ^ (~(c') * (~(a') + ~(b') * ~x)) * cosh(~(d') + ~(e') * ~x) ^ ~(m') * (~f + ~(g') * sinh(~(d') + ~(e') * ~x)) ^ ~(n'), ~x) => (~g) ^ ~n * Antiderivative((~F) ^ (~c * (~a + ~b * ~x)) * tanh((~d / 2 + (~e * ~x) / 2) - (~f * pi) / (4 * ~g)) ^ ~m, ~x) <-- FreeQ([~F, ~a, ~b, ~c, ~d, ~e, ~f, ~g], ~x) && (EqQ((~f) ^ 2 + (~g) ^ 2, 0) && (IntegersQ(~m, ~n) && EqQ(~m + ~n, 0)))
    @apply_utils Antiderivative((~F) ^ (~(c') * (~(a') + ~(b') * ~x)) * sinh(~(d') + ~(e') * ~x) ^ ~(m') * (~f + ~(g') * cosh(~(d') + ~(e') * ~x)) ^ ~(n'), ~x) => (~g) ^ ~n * Antiderivative((~F) ^ (~c * (~a + ~b * ~x)) * tanh(~d / 2 + (~e * ~x) / 2) ^ ~m, ~x) <-- FreeQ([~F, ~a, ~b, ~c, ~d, ~e, ~f, ~g], ~x) && (EqQ(~f - ~g, 0) && (IntegersQ(~m, ~n) && EqQ(~m + ~n, 0)))
    @apply_utils Antiderivative((~F) ^ (~(c') * (~(a') + ~(b') * ~x)) * sinh(~(d') + ~(e') * ~x) ^ ~(m') * (~f + ~(g') * cosh(~(d') + ~(e') * ~x)) ^ ~(n'), ~x) => (~g) ^ ~n * Antiderivative((~F) ^ (~c * (~a + ~b * ~x)) * coth(~d / 2 + (~e * ~x) / 2) ^ ~m, ~x) <-- FreeQ([~F, ~a, ~b, ~c, ~d, ~e, ~f, ~g], ~x) && (EqQ(~f + ~g, 0) && (IntegersQ(~m, ~n) && EqQ(~m + ~n, 0)))
    @apply_utils Antiderivative(((~F) ^ (~(c') * (~(a') + ~(b') * ~x)) * (~h + ~(i') * cosh(~(d') + ~(e') * ~x))) / (~f + ~(g') * sinh(~(d') + ~(e') * ~x)), ~x) => 2 * ~i * Antiderivative((~F) ^ (~c * (~a + ~b * ~x)) * (cosh(~d + ~e * ~x) / (~f + ~g * sinh(~d + ~e * ~x))), ~x) + Antiderivative((~F) ^ (~c * (~a + ~b * ~x)) * ((~h - ~i * cosh(~d + ~e * ~x)) / (~f + ~g * sinh(~d + ~e * ~x))), ~x) <-- FreeQ([~F, ~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~i], ~x) && (EqQ((~f) ^ 2 + (~g) ^ 2, 0) && (EqQ((~h) ^ 2 - (~i) ^ 2, 0) && EqQ(~g * ~h - ~f * ~i, 0)))
    @apply_utils Antiderivative(((~F) ^ (~(c') * (~(a') + ~(b') * ~x)) * (~h + ~(i') * sinh(~(d') + ~(e') * ~x))) / (~f + ~(g') * cosh(~(d') + ~(e') * ~x)), ~x) => 2 * ~i * Antiderivative((~F) ^ (~c * (~a + ~b * ~x)) * (sinh(~d + ~e * ~x) / (~f + ~g * cosh(~d + ~e * ~x))), ~x) + Antiderivative((~F) ^ (~c * (~a + ~b * ~x)) * ((~h - ~i * sinh(~d + ~e * ~x)) / (~f + ~g * cosh(~d + ~e * ~x))), ~x) <-- FreeQ([~F, ~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~i], ~x) && (EqQ((~f) ^ 2 - (~g) ^ 2, 0) && (EqQ((~h) ^ 2 + (~i) ^ 2, 0) && EqQ(~g * ~h + ~f * ~i, 0)))
    @apply_utils Antiderivative((~F) ^ (~(c') * ~u) * (~G)(~v) ^ ~(n'), ~x) => Antiderivative((~F) ^ (~c * ExpandToSum(~u, ~x)) * (~G)(ExpandToSum(~v, ~x)) ^ ~n, ~x) <-- FreeQ([~F, ~c, ~n], ~x) && (HyperbolicQ(~G) && (LinearQ([~u, ~v], ~x) && Not(LinearMatchQ([~u, ~v], ~x))))
    @apply_utils Antiderivative((~(f') * ~x) ^ ~(m') * (~F) ^ (~(c') * (~(a') + ~(b') * ~x)) * sinh(~(d') + ~(e') * ~x) ^ ~(n'), ~x) => Module([u = IntHide((~F) ^ (~c * (~a + ~b * ~x)) * sinh(~d + ~e * ~x) ^ ~n, ~x)], Dist((~f * ~x) ^ ~m, u, ~x) - ~f * ~m * Antiderivative((~f * ~x) ^ (~m - 1) * u, ~x)) <-- FreeQ([~F, ~a, ~b, ~c, ~d, ~e, ~f], ~x) && (IGtQ(~n, 0) && GtQ(~m, 0))
    @apply_utils Antiderivative((~(f') * ~x) ^ ~(m') * (~F) ^ (~(c') * (~(a') + ~(b') * ~x)) * cosh(~(d') + ~(e') * ~x) ^ ~(n'), ~x) => Module([u = IntHide((~F) ^ (~c * (~a + ~b * ~x)) * cosh(~d + ~e * ~x) ^ ~n, ~x)], Dist((~f * ~x) ^ ~m, u, ~x) - ~f * ~m * Antiderivative((~f * ~x) ^ (~m - 1) * u, ~x)) <-- FreeQ([~F, ~a, ~b, ~c, ~d, ~e, ~f], ~x) && (IGtQ(~n, 0) && GtQ(~m, 0))
    @apply_utils Antiderivative((~(f') * ~x) ^ ~m * (~F) ^ (~(c') * (~(a') + ~(b') * ~x)) * sinh(~(d') + ~(e') * ~x), ~x) => (((~f * ~x) ^ (~m + 1) / (~f * (~m + 1))) * (~F) ^ (~c * (~a + ~b * ~x)) * sinh(~d + ~e * ~x) - (~e / (~f * (~m + 1))) * Antiderivative((~f * ~x) ^ (~m + 1) * (~F) ^ (~c * (~a + ~b * ~x)) * cosh(~d + ~e * ~x), ~x)) - ((~b * ~c * log(~F)) / (~f * (~m + 1))) * Antiderivative((~f * ~x) ^ (~m + 1) * (~F) ^ (~c * (~a + ~b * ~x)) * sinh(~d + ~e * ~x), ~x) <-- FreeQ([~F, ~a, ~b, ~c, ~d, ~e, ~f, ~m], ~x) && (LtQ(~m, -1) || SumSimplerQ(~m, 1))
    @apply_utils Antiderivative((~(f') * ~x) ^ ~m * (~F) ^ (~(c') * (~(a') + ~(b') * ~x)) * cosh(~(d') + ~(e') * ~x), ~x) => (((~f * ~x) ^ (~m + 1) / (~f * (~m + 1))) * (~F) ^ (~c * (~a + ~b * ~x)) * cosh(~d + ~e * ~x) - (~e / (~f * (~m + 1))) * Antiderivative((~f * ~x) ^ (~m + 1) * (~F) ^ (~c * (~a + ~b * ~x)) * sinh(~d + ~e * ~x), ~x)) - ((~b * ~c * log(~F)) / (~f * (~m + 1))) * Antiderivative((~f * ~x) ^ (~m + 1) * (~F) ^ (~c * (~a + ~b * ~x)) * cosh(~d + ~e * ~x), ~x) <-- FreeQ([~F, ~a, ~b, ~c, ~d, ~e, ~f, ~m], ~x) && (LtQ(~m, -1) || SumSimplerQ(~m, 1))
    #= Antiderivative(((~f')*(~x))^(~m')*(~F)^((~c')*((~a')+(~b')*(~x)))*sinh((~d')+(~e')*(~x))^(~n'),(~x)*Symbol) :=  (-1)^n/2^n*Antiderivative(ExpandIntegrand((f*x)^m*F^(c*(a+b*x)),(ℯ^(-(d+e*x))-ℯ^( d+e*x))^n,x),x) <-- FreeQ([F,a,b,c,d,e,f],x) && IGtQ(n,0) =#
    #= Antiderivative(((~f')*(~x))^(~m')*(~F)^((~c')*((~a')+(~b')*(~x)))*cosh((~d')+(~e')*(~x))^(~n'),(~x)*Symbol) :=  1/2^n*Antiderivative(ExpandIntegrand((f*x)^m*F^(c*(a+b*x)),(ℯ^(-(d+e*x))+ℯ^(d+e* x))^n,x),x) <-- FreeQ([F,a,b,c,d,e,f],x) && IGtQ(n,0) =#
    @apply_utils Antiderivative((~F) ^ (~(c') * (~(a') + ~(b') * ~x)) * sinh(~(d') + ~(e') * ~x) ^ ~(m') * cosh(~(f') + ~(g') * ~x) ^ ~(n'), ~x) => Antiderivative(ExpandTrigReduce((~F) ^ (~c * (~a + ~b * ~x)), sinh(~d + ~e * ~x) ^ ~m * cosh(~f + ~g * ~x) ^ ~n, ~x), ~x) <-- FreeQ([~F, ~a, ~b, ~c, ~d, ~e, ~f, ~g], ~x) && (IGtQ(~m, 0) && IGtQ(~n, 0))
    @apply_utils Antiderivative((~x) ^ ~(p') * (~F) ^ (~(c') * (~(a') + ~(b') * ~x)) * sinh(~(d') + ~(e') * ~x) ^ ~(m') * cosh(~(f') + ~(g') * ~x) ^ ~(n'), ~x) => Antiderivative(ExpandTrigReduce((~x) ^ ~p * (~F) ^ (~c * (~a + ~b * ~x)), sinh(~d + ~e * ~x) ^ ~m * cosh(~f + ~g * ~x) ^ ~n, ~x), ~x) <-- FreeQ([~F, ~a, ~b, ~c, ~d, ~e, ~f, ~g], ~x) && (IGtQ(~m, 0) && (IGtQ(~n, 0) && IGtQ(~p, 0)))
    @apply_utils Antiderivative((~F) ^ (~(c') * (~(a') + ~(b') * ~x)) * (~G)(~(d') + ~(e') * ~x) ^ ~(m') * (~H)(~(d') + ~(e') * ~x) ^ ~(n'), ~x) => Antiderivative(ExpandTrigToExp((~F) ^ (~c * (~a + ~b * ~x)), (~G)(~d + ~e * ~x) ^ ~m * (~H)(~d + ~e * ~x) ^ ~n, ~x), ~x) <-- FreeQ([~F, ~a, ~b, ~c, ~d, ~e], ~x) && (IGtQ(~m, 0) && (IGtQ(~n, 0) && (HyperbolicQ(~G) && HyperbolicQ(~H))))
    @apply_utils Antiderivative((~F) ^ ~u * sinh(~v) ^ ~(n'), ~x) => Antiderivative(ExpandTrigToExp((~F) ^ ~u, sinh(~v) ^ ~n, ~x), ~x) <-- FreeQ(~F, ~x) && ((LinearQ(~u, ~x) || PolyQ(~u, ~x, 2)) && ((LinearQ(~v, ~x) || PolyQ(~v, ~x, 2)) && IGtQ(~n, 0)))
    @apply_utils Antiderivative((~F) ^ ~u * cosh(~v) ^ ~(n'), ~x) => Antiderivative(ExpandTrigToExp((~F) ^ ~u, cosh(~v) ^ ~n, ~x), ~x) <-- FreeQ(~F, ~x) && ((LinearQ(~u, ~x) || PolyQ(~u, ~x, 2)) && ((LinearQ(~v, ~x) || PolyQ(~v, ~x, 2)) && IGtQ(~n, 0)))
    @apply_utils Antiderivative((~F) ^ ~u * sinh(~v) ^ ~(m') * cosh(~v) ^ ~(n'), ~x) => Antiderivative(ExpandTrigToExp((~F) ^ ~u, sinh(~v) ^ ~m * cosh(~v) ^ ~n, ~x), ~x) <-- FreeQ(~F, ~x) && ((LinearQ(~u, ~x) || PolyQ(~u, ~x, 2)) && ((LinearQ(~v, ~x) || PolyQ(~v, ~x, 2)) && (IGtQ(~m, 0) && IGtQ(~n, 0))))
end

max_stars_repo_path: src/rules/6 Hyperbolic functions/6.7 Miscellaneous/6.7.7 F^(c (a+b x)) hyper(d+e x)^n.jl


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule LAPACK_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("LAPACK")
JLLWrappers.@generate_main_file("LAPACK", UUID("51474c39-65e3-53ba-86ba-03b1b862ec14"))
end  # module LAPACK_jll


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule SoapyLMS7_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("SoapyLMS7")
JLLWrappers.@generate_main_file("SoapyLMS7", UUID("e7ed14a0-13eb-5dc7-93e2-6133b2eb8bed"))
end  # module SoapyLMS7_jll


text:
using PolyChaos

mutable struct PolynomialChaosSurrogate{X,Y,L,U,C,O,N} <: AbstractSurrogate
    x::X
    y::Y
    lb::L
    ub::U
    coeff::C
    ortopolys::O
    num_of_multi_indexes::N
 end

function _calculatepce_coeff(x,y,num_of_multi_indexes,op::AbstractCanonicalOrthoPoly)
    n = length(x)
    A = zeros(eltype(x),n,num_of_multi_indexes)
    for i = 1:n
        A[i,:] = PolyChaos.evaluate(x[i],op)
    end
    return (A'*A)\(A'*y)
end

function PolynomialChaosSurrogate(x,y,lb::Number,ub::Number; op::AbstractCanonicalOrthoPoly = GaussOrthoPoly(2))
    n = length(x)
    poly_degree = op.deg
    num_of_multi_indexes = 1+poly_degree
    if n < 2+3*num_of_multi_indexes
        throw("To avoid numerical problems, it's strongly suggested to have at least $(2+3*num_of_multi_indexes) samples")
    end
    coeff = _calculatepce_coeff(x,y,num_of_multi_indexes,op)
    return PolynomialChaosSurrogate(x,y,lb,ub,coeff,op,num_of_multi_indexes)

end

function (pc::PolynomialChaosSurrogate)(val::Number)
    return sum([pc.coeff[i]*PolyChaos.evaluate(val,pc.ortopolys)[i] for i = 1:pc.num_of_multi_indexes])
end

function _calculatepce_coeff(x,y,num_of_multi_indexes,op::MultiOrthoPoly)
    n = length(x)
    d = length(x[1])
    A = zeros(eltype(x[1]),n,num_of_multi_indexes)
    for i = 1:n
        xi = zeros(eltype(x[1]),d)
        for j = 1:d
            xi[j] = x[i][j]
        end
        A[i,:] = PolyChaos.evaluate(xi,op)
    end
    return (A'*A)\(A'*y)
end

function PolynomialChaosSurrogate(x,y,lb,ub; op::MultiOrthoPoly = MultiOrthoPoly([GaussOrthoPoly(2) for j = 1:length(lb)],2))
    n = length(x)
    d = length(lb)
    poly_degree = op.deg
    num_of_multi_indexes = binomial(d+poly_degree,poly_degree)
    if n < 2+3*num_of_multi_indexes
        throw("To avoid numerical problems, it's strongly suggested to have at least $(2+3*num_of_multi_indexes) samples")
    end
    coeff = _calculatepce_coeff(x,y,num_of_multi_indexes,op)
    return PolynomialChaosSurrogate(x,y,lb,ub,coeff,op,num_of_multi_indexes)

end

function (pcND::PolynomialChaosSurrogate)(val)
    sum = zero(eltype(val[1]))
    for i = 1:pcND.num_of_multi_indexes
        sum = sum + pcND.coeff[i]*first(PolyChaos.evaluate(pcND.ortopolys.ind[i,:],collect(val),pcND.ortopolys))
    end
    return sum
end

function add_point!(polych::PolynomialChaosSurrogate,x_new,y_new)
    if length(polych.lb) == 1
        #1D
        polych.x = vcat(polych.x,x_new)
        polych.y = vcat(polych.y,y_new)
        polych.coeff = _calculatepce_coeff(polych.x,polych.y,polych.num_of_multi_indexes,polych.ortopolys)
    else
        polych.x = vcat(polych.x,x_new)
        polych.y = vcat(polych.y,y_new)
        polych.coeff = _calculatepce_coeff(polych.x,polych.y,polych.num_of_multi_indexes,polych.ortopolys)
    end
    nothing
end


text:
@testset "base" begin
    @testset "Trig" begin
        @testset "Basics" for x = (Float64(π)-0.01, Complex(π, π/2))
            test_scalar(sec, x)
            test_scalar(csc, x)
            test_scalar(cot, x)
            test_scalar(sinpi, x)
            test_scalar(cospi, x)
        end
        @testset "Hyperbolic" for x = (Float64(π)-0.01, Complex(π-0.01, π/2))
            test_scalar(sech, x)
            test_scalar(csch, x)
            test_scalar(coth, x)
        end
        @testset "Degrees" begin
            x = 45.0
            test_scalar(sind, x)
            test_scalar(cosd, x)
            test_scalar(tand, x)
            test_scalar(secd, x)
            test_scalar(cscd, x)
            test_scalar(cotd, x)
        end
        @testset "Inverses" for x = (0.5, Complex(0.5, 0.25))
            test_scalar(asec, 1/x)
            test_scalar(acsc, 1/x)
            test_scalar(acot, 1/x)
        end
        @testset "Inverse hyperbolic" for x = (0.5, Complex(0.5, 0.25),  Complex(-2.1 -3.1im))
            test_scalar(asinh, x)
            test_scalar(acosh, x + 1)  # +1 accounts for domain for real
            test_scalar(atanh, x)
            test_scalar(asech, x)
            test_scalar(acsch, x)
            test_scalar(acoth, x + 1)
        end

        @testset "Inverse degrees" for x = (0.5, Complex(0.5, 0.25))
            test_scalar(asind, x)
            test_scalar(acosd, x)
            test_scalar(atand, x)
            test_scalar(asecd, 1/x)
            test_scalar(acscd, 1/x)
            test_scalar(acotd, 1/x)
        end
    end  # Trig

    @testset "Angles" begin
        for x in (-0.1, 6.4)
            test_scalar(deg2rad, x)
            test_scalar(rad2deg, x)
        end
    end

    @testset "Unary complex functions" begin
        for x in (-4.1, 6.4)
            test_scalar(real, x)
            test_scalar(imag, x)
            test_scalar(hypot, x)
            test_scalar(adjoint, x)
        end
    end

    @testset "*(x, y) (scalar)" begin
        # This is pretty important so testing it fairly heavily
        test_points = (0.0, -2.1, 3.2, 3.7+2.12im, 14.2-7.1im)
        @testset "$x * $y; (perturbed by: $perturb)" for
            x in test_points, y in test_points, perturb in test_points

            # give small off-set so as can't slip in symmetry
            x̄ = ẋ = 0.5 + perturb
            ȳ = ẏ = 0.6 + perturb
            Δz = perturb

            frule_test(*, (x, ẋ), (y, ẏ))
            rrule_test(*, Δz, (x, x̄), (y, ȳ))
        end
    end

    @testset "matmul *(x, y)" begin
        x, y = rand(3, 2), rand(2, 5)
        z, pullback = rrule(*, x, y)

        @test z == x * y

        z̄ = rand(3, 5)
        (ds, dx, dy) = pullback(z̄)

        @test ds === NO_FIELDS

        @test extern(dx) == extern(zeros(3, 2) .+ dx)
        @test extern(dy) == extern(zeros(2, 5) .+ dy)
    end

    @testset "binary function ($f)" for f in (mod, \)
        x, Δx, x̄ = 10rand(3)
        y, Δy, ȳ = rand(3)
        Δz = rand()

        frule_test(f, (x, Δx), (y, Δy))
        rrule_test(f, Δz, (x, x̄), (y, ȳ))
    end

    @testset "x^n for x<0" begin
        x = -15*rand()
        Δx, x̄ = 10rand(2)
        y, Δy, ȳ = rand(3)
        Δz = rand()

        frule_test(^, (-x, Δx), (y, Δy))
        rrule_test(^, Δz, (-x, x̄), (y, ȳ))
    end

    @testset "identity" begin
        rrule_test(identity, randn(), (randn(), randn()))
        rrule_test(identity, randn(4), (randn(4), randn(4)))

        rrule_test(
            identity, Tuple(randn(3)),
            (Composite{Tuple}(randn(3)...), Composite{Tuple}(randn(3)...))
        )
    end

    @testset "Constants" for x in (-0.1, 6.4, 1.0+0.5im, -10.0+0im, 0+200im)
        test_scalar(one, x)
        test_scalar(zero, x)
    end

    @testset "trinary ($f)" for f in (muladd, fma)
        x, Δx, x̄ = 10randn(3)
        y, Δy, ȳ = randn(3)
        z, Δz, z̄ = randn(3)
        Δk = randn()

        frule_test(f, (x, Δx), (y, Δy), (z, Δz))
        rrule_test(f, Δk, (x, x̄), (y, ȳ), (z, z̄))
    end
end


text:
#!/usr/bin/env julia

println("Loading dependencies...")
using Plots
using LibSerialPort
using JSON

# Plot backend
gr(size=(900, 500))

function readmcu(mcu::SerialPort)
    datadict = Dict()
    msg = readstring(mcu)
    if startswith(msg, "{") && endswith(msg, "}")
        try
            datadict = JSON.parse(msg)
        catch
            println("skipping $msg")
        end
    elseif length(msg) > 0
        println(msg) # Pass through non-data messages
    end
    return datadict
end

# Manage arrays as FIFO queues
function bump_queues(datadict, times, setpoints, positions)

    for key in ["time", "setpoint", "x"]
        haskey(datadict, key) || return false
    end

    push!(times, datadict["time"]/1000) # ms to seconds
    push!(setpoints, datadict["setpoint"])
    push!(positions, datadict["x"])
    for v in (times, setpoints, positions)
        shift!(v)
    end
    return true
end

function rtplot(times, setpoints, positions)
    plot(times, [setpoints, positions],
         xaxis="time [s]",
         yaxis=("position [% max]", (-110, 110)), # title, limits
         labels=(["setpoint","position"]), # legend
         line=([1 2], [:steppre :path], ["gray" "crimson"], [1.0 0.5]), # lineweights, linetypes, colors, alphas
         fill=(0, [0.1 0], "dodgerblue")) # y-origin, alpha, colors

    # Update the plot window
    gui()
end

let
    counter = 0
    global process_mcu_stream

    function process_mcu_stream(mcu, times, setpoints, positions)
        datadict = readmcu(mcu)
        counter += 1
        if counter > 100
            bump_queues(datadict, times, setpoints, positions) || return
            rtplot(times, setpoints, positions)
        end
    end
end

function main()

    # Number of points to display
    const N = 500
    const address = length(ARGS) == 1 ? ARGS[1] : "/dev/cu.usbmodem2370891"
    times, setpoints, positions = collect(linspace(0,10,N)), zeros(N), zeros(N)

    println("Opening plot window...")
    rtplot(times, setpoints, positions)
    sleep(3)

    println("Connecting to MCU...")
    mcu = open(address, 115200)

    println("Starting I/O loop.")
    while true
        @async begin
            a = readline(STDIN)
            a == "q\n" && quit()
            write(mcu, "$a")
        end
        @async begin
            process_mcu_stream(mcu, times, setpoints, positions)
        end
        # Give the queued tasks a chance to run
        sleep(0.0001)
    end
end

main()

max_stars_repo_name: andrewadare/inverted-pendulum

max_issues_repo_name: andrewadare/inverted-pendulum

max_forks_repo_name: andrewadare/inverted-pendulum


text:
using Pkg

Pkg.instantiate()

using Test, UUIDs, HTTP

const REGISTRY_PATH = joinpath(@__DIR__, "..", "Registry.toml")
const VALID_METHODS = ("hosted", "git-repo", "vendored")

isurlish(url) = occursin(r"(http|https|git)\://.*", url)

function testkeyvalidity(registry, key)
    uuid = UUID(key)

    entry = registry[key]
    name = entry["name"]

    method = ""
    if haskey(entry, "method")
        method = entry["method"]
        if !(method in VALID_METHODS)
            @error("Invalid method `$(method)` specified for `$(name)` ($(uuid)).")
            return false
        end
    else
        @error("No `method` specified for `$(name)` ($(uuid)).")
        return false
    end

    if haskey(entry, "location")
        uri = entry["location"]
        if method in ("git-repo", "hosted")
            if isurlish(uri)
                try
                    req = HTTP.request("GET", uri)
                    if req.status > 400
                        @error("""
                              `$(uri)` returned `$(req.status)` for `$(name)` ($(uuid)).
                              Please double check the URL.
                              """)
                        return false
                    end
                catch err
                    @error("""
                          `$(uri)` request failed for `$(name)` ($(uuid)).
                          Please double check the URL.
                          """)
                    return false
                end
            else
                @error("Invalid URL `$(uri)` for `$(name)` ($(uuid)).")
                return false
            end
        elseif method == "vendored"
            if isurlish(uri)
                @error("`location` should not be a URL for vendored docs in `$(name)` ($(uuid)).")
                return false
            end
        end
    else
        @error("No `location` specified for `$(name)` ($(uuid)).")
        return false
    end

    return true
end

@testset "Registry Validity" begin
    @test isfile(REGISTRY_PATH)
    toml = Pkg.TOML.parsefile(REGISTRY_PATH)
    @test collect(unique(keys(toml))) == collect(keys(toml))
    @testset "Package: $(toml[key]["name"]) ($(key))" for key in keys(toml)
        isvalid = testkeyvalidity(toml, key)
        @test isvalid
    end
end

max_forks_repo_name: pfitzseb/DocumentationGeneratorRegistry


text:
export RMSE, MAE

"""
    RMSE

Root Mean Squared Error.

    measure(
        metric::RMSE,
        truth::AbstractVector,
        pred::AbstractVector
    )
"""
struct RMSE <: AccuracyMetric end
function measure(metric::RMSE, truth::AbstractVector, pred::AbstractVector)
    sqrt(sum((truth - pred).^2) / length(truth))
end

"""
    MAE

Mean Absolute Error.

    measure(
        metric::MAE,
        truth::AbstractVector,
        pred::AbstractVector
    )
"""
struct MAE <: AccuracyMetric end
function measure(metric::MAE, truth::AbstractVector, pred::AbstractVector)
    sum(abs.(truth - pred)) / length(truth)
end


text:
# wrap_cublas.jl
#
# Script to generate CUBLAS wrapper code. Generated code in libcublas_types.jl
# file needs to be modified before including in the CUBLAS module. Code is
# generated here and manually moved to the ../src directory.
#
# This script works with Clang.jl v0.0.1.  To install and use on Fedora 20:
#   julia> ENV["LLVM_CONFIG"]="/usr/bin/llvm-config"
#   julia> Pkg.add("Clang")
#   julia> Pkg.pin("Clang",v"0.0.1")
#   julia> Pkg.build("Clang")
#
# To run the script execute the following from a shell in this directory:
#   $ julia wrap_cublas.jl
#
# Generated files:
#   libcublas_types.jl: contains julia representations of enums and types in
#                       cublas_v2.h
#   libcublas.jl: contains julia wrappers to functions in libcublas.so
#
# Note: this script needs to be changed for Clang.jl v0.0.2, which is in
# METADATA.jl at the time of this writing.
#
# Author: Nick Henderson <nwh@stanford.edu>
# Created: 2014-08-26
# License: MIT
#

using Clang, Clang.cindex
import Clang.wrap_c: repr_jl, rep_type, rep_args, name_safe

# main cublas header
headers = ["/usr/local/cuda/include/cublas_v2.h"]

# may need to be modified to host system
includes = ["/usr/include",
            "/usr/lib/gcc/x86_64-redhat-linux/4.8.3"]

# skip cursors with empty names, this avoids anonymous enums
function check_anon(cursorname, cursor)
    if isempty(cursorname)
        return false
    end
    return true
end

# overload wrap_c.wrap to ignore CUBLAS macro definitions
function wrap_c.wrap(strm::IO, md::cindex.MacroDefinition)
    warn("macro definition")
end

# Customize the wrap function for functions. This was copied
# from Clang/src/wrap_c.jl, with the following customizations:
#   - error-check functions that return a cublasStatus_t
#   - omit types from function prototypes
# This function was copied from CUDArt.jl/gen/wrap_cuda.jl
skipcheck = []
function wrap_c.wrap(buf::IO, funcdecl::FunctionDecl, libname::ASCIIString)
    function print_args(buf::IO, cursors, types)
        i = 1
        for c in cursors
            print(buf, name_safe(c))
            (i < length(cursors)) && print(buf, ", ")
            i += 1
        end
    end

    cu_spelling = spelling(funcdecl)

    funcname = spelling(funcdecl)
    arg_types = cindex.function_args(funcdecl)
    args = [x for x in search(funcdecl, ParmDecl)]
    arg_list = tuple( [repr_jl(x) for x in arg_types]... )
    ret_type = repr_jl(return_type(funcdecl))

    print(buf, "function ")
    fname = spelling(funcdecl)
    print(buf, fname)
    print(buf, "(")
    print_args(buf, args, [myrepr_jl(x) for x in arg_types])
    println(buf, ")")
    print(buf, "  ")
    checkerr = ret_type == "cublasStatus_t" && !in(fname,skipcheck)
    checkerr && print(buf, "statuscheck(")
    print(buf, "ccall( (:", funcname, ", ", libname, "), ")
    print(buf, rep_type(ret_type))
    print(buf, ", ")
    print(buf, rep_args(arg_list), ", ")
    for (i,arg) in enumerate(args)
        print(buf, name_safe(arg))
        (i < length(args)) && print(buf, ", ")
    end
    checkerr && print(buf, ")")
    println(buf, ")")
    println(buf, "end")
end

function myrepr_jl(x)
    str = repr_jl(x)
    return str
    #return (str == "Ptr{Cint}") ? "Array{Cint}" : str
end

# initialize wrap context
context = wrap_c.init(output_file="libcublas.jl",
                      common_file="libcublas_types.jl",
                      header_library=x->"libcublas",
                      clang_includes=includes,
                      cursor_wrapped=check_anon,
                      header_wrapped=(x,y)->(contains(y,"cublas")))

# generate the wrapper
wrap_c.wrap_c_headers(context,headers)


text:
################################################################################
#
#  NfFactorBase : Factor bases for number fields
#  A factor basis is mostly a collection of prime ideals, designed,
#  if possible, to allow for rapid testing if elements are smooth.
#
################################################################################

function NfFactorBase(O::NfOrd, B::Int, F::Function, complete::Bool = false, degree_limit::Int = 0)
  @vprint :ClassGroup 2 "Splitting the prime ideals ...\n"
  lp = prime_ideals_up_to(O, B, F, complete = complete, degree_limit = degree_limit)
  @vprint :ClassGroup 2 " done \n"
  return NfFactorBase(O, lp)
end

function NfFactorBase(O::NfOrd, lp::AbstractVector{Int}, degree_limit::Int = 0)
  @vprint :ClassGroup 2 "Splitting the prime ideals ...\n"
  lP = prime_ideals_over(O, lp, degree_limit = degree_limit)
  @vprint :ClassGroup 2 " done \n"
  return NfFactorBase(O, lP)
end

function NfFactorBase(O::NfOrd, B::Int;
                        complete::Bool = true, degree_limit::Int = 5)
  @vprint :ClassGroup 2 "Splitting the prime ideals ...\n"
  lp = prime_ideals_up_to(O, B, complete = complete, degree_limit = degree_limit)
  @vprint :ClassGroup 2 " done \n"
  return NfFactorBase(O, lp)
end

function NfFactorBase(O::NfOrd, lp::Vector{NfOrdIdl})
  lp = sort(lp, lt = function(a,b) return norm(a) > norm(b); end)
  FB = NfFactorBase()
  FB.size = length(lp)
  FB.ideals = lp

  # Magic constant 20?
  FB.rw = Array{Int}(undef, 20)
  FB.mx = 20

  fb = Dict{fmpz, Vector{Tuple{Int, NfOrdIdl}}}()

  for i = 1:length(lp)
    if !haskey(fb, minimum(lp[i]))
      fb[minimum(lp[i])] = [(i, lp[i])]
    else
      push!(fb[minimum(lp[i])], (i, lp[i]))
    end
  end

  FB.fb = Dict{fmpz, FactorBaseSingleP}()
  for (p, v) in fb
    if fits(Int, p)
      FB.fb[p] = FactorBaseSingleP(Int(p), v)
    else
      FB.fb[p] = FactorBaseSingleP(p, v)
    end
  end

  FB.fb_int = FactorBase(Set(keys(FB.fb)); check = false)

  return FB
end

function ring(F::NfFactorBase)
  return order(F.ideals[1])
end

order(F::NfFactorBase) = ring(F)

################################################################################
#
#  Factor number field element over factor base. Put valuations into row i of
#  the relation matrix M. The matrix M needs to have at least as many columns
#  as the FB has ideals.
#
################################################################################

function factor!(M::SMat{T}, i::Int, FB::NfFactorBase, a::nf_elem;
                 error = true, n = abs(norm(a))) where T
  fl, res = _factor(FB, a, error=error, n=n)
  if fl
    push!(M, res)
  end
  return fl
end

function _factor!(FB::NfFactorBase, a::nf_elem,
                    error::Bool = true, n::fmpq = abs(norm(a)), integral::Bool = true)
  T = fmpz
  O = order(FB.ideals[1])
  n = deepcopy(n)

  if integral
    df =numerator(n)
  else
    df = numerator(n)*denominator(a, O)
  end
  if isone(df)
    return true, sparse_row(FlintZZ)
  end

  d = factor(FB.fb_int, df, error)  #careful: if df is non-int-smooth, then error is ignored

  rw = FB.rw
  r = Vector{Tuple{Int, Int}}()
  ret = true
  for p in keys(d)
    vp = valuation!(n, p)
#    s::Vector{Tuple{Int, Int}}, vp::Int = FB.fb[p].doit(a, vp)
    s::Vector{Tuple{Int, Int}}, vp::Int = fb_doit(a, vp, FB.fb[p], fmpq(p)^vp)
    if !iszero(vp)
      ret = false
      if error
        @assert iszero(vp)
      end
    end
    r = vcat(r, s)
  end
  lg::Int = length(r)
  if lg > 0
    if length(rw) > FB.mx
      FB.mx = length(rw)
    end
    sort!(r, lt=function(a,b) return a[1] < b[1]; end)
    if ret
      @hassert :ClassGroup 9000  ideal(O, a) == prod([FB.ideals[i]^j for (i, j) in r])
    end
    @hassert :ClassGroup 1 length(r) > 0
    return ret, sparse_row(FlintZZ, r)
  else
    # factor failed or I have a unit.
    # sparse rel mat must not have zero-rows.
    return false, sparse_row(FlintZZ)
  end
end

function factor(FB::NfFactorBase, a::nf_elem)
  return _factor!(FB, a, true, abs(norm(a)), false)[2]
end

function _factor!(FB::Hecke.NfFactorBase, A::Hecke.NfOrdIdl,
                    error::Bool = true)
  T = fmpz
  O = order(A)

  n = norm(A)

  # If the ideal is the trivial ideal, return true, and the zero row
  # Otherwise factor will choke
  if isone(n)
    return true, sparse_row(FlintZZ)
  end

  d = factor(FB.fb_int, n) # as above: fails - even if error is false -
  # if the norm is not smooth

  rw = FB.rw
  r = Vector{Tuple{Int, Int}}()
  for p in keys(d)
    vp = valuation(n, p)
    s = Vector{Tuple{Int, Int}}()
    for P=FB.fb[p].lp
      v = valuation(A, P[2])
      if v != 0
        push!(s, (P[1], v))
        vp -= v*P[2].splitting_type[2]
        if iszero(vp)
          break
        end
      end
    end
    if vp != 0
      if error
        @assert vp == 0
      end
      return false, sparse_row(FlintZZ)
    end
    r = vcat(r, s)
  end
  lg::Int = length(r)
  if lg > 0
    if length(rw) > FB.mx
      FB.mx = length(rw)
    end
    sort!(r, lt = (a,b) -> a[1] < b[1])
    res = sparse_row(FlintZZ, r)
    @hassert :ClassGroup 9000 A == prod([FB.ideals[i]^j for (i, j) in r])
    @hassert :ClassGroup 1 length(r) > 0
    return true, res
  else
    # factor failed or I have a unit.
    # sparse rel mat must not have zero-rows.
    return false, sparse_row(FlintZZ)
  end
end


text:
# implementations of Graphs.jl interface methods

using Graphs
using Iterators

Graphs.edge_index(pm::PlayerMove) = error("not implemented")
Graphs.source(pm::PlayerMove) = pm.from
Graphs.target(pm::PlayerMove) = pm.to

# basic
Graphs.is_directed(game::Game) = true
Graphs.vertex_type(game::Game) = Loc
Graphs.edge_type(game::Game) = PlayerMove

# vertex map
Graphs.vertex_index(loc::Loc, game::Game) = 9*(loc.x-1) + loc.y
Graphs.implements_vertex_map(game::Game) = true

# vertex list
Graphs.num_vertices(game::Game) = 81
Graphs.vertices(game::Game) = map(Loc, product(1:9, 1:9))
Graphs.implements_vertex_list(game::Game) = true

# # edge map
Graphs.edge_index(edge::PlayerMove, game::Game) = error("not implemented")
Graphs.implements_edge_map(game::Game) = true

# incidence_list
Graphs.out_degree(loc::Loc, game::Game) = length(out_edges(loc, game))
function Graphs.out_edges(loc::Loc, game::Game)
    ret = PlayerMove[]
    for dir in DIRS  # this comes from loc.jl just fyi
        adj = dir(loc)
        if neighbors(game, loc, adj)
            if !occupied(game, adj) && inbounds(adj)
                push!(ret, PlayerMove(loc, adj))
            else
                # consider jumps
                #
                # first we need to check if it's possible to jump
                # directly over the piece occupying that space
                #
                # one thing to note is that we don't need to consider
                # the possibility of a space we're jumping to being
                # occupied (since we already know where both players
                # are)
                straight_jump = dir(adj)
                adj, straight_jump
                if (neighbors(game, adj, straight_jump) &&
                    inbounds(straight_jump))
                    push!(ret, PlayerMove(loc, straight_jump))
                else
                    # consider each of the two possible perpendicular jumps
                    for perp_dir in PERPENDICULAR[dir]
                        perp_jump = perp_dir(adj)
                        dir, perp_dir, adj, perp_jump
                        if (neighbors(game, adj, perp_jump) &&
                            inbounds(perp_jump))
                            push!(ret, PlayerMove(loc, perp_jump))
                        end
                    end
                end
            end
        end
    end
    return ret
end
Graphs.source(pm::PlayerMove, game::Game) = pm.from
Graphs.target(pm::PlayerMove, game::Game) = pm.to
Graphs.implements_incidence_list(game::Game) = true


text:
using  ApproxFunRational, AbstractIterativeSolvers
### Vector case
tol = 1.e-9
Ds1 = 0.5
Ds2 = 1.
Δs = Ds2-Ds1
Dn1 = 1.
Dn2 = 0.5
Δn = Dn2 - Dn1
rn = 1.
s0n = 1/sqrt(2.)
rs = 1.
s0s = 1/sqrt(2.)
pnn = 1
pmm = .2
L = 1.
D = 1. # overall delay

fPss = z -> 2*rs*s0s^2.0/(rs^2+z.^2)
#fPss = z -> .1*exp(-z^2/2)
fPnn = z -> pnn
fPmm = z -> pmm
Pss = Fun(zai(fPss),OscLaurent(0.0,L))
One = pad(Fun(1.,OscLaurent(0.0,L)),3)
Pnn = pnn*One
Pmm = pmm*One

H = zeros(Fun,2,2)
a = Pss + Pnn + Pmm - One
b = Pss*Fun(1.,OscLaurent(-Δs,L)) + Pnn*Fun(1.,OscLaurent(-Δn,L))
bt = Pss*Fun(1.,OscLaurent(Δs,L)) + Pnn*Fun(1.,OscLaurent(Δn,L))
H[1,1] = copy(a)
H[2,2] = copy(a)
H[1,2] = copy(b)
H[2,1] = copy(bt)
G = map(SumFun,copy(H))

# Set up right-hand side
b1 = Pss*Fun(1.,OscLaurent(- D + Ds1,L))
b2 = Pss*Fun(1.,OscLaurent(- D + Ds2,L))
h = map(SumFun,[b1,b2])

𝓒 = Cauchy(-1)
𝓢 = Cauchy(1)
simp(f) = combine!(chop!(f,tol))
inner(a,b) = ⋅(simp(a),simp(b))

function op(x::Array{T,1}) where T<:SumFun
    println("Apply Cauchy:")
    @time y = 𝓒*x
    println("Apply G:")
    @time y = G*y
    println("Subtract:")
    @time y = (x - y)
    return y
end

sop = x -> ( x - G*(𝓒*x))
out = GMRES(op,h,⋅,10*tol,30,simp)

u = sum([out[2][i]*out[1][i] for i=1:length(out[2])])

si_op(out[1][end])

𝓕 = FourierTransform(-1.0)
U = 𝓕*u

x = 0:.01:10
y1 = real(map(U[1],x))
y2 = real(map(U[2],x))

for i = 1:length(y1)-1
    if isinf(y1[i])
        y1[i] = y1[i+1]
    end
    if isinf(y2[i])
        y2[i] = y2[i+1]
    end
end

plot(x,y1)
plot!(x,y2)


text:
using VectorAutoregressions
using CSV, DataFrames
using Test

cd(@__DIR__)

@testset "VAR" begin
    Y = DataFrame(CSV.File("testdata1.csv"))
    v = VAR(Y, 2)

    @test coef(v) == v.B
    @test stderror(v) == v.seB
    @test residuals(v) == v.U
    @test fitted(v) == v.Yhat

    Bhat = [-0.017  0.016  0.013
            -0.320  0.044 -0.002
             0.146 -0.153  0.225
             0.961  0.289 -0.264
            -0.161  0.050  0.034
             0.115  0.019  0.355
             0.934 -0.010 -0.022]
    @test isapprox(coef(v), Bhat, atol=2e-3)

    sigu = [0.00212963  7.16167e-5   0.00012324
            7.16167e-5  0.000137338  6.14587e-5
            0.00012324  6.14587e-5   8.92035e-5]
    @test isapprox(v.ΣU, sigu, atol=1e-8)
end


text:
function ct_kalman_simple(T, Z, Q, E, mean_0, var_0, data_y, dt)
    # Compute Kalman filter of a continuous time process sampled at discrete time
    # ported from SeHyoun Ahn
    #
    #   & dx_t = Tx_t dt + R dW_t \\
    #   & y_{t_k} = Zx_{t_k} + v_{t_k} \\
    #   & x_0 \sim N(mean_0, var_0) \\
    #   & k\in\{1, 2, \ldots, n\} \\
    #   & W_t = \text{Brownian Motion} \\
    #   & v_k \sim N(0, E) \\
    #   & dt_i = t_{i} - t_{i-1}
    #
    # Args:
    #   T (matrix): see equation above
    #   Z (matrix): see equation above
    #   Q (matrix): :math:`Q = RR'` i.e., variance
    #   E (matrix): see equation above
    #   mean_0 (vector): see equation above
    #   var_0 (matrix): see equation above
    #   data_y (vector): :math:`y_k` above
    #   dt (vector): time period between data measurement (see above)
    #
    # Note:
    #   - The data does not have to come with uniform step size.
    #   - There is a `kalman <https://www.mathworks.com/help/control/ref/kalman.html>`_ function in MATLAB as well, but it is under control system toolbox.
    #
    # Returns:
    #   : [mean_out, var_out, prob_out]
    #
    #   - **mean_out** (vector): :math:`E(x_{t_n}| y_{t_1},\ldots, y_{t_n})`
    #   - **var_out** (matrix): :math:`var(x_{t_n}| y_{t_1},\ldots, y_{t_n})`
    #   - **prob_out** (double): log-likelihood of the data
    #
    # Reference:
    #
    #   - `Wikipedia article on Kalman Filtering <https://en.wikipedia.org/wiki/Kalman_filter#Hybrid_Kalman_filter>`_
    #   - :cite:`grewal2011kalman`
    #
    # .. Syntax:
    #
    # ..   [mean_out, var_out, prob_out] = kalman_ct(T, Z, Q, E, mean_0, var_0, data_y, dt)
    #

    # Get data size
    n_data = size(data_y, 1)
    if (length(dt) < n_data)
        dt = repeat(dt, n_data)
    end
    n_dim = length(mean_0)
    prob_out = 0
    mean_out = mean_0
    var_out = var_0

    # Step through data points

    for iter_data = 1:n_data

        #       tspan = (0.0,dt[iter_data][1])
        # s0 = mean_0 #initial values
        # P0 = var_out
        #        f(u,p,t) = T*u
        #     g(u,p,t) = T*u + u*T' + Q
        #           s_prob = ODEProblem(f, mean_out, tspan)
        #     P_prob = ODEProblem(g, P0, tspan)
        # sol = solve(s_prob,Tsit5(),reltol=1e-4,abstol=1e-4)

        ## forecasting state (this is different for ct)
        mean_out = exp(dt[iter_data][1]*T)*mean_out
        AB = exp(dt[iter_data][1]*[hcat(T, Q); hcat(zeros(n_dim,n_dim), -T')])*[var_out; Matrix{Float64}(I,n_dim,n_dim)]
        var_inter = AB[1:n_dim, :]/AB[n_dim+1:end, :]
        ##########

        ## the rest is the same
        var_obs = Z*var_inter*Z' .+ E

        v_t = data_y[iter_data,:] .- Z*mean_out
        prob_out = prob_out - .5*n_dim*log(2*pi) -.5*log(abs(det(var_obs))) - 0.5* v_t'*(var_obs\v_t)
        mean_out = mean_out .+ var_inter*Z'*(var_obs\v_t)
        var_out = var_inter .- var_inter*Z'*(var_obs\Z)*var_inter

    end

    return prob_out
end

#=

function [mean_out, var_out, prob_out] = kalman_ct(F, H, Q, R, mean_0, var_0, data_y, dt)
  % Compute Kalman filter of a continuous time process sampled at discrete time
  %
  % .. math::
  %
  %   & dx_t = Fx_t dt + J dW_t \\
  %   & y_{t_k} = Hx_{t_k} + v_{t_k} \\
  %   & x_0 \sim N(mean_0, var_0) \\
  %   & k\in\{1, 2, \ldots, n\} \\
  %   & W_t = \text{Brownian Motion} \\
  %   & v_k \sim N(0, R) \\
  %   & dt_i = t_{i} - t_{i-1}
  %
  % Args:
  %   F (matrix): see equation above
  %   H (matrix): see equation above
  %   Q (matrix): :math:`Q = JJ'` i.e., variance
  %   R (matrix): see equation above
  %   mean_0 (vector): see equation above
  %   var_0 (matrix): see equation above
  %   data_y (vector): :math:`y_k` above
  %   dt (vector): time period between data measurement (see above)
  %
  % Note:
  %   - The data does not have to come with uniform step size.
  %   - There is a `kalman <https://www.mathworks.com/help/control/ref/kalman.html>`_ function in MATLAB as well, but it is under control system toolbox.
  %
  % Returns:
  %   : [mean_out, var_out, prob_out]
  %
  %   - **mean_out** (vector): :math:`E(x_{t_n}| y_{t_1},\ldots, y_{t_n})`
  %   - **var_out** (matrix): :math:`var(x_{t_n}| y_{t_1},\ldots, y_{t_n})`
  %   - **prob_out** (double): log-likelihood of the data
  %
  % Reference:
  %
  %   - `Wikipedia article on Kalman Filtering <https://en.wikipedia.org/wiki/Kalman_filter#Hybrid_Kalman_filter>`_
  %   - :cite:`grewal2011kalman`
  %
  % .. Syntax:
  %
  % ..   [mean_out, var_out, prob_out] = kalman_ct(F, H, Q, R, mean_0, var_0, data_y, dt)
  %

  % Get data size
  n_data = size(data_y, 2);
  if (length(dt) < n_data) & (length(n_data) == 1)
    dt = repmat(dt, n_data, 1);
  end
  n_dim = length(mean_0);

  prob_out = 0;
  mean_out = mean_0;
  var_out = var_0;

  % Step through data points
  for iter_data = 1:n_data
    mean_out = expm(dt(iter_data)*F)*mean_out;
    AB = expm(dt(iter_data)*[F, Q; H'*(R\H), -F])*[var_out; speye(n_dim)];
    var_inter = AB(1:n_dim, :)*inv(AB((n_dim+1):end, :));

    var_obs = H*var_inter*H' + R;
    v_t = data_y(:, iter_data) - H*mean_out;

    prob_out = prob_out - n_dim*log(2*pi)/2 ...
               -(log(abs(det(var_obs))) + v_t'*(var_obs\v_t))/2;

    mean_out = mean_out + var_inter*H'*(var_obs\v_t);
    var_out = var_inter - (var_inter*H')*(var_obs\H)*var_inter;
  end
end

=#


text:
const PROJ_ROOT = dirname(dirname(@__DIR__))

const DATA_DIR = joinpath(PROJ_ROOT, "data")

    const RAW_DATA_DIR = joinpath(DATA_DIR, "raw")
    const PROCESSED_DATA_DIR = joinpath(DATA_DIR, "processed")
    const FIGURES_DATA_DIR = joinpath(DATA_DIR, "figures")

for dir in [DATA_DIR, RAW_DATA_DIR,
            PROCESSED_DATA_DIR, FIGURES_DATA_DIR]
    if !isdir(dir)
        mkpath(dir)
        println("created $dir")
    end
end