text:
import Base: in, +, -, *, ==, convert, isempty, isequal, length, push!, sizeof, append!
import Base: lastindex, firstindex, getindex, setindex!, sort, sort!
import SeisIO: GeoLoc, PZResp, InstrumentResponse
export RawData, FFTData, CorrData, gpu, cpu, sort, sort!, append!
export isequal, isempty, length, +, ==

# Many thanks to SeisIO for inspiration for NoiseData objects

const datafields = (:x,:fft,:corr)

abstract type NoiseData end

"""
  RawData

A structure for raw ambient noise data.

## Fields: RawData
| **Field** | **Description** |
|:--------------------|:--------------------|
| :name       | Freeform channel names |
| :id         | Channel ids. use NET.STA.LOC.CHAN format when possible. |
| :loc        | Location (position) object |
| :fs         | Sampling frequency in Hz. |
| :gain       | Scalar gain; divide data by the gain to convert to units  |
| :freqmin    | Minimum frequency from instrument response.  |
| :freqmax    | Maximum frequency from instrument response. |
| :cc_len     | Raw data window length in seconds. |
| :cc_step    | Spacing between windows in seconds. |
| :time_norm  | Apply one-bit whitening with "one_bit". |
| :resp       | SeisIO InstrumentResponse object |
| :misc       | Dictionary for non-critical information. |
| :notes      | Timestamped notes; includes automatically-logged acquisition and |
|             | processing information. |
| :t          | Starttime of each window. |
| :x          | Raw data stored in columns. |
"""
mutable struct RawData <: NoiseData
   name::String                                # name [Net.Sta.Loc.Chan]
   id::String                                  # id [Y-mm-dd] this is date of raw
   loc::GeoLoc                                 # loc
   fs::Float64                                 # sampling rate [Hz]
   gain::Float64                               # gain
   freqmin::Float64                            # minumum frequency [Hz]
   freqmax::Float64                            # maximum frequency [Hz]
   cc_len::Float64                             # window_length [S]
   cc_step::Float64                            # step between windows [S]
   whitened::Bool                              # whitening applied
   time_norm::String                           # time normaliation
   resp::InstrumentResponse                    # response poles/zeros
   misc::Dict{String,Any}                      # misc
   notes::Array{String,1}                      # notes
   t::Array{Float64,1}                         # time
   x::AbstractArray{<:AbstractFloat,2}           # raw data

 function RawData(
     name     ::String,
     id       ::String,
     loc      ::GeoLoc,
     fs       ::Float64,
     gain     ::Float64,
     freqmin  ::Float64,
     freqmax  ::Float64,
     cc_len   ::Float64,
     cc_step  ::Float64,
     whitened ::Bool,
     time_norm::String,
     resp     ::InstrumentResponse,
     misc     ::Dict{String,Any},
     notes    ::Array{String,1},
     t        ::Array{Float64,1},
     x        ::AbstractArray{<:AbstractFloat,2}
     )

     return new(name, id, loc, fs, gain, freqmin, freqmax, cc_len, cc_step,
             whitened, time_norm, resp, misc, notes, t, x)
   end

   function RawData(S::SeisData,cc_len::Real,cc_step::Real)
     merge!(S)
     # check cc_len & fs compatibility
     if mod(cc_len * S.fs[1], 1.) !== 0.
          throw(DomainError(cc_len, "cc_len * fs must be an integer for slicing."))
     end

     # check cc_step and sampling rate
     if mod(cc_step * S.fs[1], 1.) !== 0.
          throw(DomainError(cc_step, "cc_step * fs must be an integer for slicing."))
     end
     ungap!(S)

     # convert input to float 64
     cc_len = Float64(cc_len)
     cc_step = Float64(cc_step)

     # check if waveform length is < cc_len
     if length(S.x[1]) < cc_len * S.fs[1]
         throw(DomainError(cc_len, "cc_len must be less than length of data."))
     end

     # phase shift data onto exaclty the sampling rate
     phase_shift!(S)
     # subset by time
     startslice, endslice = nearest_start_end(S[1],cc_len, cc_step)
     # get starting and ending indices
     startind, endind = slide_ind(startslice, endslice, S.fs[1], S.t[1])
     x, starts = slide(@view(S.x[1][startind:endind]), cc_len, cc_step, S.fs[1], startslice,endslice)
     return new(S[1].id,Dates.format(u2d(S.t[1][1,2]*1e-6),"Y-mm-dd"),S.loc[1],
                S.fs[1],S.gain[1],1. / cc_len,S[1].fs/2,cc_len,cc_step,false,
                "",S[1].resp,S.misc[1],S.notes[1],starts,x)
   end

   function RawData(C::SeisChannel,cc_len::Real,cc_step::Real)
    # check cc_len and sampling rate
    if mod(cc_len * C.fs, 1.) !== 0.
         throw(DomainError(cc_len, "cc_len * fs must be an integer for slicing."))
    end

    # check cc_step and sampling rate
    if mod(cc_step * C.fs, 1.) !== 0.
         throw(DomainError(cc_step, "cc_step * fs must be an integer for slicing."))
    end

    ungap!(C)
    # check if waveform length is < cc_len
    if length(C.x) < cc_len * C.fs
        throw(DomainError(cc_len, "cc_len must be less than length of data."))
    end

    # convert input to float 64
    cc_len = Float64(cc_len)
    cc_step = Float64(cc_step)

    # phase shift data onto exaclty the sampling rate
    phase_shift!(C)

    # subset by time
    startslice, endslice = nearest_start_end(C,cc_len, cc_step)

    # get starting and ending indices
    startind, endind = slide_ind(startslice, endslice, C.fs,C.t)
    x, starts = slide(@view(C.x[startind:endind]), cc_len, cc_step, C.fs, startslice,endslice)
    return new(C.id,Dates.format(u2d(C.t[1,2]*1e-6),"Y-mm-dd"),C.loc,C.fs,
               C.gain,1. / cc_len,C.fs/2,cc_len,cc_step,false,"", C.resp,
               C.misc,C.notes,starts,x)
   end
end

RawData(;
         name     ::String                    = "",
         id       ::String                    = "",
         loc      ::GeoLoc                    = GeoLoc(),
         fs       ::Float64                   = zero(Float64),
         gain     ::Float64                   = one(Float64),
         freqmin  ::Float64                   = zero(Float64),
         freqmax  ::Float64                   = zero(Float64),
         cc_len   ::Float64                   = zero(Float64),
         cc_step  ::Float64                   = zero(Float64),
         whitened ::Bool                      = false,
         time_norm::String                    = "",
         resp     ::InstrumentResponse        = PZResp(),
         misc     ::Dict{String,Any}          = Dict{String,Any}(),
         notes    ::Array{String,1}           = Array{String,1}(undef, 0),
         t        ::Array{Float64,1}          = Array{Float64,1}(undef, 0),
         x        ::AbstractArray{<:AbstractFloat,2} = Array{Float32,2}(undef, 0, 2)
         ) = RawData(name, id, loc, fs, gain, freqmin, freqmax, cc_len, cc_step,
                whitened, time_norm, resp, misc, notes, t, x)

"""
FFTData

A structure for fourier transforms (FFT) of ambient noise data.

## Fields: FFTData
| **Field** | **Description** |
|:--------------------|:--------------------|
| :name       | Freeform channel names |
| :id         | Channel ids. use NET.STA.LOC.CHAN format when possible. |
| :loc        | Location (position) object |
| :fs         | Sampling frequency in Hz. |
| :gain       | Scalar gain; divide data by the gain to convert to units  |
| :freqmin    | Minimum frequency for whitening.  |
| :freqmax    | Maximum frequency for whitening. |
| :cc_len     | Raw data window length in number of points. |
| :cc_step    | Spacing between windows in number of points. |
| :whitened   | Whitening applied.
| :time_norm  | Apply one-bit whitening with "one_bit". |
| :resp       | SeisIO InstrumentResponse object |
| :misc       | Dictionary for non-critical information. |
| :notes      | Timestamped notes; includes automatically-logged acquisition and |
|             | processing information. |
| :t          | Starttime of each FFT. |
| :fft        | FFTs stored in columns. |
"""
mutable struct FFTData <: NoiseData
   name::String                                # name [Net.Sta.Loc.Chan]
   id::String                                  # id [Y-mm-dd] this is date of fft
   loc::GeoLoc                                 # loc
   fs::Float64                                 # sampling rate [Hz]
   gain::Float64                               # gain
   freqmin::Float64                            # minumum frequency [Hz]
   freqmax::Float64                            # maximum frequency [Hz]
   cc_len::Float64                             # window_length [S]
   cc_step::Float64                            # step between windows [S]
   whitened::Bool                              # whitening applied
   time_norm::String                           # time normaliation
   resp::InstrumentResponse                    # response poles/zeros
   misc::Dict{String,Any}                      # misc
   notes::Array{String,1}                      # notes
   t::Array{Float64,1}                         # time
   fft::AbstractArray{<:Complex{<:AbstractFloat},2}  # fft data

function FFTData(
   name     ::String,
   id       ::String,
   loc      ::GeoLoc,
   fs       ::Float64,
   gain     ::Float64,
   freqmin  ::Float64,
   freqmax  ::Float64,
   cc_len   ::Float64,
   cc_step  ::Float64,
   whitened ::Bool,
   time_norm::String,
   resp     ::InstrumentResponse,
   misc     ::Dict{String,Any},
   notes    ::Array{String,1},
   t        ::Array{Float64,1},
   fft      ::AbstractArray{<:Complex{<:AbstractFloat},2}
   )

   return new(name, id, loc, fs, gain, freqmin, freqmax, cc_len, cc_step,
           whitened, time_norm, resp, misc, notes, t, fft)
 end
end

FFTData(;
       name     ::String                    = "",
       id       ::String                    = "",
       loc      ::GeoLoc                    = GeoLoc(),
       fs       ::Float64                   = zero(Float64),
       gain     ::Float64                   = one(Float64),
       freqmin  ::Float64                   = zero(Float64),
       freqmax  ::Float64                   = zero(Float64),
       cc_len   ::Float64                   = zero(Float64),
       cc_step  ::Float64                   = zero(Float64),
       whitened ::Bool                      = false,
       time_norm::String                    = "",
       resp     ::InstrumentResponse        = PZResp(),
       misc     ::Dict{String,Any}          = Dict{String,Any}(),
       notes    ::Array{String,1}           = Array{String,1}(undef, 0),
       t        ::Array{Float64,1}          = Array{Float64,1}(undef, 0),
       fft      ::AbstractArray{<:Complex{<:AbstractFloat},2} =
                  Array{Complex{Float32},2}(undef, 0, 2)
       ) = FFTData(name, id, loc, fs, gain, freqmin, freqmax, cc_len, cc_step,
             whitened, time_norm, resp, misc, notes, t, fft)

# FFTData(C::SeisChannel,freqmin::Float64, freqmax::Float64,cc_len::Int,
#      cc_step::Int, whitened::Bool,time_norm::Union{Bool,String},
#      t::Array{Float64,1},fft::Array{<:Union{Complex{Float32},Complex{Float64}},2}
#     ) = FFTData(C.id, Dates.format(u2d(C.t[1,2]*1e-6),"Y-mm-dd"), C.loc, C.fs, C.gain, freqmin, freqmax, cc_len,
#                 cc_step, whitened, time_norm, C.resp, C.misc, C.notes, t, fft)

"""
   CorrData

A structure for cross-correlations of ambient noise data.

## Fields: CorrData
| **Field** | **Description** |
|:------------|:-------------|
| :name       | Freeform channel name in NET1.STA1.LOC1.CHAN1.NET2.STA2.LOC2.CHAN2 form |
| :id         | Channel ids. use NET.STA.LOC.CHAN format when possible. |
| :loc        | Location (position) object for station 1. |
| :comp       | Component of cross-correlation (e.g ZZ, RT, etc..). |
| :rotated    | Boolean for rotation. |
| :fs         | Sampling frequency in Hz. |
| :gain       | Scalar gain; divide data by the gain to convert to units.  |
| :freqmin    | Minimum frequency for whitening.  |
| :freqmax    | Maximum frequency for whitening. |
| :cc_len     | Raw data window length in number of points. |
| :cc_step    | Spacing between windows in number of points. |
| :whitened   | Whitening applied.
| :time_norm  | Apply one-bit whitening with "one_bit". |
| :resp       | SeisIO InstrumentResponse object |
| :misc       | Dictionary for non-critical information. |
| :notes      | Timestamped notes; includes automatically-logged acquisition and |
|             | processing information. |
| :dist       | Distance between station 1 and station 2 in Km. |
| :azi        | Azimuth from station 1 and station 2 in degrees. |
| :baz        | Back azimuth between station 1 and station 2 in degrees. |
| :maxlag     | Maximum lag time in seconds to keep from correlations. |
| :t          | Starttime of each correlation. |
| :corr       | Correlations stored in columns. |
"""
mutable struct CorrData <: NoiseData
  name::String                                # name [Net1.Sta1.Loc1.Chan1.Net2.Sta2.Loc2.Chan2]
  id::String                                  # id [Y-mm-dd] this is date of corr
  loc::GeoLoc                                 # loc
  comp::String                                # 1st channel, 2nd channel [ZZ,RT,..]
  rotated::Bool                               # rotation applied?
  corr_type::String                           # "corr", "deconv", "coh"
  fs::Float64                                 # sampling rate [Hz]
  gain::Float64                               # gain
  freqmin::Float64                            # minumum frequency [Hz]
  freqmax::Float64                            # maximum frequency [Hz]
  cc_len::Float64                             # window_length [S]
  cc_step::Float64                            # step between windows [S]
  whitened::Bool                              # whitening applied
  time_norm::String                           # time normaliation
  resp::InstrumentResponse                    # response poles/zeros
  misc::Dict{String,Any}                      # misc
  notes::Array{String,1}                      # notes
  dist::Float64                               # distance between stations [km]
  azi::Float64                                # azimuth between stations [degrees]
  baz::Float64                                # back azimuth between stations [degrees]
  maxlag::Float64                             # maximum lag time [s]
  t::Array{Float64,1}                         # time
  corr::AbstractArray{<:AbstractFloat,2}     # corr data

  function CorrData(
      name     ::String,
      id       ::String,
      loc      ::GeoLoc,
      comp     ::String,
      rotated  ::Bool,
      corr_type::String,
      fs       ::Float64,
      gain     ::Float64,
      freqmin  ::Float64,
      freqmax  ::Float64,
      cc_len   ::Float64,
      cc_step  ::Float64,
      whitened ::Bool,
      time_norm::String,
      resp     ::InstrumentResponse,
      misc     ::Dict{String,Any},
      notes    ::Array{String,1},
      dist     ::Float64,
      azi      ::Float64,
      baz      ::Float64,
      maxlag   ::Float64,
      t        ::Array{Float64,1},
      corr     ::AbstractArray{<:AbstractFloat,2}
      )

      return new(name, id, loc, comp, rotated, corr_type, fs, gain, freqmin,
              freqmax, cc_len, cc_step, whitened, time_norm, resp, misc,
              notes, dist, azi, baz, maxlag, t, corr)
    end
end

CorrData(;
          name     ::String                    = "",
          id       ::String                    = "",
          loc      ::GeoLoc                    = GeoLoc(),
          comp     ::String                    = "",
          rotated  ::Bool                      = false,
          corr_type::String                    = "",
          fs       ::Float64                   = zero(Float64),
          gain     ::Float64                   = one(Float64),
          freqmin  ::Float64                   = zero(Float64),
          freqmax  ::Float64                   = zero(Float64),
          cc_len   ::Float64                   = zero(Float64),
          cc_step  ::Float64                   = zero(Float64),
          whitened ::Bool                      = false,
          time_norm::String                    = "",
          resp     ::InstrumentResponse        = PZResp(),
          misc     ::Dict{String,Any}          = Dict{String,Any}(),
          notes    ::Array{String,1}           = Array{String,1}(undef, 0),
          dist     ::Float64                   = zero(Float64),
          azi      ::Float64                   = zero(Float64),
          baz      ::Float64                   = zero(Float64),
          maxlag   ::Float64                   = zero(Float64),
          t        ::Array{Float64,1}          = Array{Float64,1}(undef, 0),
          corr     ::AbstractArray{<:AbstractFloat,2} = Array{Float32,2}(undef, 0, 2)
          ) = CorrData(name, id, loc, comp, rotated, corr_type, fs, gain,
                 freqmin, freqmax, cc_len, cc_step, whitened, time_norm,
                 resp, misc, notes, dist, azi, baz, maxlag, t, corr)

CorrData(F1::FFTData, F2::FFTData, comp::String, rotated::Bool, corr_type::String,
       maxlag::Float64, t::Array{Float64,1}, corr::AbstractArray{<:AbstractFloat,2}
      ) = CorrData(F1.name * '.' * F2.name, F1.id, F1.loc, comp, rotated,
                 corr_type, F1.fs, F1.gain, F1.freqmin, F1.freqmax, F1.cc_len,
                 F1.cc_step, F1.whitened, F1.time_norm, F1.resp, F1.misc,
                 F1.notes, get_dist(F1,F2), get_azi(F1,F2), get_baz(F1,F2),
                 maxlag, t, corr)

in(s::String, A::NoiseData) = A.id==s
length(N::NoiseData) = length(getfield(N,:t))
firstindex(N::NoiseData) = 1
lastindex(N::NoiseData) = length(N)

function getindex(N::NoiseData, J::Array{Int,1})
  n = length(N)
  T = typeof(N)
  U = T()
  F = fieldnames(T)
  for f in F
     if f in datafields
        setfield!(U, f, getindex(getfield(N,f),:,J))
     elseif f == :t
        setfield!(U, f, getindex(getfield(N,f),J))
     else
        setfield!(U, f, getfield(N,f))
     end
  end
  return U
end

function getindex(N::NoiseData, J::Array{DateTime,1})
   n = length(N)
   T = typeof(N)
   U = T()
   F = fieldnames(T)
   t = d2u.(J)
   ind = findall(in(t),N.t)
   for f in F
      if f in datafields
         setfield!(U, f, getindex(getfield(N,f),:,ind))
      elseif f == :t
         setfield!(U, f, getindex(getfield(N,f),ind))
      else
         setfield!(U, f, getfield(N,f))
      end
   end
   return U
 end
getindex(N::NoiseData, J::AbstractRange) = getindex(N, collect(J))

function setindex!(N::T, U::T, J::Array{Int,1}) where {T<:NoiseData}
  typeof(N) == typeof(U) || throw(MethodError)
  length(J) == length(U) || throw(BoundsError)
  length(N) >= maximum(J) || throw(BoundsError)
  F = fieldnames(T)
  for f in F
     if f in datafields
        setindex!(getfield(N, f), getfield(U, f), :,J)
     elseif f == :t
        setindex!(getfield(N,f),getfield(U,f),J)
     end
  end
  # ([(getfield(S, f))[J] = getfield(U, f) for f in datafields])
  return nothing
end
setindex!(N::NoiseData, U::NoiseData, J::AbstractRange) = setindex!(N, U, collect(J))

isempty(R::RawData) = min(R.name=="",R.id == "",isempty(R.loc),
                          R.fs == zero(typeof(R.fs)),
                          R.gain == one(typeof(R.gain)),
                          R.freqmin == zero(typeof(R.freqmin)),
                          R.freqmax == zero(typeof(R.freqmax)),
                          R.cc_len == zero(typeof(R.cc_len)),
                          R.cc_step == zero(typeof(R.cc_step)),
                          R.time_norm == "",
                          isempty(R.resp), isempty(R.misc),isempty(R.notes),
                          isempty(R.t),isempty(R.x))

isempty(F::FFTData) = min(F.name=="",F.id == "",isempty(F.loc),
                         F.fs == zero(typeof(F.fs)),
                         F.gain == one(typeof(F.gain)),
                         F.freqmin == zero(typeof(F.freqmin)),
                         F.freqmax == zero(typeof(F.freqmax)),
                         F.cc_len == zero(typeof(F.cc_len)),
                         F.cc_step == zero(typeof(F.cc_step)),
                         F.whitened == false,
                         F.time_norm == "",
                         isempty(F.resp), isempty(F.misc),isempty(F.notes),
                         isempty(F.t),isempty(F.fft))

isempty(C::CorrData) = min(C.name=="",C.id == "",isempty(C.loc), C.comp=="",
                        C.rotated==false,C.fs == zero(typeof(C.fs)),
                        C.corr_type=="",
                        C.gain == one(typeof(C.gain)),
                        C.freqmin == zero(typeof(C.freqmin)),
                        C.freqmax == zero(typeof(C.freqmax)),
                        C.cc_len == zero(typeof(C.cc_len)),
                        C.cc_step == zero(typeof(C.cc_step)),
                        C.whitened == false,
                        C.time_norm == "",
                        isempty(C.resp), isempty(C.misc),isempty(C.notes),
                        C.maxlag == zero(typeof(C.maxlag)), isempty(C.t),
                        isempty(C.corr))

function isequal(N::T,U::T) where {T<:NoiseData}
   return minimum([hash(getfield(N,i))==hash(getfield(U,i)) for i in fieldnames(T)]::Array{Bool,1})
end
==(N::T,U::T) where {T<:NoiseData} = isequal(N,U)

function append!(N::T, U::T) where {T<:NoiseData}
   F = fieldnames(T)
   if N.name == U.name
      d = intersect(datafields,F)[1]
      setfield!(N, :t, vcat(getfield(N,:t), getfield(U,:t)))
      setfield!(N, d, hcat(getfield(N,d), getfield(U,d)))
   elseif isempty(N)
      [setfield!(N, f, getfield(U,f)) for f in F]
   end
   return N
end

function sort!(N::T; rev=false::Bool) where {T<:NoiseData}
  j = sortperm(getfield(N, :t), rev=rev)
  F = fieldnames(T)
  d = intersect(datafields,F)[1]
  setfield!(N, d, getfield(N,d)[:,j])
  setfield!(N, :t, getfield(N,:t)[j])
  return nothing
end

function sort(N::T; rev=false::Bool) where {T<:NoiseData}
  U = deepcopy(N)
  sort!(U, rev=rev)
  return U
end

+(N::T, U::T) where {T<:NoiseData} = (Z = deepcopy(N); return append!(Z, U))

# Allow for transfer from CPU to GPU

Adapt.adapt_structure(to, R::RawData) = RawData(R.name, R.id, R.loc, R.fs, R.gain,
                                                R.freqmin, R.freqmax, R.cc_len,
                                                R.cc_step, R.whitened, R.time_norm,
                                                R.resp, R.misc, R.notes, R.t,
                                                adapt(to, R.x))
Adapt.adapt_structure(to, F::FFTData) = FFTData(F.name, F.id, F.loc, F.fs, F.gain,
                                                F.freqmin, F.freqmax, F.cc_len,
                                                F.cc_step, F.whitened, F.time_norm,
                                                F.resp, F.misc, F.notes, F.t,
                                                adapt(to, F.fft))
Adapt.adapt_structure(to, C::CorrData) = CorrData(C.name, C.id, C.loc, C.comp,
                                                  C.rotated, C.corr_type, C.fs,
                                                  C.gain, C.freqmin, C.freqmax,
                                                  C.cc_len, C.cc_step, C.whitened,
                                                  C.time_norm, C.resp, C.misc,
                                                  C.notes, C.dist,C.azi,C.baz,
                                                  C.maxlag,C.t,
                                                  adapt(to, C.corr))

# CPU/GPU movement conveniences
cpu(m::NoiseData) = adapt(Array,m)
gpu(x::NoiseData) = use_cuda[] ? cu(x) : x


text:
function quantnormal(
    obs::Vector{Float64}, mu::Float64, sigma::Float64, q::Float64, sample_size::Int,
    config::Config
)::NTuple{4, Vector{Float64}}

    weights, mus, sigmas = adaptive_parzen_normal(obs, mu, sigma, config)
    post = GMM.GMM1(weights, mus, sigmas, q, sample_size)

    return post, weights, mus, sigmas
end

ext: jl


text:
"""
This is a library of functions to help out in the Robotlab at LTH
The module imports the following files\n
include("utils.jl")\n
include("DH.jl")\n
include("kinematics.jl")\n
include("robotplot.jl")\n
include("Frames.jl")\n
include("read_log.jl")\n
include("csv2mat.jl")\n
include("posDepFric.jl")\n

The module includes a submodule, Frames, which is aimed at replacing the Nikon K600 software. It supports creation of frames, simple projections, fitting of planes, lines etc. and has a number of plotting options. It must be separately imported with `using Robotlib.Frames`

The module includes a submodule, Calibration, which includes a number of calibration routines. It must be separately imported with `using Robotlib.Calibration`

Usage:
fkine, ikine, jacobian = get_kinematic_functions("yumi")
data = orcalog2mat(pathopen, pathsave)
q = getData(\"robot_0.*posRawAbs\", data, 1, removeNaN = false)

For YuMi, joint angles `q` must be converted to logical order using e.g. abb2logical!(q)
You must also consider the base transform of YuMi
"""
module Robotlib
using LinearAlgebra, Statistics, StaticArrays, SparseArrays
const I4 = SMatrix{4, 4, Float64, 16}(Matrix{Float64}(I, 4, 4))
const I3 = SMatrix{3, 3, Float64, 9}(Matrix{Float64}(I, 3, 3))
include("DH.jl")
include("utils.jl")
include("kinematics.jl")
include("Frames.jl") # Must come before robotplot
include("robotplot.jl")
# https://github.com/JuliaIO/MAT.jl/issues/90
include("read_log.jl")# Re-enable when MAT.jl builds under Julia v1.0. Then add MAT to require
include("csv2mat.jl") # Re-enable when MAT.jl builds under Julia v1.0. Then add MAT to require
include("posDepFric.jl")

module Calibration
using TotalLeastSquares
using ..Robotlib
using ..Robotlib: T2R,Rt2T, T2R, T2t, skewcoords, twistcoords, I3, I4
using LinearAlgebra, Statistics, SparseArrays
import Robotlib.xii
import Robotlib.Ai
include("calibLPOE.jl")
include("calibForce.jl")
include("calibNAXP.jl")
export calibLPOE, calibLPOEdual, calibForce, calibPOE_offsets_from_points, calibNAXP

end

export skewcoords, twistcoords, skew, skew4, expω, expξ, expξ2, expξ!, logT, logR, trinv, isrot, isse3, Rangle, conformize, DH2twistsPOE, DH2twistsLPOE, dh2Tn, toOrthoNormal, toOrthoNormal!, rpy2R, Quaternion, xyθ, smartDiff, R2rpy, traj2quat, centraldiff
export fkinePOE, fkineLPOE, ikinePOE, dh2Tn, jacobianPOE, jacobianPOEikine, jacobian, get_kinematic_functions
export trajplot, trajplot3, plot3smart
export DH, DH2400, DHYuMi, DH7600, DHtest, abb2logical!, logical2abb!, abb2logical, logical2abb
export csv2mat, orcalog2mat, getData, readmat
export frictionRBFN, getCenters

dir(x...) = joinpath(dirname(pathof(Robotlib)),x...)

precompile(get_kinematic_functions, (String,))
precompile(fkineLPOE, (Array{Float64,3},Array{Float64,2},Array{Float64,1}))
precompile(jacobianPOE, (Array{Float64,1},Array{Float64,2}))
end


text:
#spherecomponent.jl
#Algorithm 2.3 of Random Eigenvalues by Alan Edelman

#Experiment:    Random components of a vector on a sphere
#Plot:          Histogram of a single component
#Theory:        Beta distribution

## Parameters
t = 100000      # trials
n = 12          # dimensions of sphere
dx = 0.1        # bin size

## Experiment
v = randn(n, t)
v = (v[1,:] ./ sqrt(sum(v.^2,1)))'
grid, count = hist(v, -1:dx:1)

## Theory
x = grid
y = (gamma(n-1)/(2^(n-2)*gamma((n-1)/2)^2))*(1-x.^2).^((n-3)/2)

## Plot
using Winston
p = FramedPlot()
h = Histogram(count/(t*dx), step(grid))
h.x0 = first(grid)
add(p, h)
add(p, Curve(x, y, "color", "blue", "linewidth", 2))
if isinteractive()
    Winston.display(p)
else
    file(p, "spherecomponent.png")
end


text:
using ITensors
using ITensorNetworkAD
using Test

using ITensorNetworkAD.ITensorNetworks: itensor_network, neighbors

@testset "itensor_network.jl" begin
  @testset "itensor_network from dims" begin
    d = (3, 3)
    tn = itensor_network(d...; linkdims=3)
    @test tn isa Matrix{ITensor}
    @test size(tn) == (3, 3)
    @test all(ITensors.isemptystorage, tn)
    @test hascommoninds(tn[1, 1], tn[1, 2])
    @test hascommoninds(tn[1, 1], tn[2, 1])
    @test !hascommoninds(tn[1, 1], tn[2, 2])
    @test hascommoninds(tn[1, 1], tn[3, 1])
    @test hascommoninds(tn[1, 1], tn[1, 3])
    @test isempty(uniqueinds(tn[2, 2], tn[2, 1], tn[2, 3], tn[1, 2], tn[3, 2]))
    @test issetequal(
      neighbors(tn, CartesianIndex(1, 2)),
      [
        CartesianIndex(1, 1),
        CartesianIndex(2, 2),
        CartesianIndex(3, 2),
        CartesianIndex(1, 3),
      ],
    )
  end
  @testset "itensor_network from siteinds" begin
    d = (3, 3)
    s = siteinds("S=1/2", d...)
    tn = itensor_network(s; linkdims=3)
    @test tn isa Matrix{ITensor}
    @test size(tn) == (3, 3)
    @test all(ITensors.isemptystorage, tn)
    @test hascommoninds(tn[1, 1], tn[1, 2])
    @test hascommoninds(tn[1, 1], tn[2, 1])
    @test !hascommoninds(tn[1, 1], tn[2, 2])
    @test hascommoninds(tn[1, 1], tn[3, 1])
    @test hascommoninds(tn[1, 1], tn[1, 3])
    @test uniqueinds(tn[2, 2], tn[2, 1], tn[2, 3], tn[1, 2], tn[3, 2]) == [s[2, 2]]
    @test issetequal(
      neighbors(tn, CartesianIndex(1, 2)),
      [
        CartesianIndex(1, 1),
        CartesianIndex(2, 2),
        CartesianIndex(3, 2),
        CartesianIndex(1, 3),
      ],
    )
  end
end


text:
module BCMFBPTests

using BinaryCommitteeMachineFBP
using Test
using GZip

const N = 21
const K = 5
const α = 0.1
const kw = Dict(:randfact=>0.1, :seed=>135, :max_iters=>1, :damping=>0.5, :quiet=>false);

patternsfile = "patterns.txt.gz"

fps = [Scoping(0:0.1:10, 21), PseudoReinforcement(0:0.01:1), FreeScoping([(atanh(√ρ), (2-ρ)/(1-ρ), Inf) for ρ in 0:0.01:1])]

# write your own tests here
function tst()
    errs, messages, patterns = focusingBP(N, K, α; kw...)
    @test errs == 0

    kw[:quiet] = true
    errs, messages, patterns = focusingBP(N, K, α; kw...)
    @test errs == 0
    errs, messages, patterns = focusingBP(N, K, patterns; kw...)
    @test errs == 0
    errs, messages, patterns = focusingBP(N, K, (patterns.X, float(patterns.output)); kw...)
    @test errs == 0
    errs, messages, patterns = focusingBP(N, K, patternsfile; kw...)
    @test errs == 0
    errs, messages, patterns = focusingBP(N, K, α; kw..., accuracy2=:accurate)
    @test errs == 0
    errs, messages, patterns = focusingBP(N, K, α; kw..., accuracy1=:none)
    @test errs == 0
    for fp in fps
        errs, messages, patterns = focusingBP(N, K, α; kw..., fprotocol=fp)
        @test errs == 0
    end
    errs, messages, patterns = focusingBP(N, K, α; kw..., messfmt=:plain)
    @test errs == 0
    errs, messages, patterns = focusingBP(N, K, α; kw..., initmess=messages)
    @test errs == 0

    f = tempname()
    isfile(f) && rm(f)
    try
        write_messages(f, messages)
        errs, messages, patterns = focusingBP(N, K, α; kw..., initmess=f)
        @test errs == 0
        m2 = read_messages(f, MagT64)
        errs, messages, patterns = focusingBP(N, K, α; kw..., initmess=m2)
        @test errs == 0
        m2 = read_messages(f, MagP64)
        errs, messages, patterns = focusingBP(N, K, α; kw..., initmess=m2)
        @test errs == 0
    finally
        isfile(f) && rm(f)
    end

    kw[:max_iters] = 100

    f = tempname()
    isfile(f) && rm(f)
    try
        errs, messages, patterns = focusingBP(N, K, α; kw..., fprotocol=PseudoReinforcement(0:0.01:0.99, 0.992:0.002:0.998),
                                              outatzero=false, outfile=f)
        @test errs == 0
        errs, messages, patterns = focusingBP(N, K, α; kw..., fprotocol=PseudoReinforcement(0:0.01:0.99, 0.992:0.002:0.998),
                                              outatzero=false, outfile=f, messfmt=:plain)
        @test errs == 0
    finally
        isfile(f) && rm(f)
    end

    f = tempname()
    isfile(f) && rm(f)
    ft = f * ".%gamma%"
    try
        errs, messages, patterns = focusingBP(N, K, α; kw..., fprotocol=PseudoReinforcement(0:0.01:0.99, 0.992:0.002:0.998),
                                              outatzero=false, outmessfiletmpl=ft)
        @test errs == 0
        errs, messages, patterns = focusingBP(N, K, α; kw..., fprotocol=PseudoReinforcement(0:0.01:0.99, 0.992:0.002:0.998),
                                              outatzero=false, outmessfiletmpl=ft, messfmt=:plain)
        @test errs == 0
    finally
        for fg in readdir(dirname(ft))
            startswith(fg, basename(f)) || continue
            dfg = joinpath(dirname(ft), fg)
            isfile(dfg) && rm(dfg)
        end
    end
end

@time tst()

end # module


text:
@testcase "#109: matching macrocall" begin
    @test 20 == @match raw"123" begin
        r"\Gaaa$" => 10
        raw"123" => 20
    end

    @lift @as_record struct Name
        n
    end

    @lift macro some_macro()
        n = esc(:Name)
        :(Name($n))
    end

    f(Name) = @match Name begin
        @some_macro() => Name
    end

    @test f(Name(1)) == 1
end


text:
struct CartPoleParams{T} <:Any where {T<:Real}
    m::T # mass of pole
    l::T # length of pole
    mc::T # mass of cart
    muc::T # some constant?
    mup::T # some constant?
    fmag::T # magnitude of force applied
	g::T  # gravity (not signed)
    CartPoleParams() = new{Float64}(0.1, 0.5, 1., 0.0005, 0.000002, 10., 9.8)
    CartPoleParams(T::Type) = new{T}(T(0.1), T(0.5), T(1.), T(0.0005), T(0.000002), T(10.), T(9.8))
	CartPoleParams(T::Type, m, l, mc, g) = new{T}(m, l, mc, T(0.0005), T(0.000002), T(10.), g)
end

function cartpole_sim(state, constants::CartPoleParams, u, dt)
	x, xDot, theta, thetaDot = view(state, 1:4)
	omegaDot = 0.
	vDot = 0.
	m = constants.m     # mass of pole
    l = constants.l     # length of pole
    mc = constants.mc   # mass of cart
    muc = constants.muc # some constant?
    mup = constants.mup # some constant?
	g = constants.g
	# for i in 1:steps
	omegaDot = (g * sin(theta) + cos(theta) * (muc * sign(xDot) - u - m * l * thetaDot^2 * sin(theta)) / (m + mc) - mup * thetaDot / (m * l)) / (l * (4.0 / 3.0 - m / (m + mc) * cos(theta)^2))
    vDot = (u + m * l * (thetaDot^2 * sin(theta) - omegaDot*cos(theta)) - muc*sign(xDot)) / (m + mc)
    theta += dt * thetaDot
    thetaDot += dt * omegaDot
    x += dt * xDot
    xDot += dt * vDot

    theta = mod(theta + π, 2 * π) - π

	return x, xDot, theta, thetaDot
end

function cartpole_compute_torque(p::CartPoleParams, action::Int)
	if action <= 0 || action > 2
		error("Action needs to be an integer in [1, 2]")
	end
	u = 0.
	if action == 1
		u = -p.fmag
	else
		u =  p.fmag
	end
	return u
end

function cartpole_compute_torque(p::CartPoleParams, action)
	return clamp(action, -1.0, 1.0) * p.fmag
end

function cartpole_step(state, action, params, dt, maxT)
	t,x = state
	u = cartpole_compute_torque(params, action)
	x .= cartpole_sim(x, params, u, dt)

	t += dt
	γ = 1.0
	if cartpole_terminal(x) || (t > maxT - 1e-8)
		fill!(x, 0.0)
		t = 0.0
		γ = 0.0
	end
	r = 1.0
	return t,x,r, γ
end

function cartpole_finitetime(;randomize=false, tMax=20.0, dt=0.02, Atype=:Discrete, droptime=true)
	if randomize
		params = random_cartpole_params()
	else
		params = CartPoleParams()
	end
	return create_finitetime_cartpole(params; tMax=tMax, dt=dt, Atype=Atype, droptime=droptime)
end

function create_finitetime_cartpole(params::CartPoleParams; tMax=20.0, dt=0.02, Atype=:Discrete, droptime=true)
	S = ([0. tMax],				# time range
		[	-2.4 2.4 ;       	# x range
			-10. 10. ;       	# xDot range
			-π/12.0 π/12.0 ; 	# theta range
			-π π ])           	# thetaDot range

	if Atype==:Discrete
		A = 1:2
	else
		A = [-1.0 1.0]
	end

	if droptime
		X = S[2]
		function get_outcome1(s,a,params,dt,tMax)
			t, x, r, γ = cartpole_step(s,a, params, dt, tMax)
			s = (t,x)
			return s,x,r,γ
		end
		p = (s,a)->get_outcome1(s,a,params,dt,tMax)
		function sample_initialdrop()
			x = zeros(4)
			return (0.0, x), x
		end
		d0 = sample_initialdrop
	else
		X = S
		function get_outcome2(s,a,params,dt,tMax)
			t, x, r, γ = cartpole_step(s,a, params, dt, tMax)
			s = (t,x)
			return s,s,r,γ
		end
		p = (s,a)->get_outcome2(s,a,params,dt,tMax)
		function sample_initial()
			x = zeros(4)
			return (0.0, x), (0.0, x)
		end
		d0 = sample_initial
	end

	meta = Dict{Symbol,Any}()
    meta[:minreward] = 1.0
    meta[:maxreward] = 1.0
    meta[:minreturn] = 9.0
    meta[:maxreturn] = ceil(tMax / dt)
    meta[:stochastic] = false
    meta[:minhorizon] = 9.0
    meta[:maxhorizon] = ceil(tMax / dt)
	meta[:discounted] = false
	meta[:episodes] = 100
	meta[:threshold] = 900

	render = (state,clearplot=false)->cartpoleplot(state, params)

	m = SequentialProblem(S,X,A,p,d0,meta,render)
	return m
end

function random_cartpole_params(::Type{T}=Float64) where {T}
	m = rand(Uniform(T(0.025), T(.25)))  # mass of pole
	l = rand(Uniform(T(0.1), T(1.0)))    # length of pole
	mc = rand(Uniform(T(0.1), T(5.)))    # mass of cart
	g = rand(Uniform(T(8.), T(10.)))     # gravity
	params = CartPoleParams(T, m, l, mc, g)
	return params
end

function cartpole_terminal(state)
	polecond = abs(state[3]) > (π / 15.0)
	cartcond = abs(state[1]) ≥ 2.4
	done = polecond | cartcond
	return done
end

# CartPole Rendering originally from JuliaML/Reinforce.jl and has only be modified slighlty
# https://github.com/JuliaML/Reinforce.jl/blob/master/src/envs/cartpole.jl
@userplot CartPolePlot
@recipe function f(cp::CartPolePlot)
	state, params = cp.args
	t, (x, xvel, θ, θvel) = state
	legend := false
	xlims := (-2.4, 2.4)
	l = params.l
	ylims := (-Inf, 2l)
	grid := false
	ticks := nothing

	# pole
	@series begin
		linecolor := :red
		linewidth := 10
		[x, x + 2l * sin(θ)], [0.0, 2l * cos(θ)]
	end

	# cart
	@series begin
		seriescolor := :black
		seriestype := :shape
		hw = 0.5
		l, r = x-hw, x+hw
		t, b = 0.0, -0.1
		[l, r, r, l], [t, t, b, b]
	end
end


text:
function SeisSmoothGathers(in, out,adj;Nsmooth=3,Nrepeat=1)

	if (adj==false)
		ngather = size(in,3)
		for igather = 1 : ngather
			A = in[:,:,igather];
			A = smooth_angles(A;Nsmooth=Nsmooth,Nrepeat=Nrepeat)
			out[:,:,igather] = A
		end
	else
		ngather = size(out,3)
		for igather = 1 : ngather
			A = out[:,:,igather];
			A = smooth_angles(A;Nsmooth=Nsmooth,Nrepeat=Nrepeat)
			in[:,:,igather] = A
		end
	end
return in,out
end

function smooth_angles(A;Nsmooth=3,Nrepeat=1)

	for iz = 1 : size(A,1)
		a = A[iz,:]
		for irepeat = 1 : Nrepeat
			a = mean_filter(a,Nsmooth)
		end
		A[iz,:] = a
	end

	return A;
end

function mean_filter(a,nw)

	b = vec(zeros(length(a),1))
	for ix = 1 : length(a)
		sum  = 0.
		nsum = 0.
		for iw = 1 : nw
			index1 = ix - Int(floor(nw/2)) - 1 + iw
			if (index1>0 && index1<length(a))
				sum += a[index1]
				nsum += 1.
			end
		end
		b[ix] = sum/nsum
	end
	return b
end


text:
module PoolOpt

export Haverly_1,Haverly_2,Haverly_3,BenTal_4,BenTal_5,Adhya_1,Adhya_2,Adhya_3,Adhya_4,Foulds_2,Foulds_3
export Foulds_4,DeyGupte_4,RT2,sppA0,sppA1,sppA2,sppA3,sppA4,sppA5,sppA6,sppB0
export elimination_equality, with_equality,pooling_with_eq_BSOS,pooling_without_eq_BSOS
export pooling_with_eq_Sparse_BSOS,pooling_without_eq_Sparse_BSOS,pooling_with_eq_Merge_Sparse_BSOS
export pooling_without_eq_Merge_Sparse_BSOS
include("pooling_data.jl")
include("pooling_problem_hierarchy.jl")

end
# module

max_stars_repo_name: shudianzhao/PoolOpt

max_issues_repo_name: shudianzhao/PoolOpt

max_forks_repo_name: shudianzhao/PoolOpt


text:
###################################
# Agent demographic profile generator
###################################

"""
Demographic profile generator

Creates socio-demographic profile of an agent based on demostats distributions per dissemination area (DA)

**Arguments**

* `DA_home` : unique id of DA centroid home location selected for an agent
* `DA_demostat` : dictionary  with population statistics for specific DA
* `demographic_categories` : dictionary with demographic categories used in generating of agent's profile

"""
function demographic_profile(DA_home::Int, DA_demostat::Nothing;
							demographic_categories = demographic_categories)::DataFrames.DataFrame
	profile = Dict()
	profile[:DA_home] = DA_home
	profile[:DA_work] = 0
	for key in keys(demographic_categories)
		categories = demographic_categories[key]
		category = rand(keys(categories))
		value = categories[category]
		if isa(value,UnitRange)
			profile[key] = rand(value)
		else
			profile[key] = value
		end
	end
	return DataFrames.DataFrame(profile)
end


text:
export zoo_jpc, zoo_jpc_file,
    JpcFormat, JpcVtreeFormat

struct JpcFormat <: FileFormat end

const JpcVtreeFormat = Tuple{JpcFormat,VtreeFormat}
Tuple{JpcFormat,VtreeFormat}() = (JpcFormat(),VtreeFormat())

##############################################
# Read JPC (Juice Probabilistic Circuit)
##############################################

zoo_jpc_file(name) =
    artifact"circuit_model_zoo" * zoo_version * "/jpcs/$name"

"""
    zoo_jpc(name)

Loads JPC file with given name from model zoo. See https://github.com/UCLA-StarAI/Circuit-Model-Zoo.
"""
zoo_jpc(name) =
    read(zoo_jpc_file(name), ProbCircuit, JpcFormat())

const jpc_grammar = raw"""
    start: header (_NL node)+ _NL?

    header : "jpc" _WS INT

    node : "L" _WS INT _WS INT _WS SIGNED_INT -> literal_node
         | "P" _WS INT _WS INT _WS INT child_nodes -> prod_node
         | "S" _WS INT _WS INT _WS INT weighted_child_nodes -> sum_node

    child_nodes : (_WS INT)+
    weighted_child_nodes: (_WS INT _WS LOGPROB)+

    %import common.INT
    %import common.SIGNED_INT
    %import common.SIGNED_NUMBER -> LOGPROB
    %import common.WS_INLINE -> _WS
    %import common.NEWLINE -> _NL
    """ * dimacs_comments

const jpc_parser = Lark(jpc_grammar)

abstract type JpcParse <: JuiceTransformer end

@inline_rule header(t::JpcParse, x) =
    Base.parse(Int,x)

@rule start(t::JpcParse, x) = begin
    @assert num_nodes(x[end]) == x[1]
    x[end]
end

@rule child_nodes(t::JpcParse, x) =
    map(id -> t.nodes[id], x)

@rule weighted_child_nodes(t::JpcParse, x) = begin
    children = map(id -> t.nodes[id], x[1:2:end])
    log_probs = Base.parse.(Float64,x[2:2:end])
    (children, log_probs)
end

#  parse unstructured
struct PlainJpcParse <: JpcParse
    nodes::Dict{String,PlainProbCircuit}
    PlainJpcParse() = new(Dict{String,PlainProbCircuit}())
end

@rule literal_node(t::PlainJpcParse, x) =
    t.nodes[x[1]] = PlainProbLiteralNode(Base.parse(Lit,x[3]))

@rule prod_node(t::PlainJpcParse,x) = begin
    @assert length(x[4]) == Base.parse(Int,x[3])
    t.nodes[x[1]] = PlainMulNode(x[4])
end

@rule sum_node(t::PlainJpcParse,x) = begin
    @assert length(x[4][1]) == length(x[4][2]) == Base.parse(Int,x[3])
    t.nodes[x[1]] = PlainSumNode(x[4][1], x[4][2])
end

function Base.parse(::Type{PlainProbCircuit}, str, ::JpcFormat)
    ast = Lerche.parse(jpc_parser, str)
    Lerche.transform(PlainJpcParse(), ast)
end

Base.read(io::IO, ::Type{PlainProbCircuit}, ::JpcFormat) =
    parse(PlainProbCircuit, read(io, String), JpcFormat())

#  parse structured
struct StructJpcParse <: JpcParse
    id2vtree::Dict{String,<:Vtree}
    nodes::Dict{String,StructProbCircuit}
    StructJpcParse(id2vtree) =
        new(id2vtree,Dict{String,StructProbCircuit}())
end

@rule literal_node(t::StructJpcParse, x) = begin
    lit = Base.parse(Lit,x[3])
    vtree = t.id2vtree[x[2]]
    t.nodes[x[1]] = StructProbLiteralNode(lit, vtree)
end

@rule prod_node(t::StructJpcParse,x) = begin
    @assert length(x[4]) == Base.parse(Int,x[3]) == 2
    vtree = t.id2vtree[x[2]]
    t.nodes[x[1]] = StructMulNode(x[4]..., vtree)
end

@rule sum_node(t::StructJpcParse,x) = begin
    @assert length(x[4][1]) == length(x[4][2]) == Base.parse(Int,x[3])
    vtree = t.id2vtree[x[2]]
    t.nodes[x[1]] = StructSumNode(x[4][1], x[4][2], vtree)
end

function Base.parse(::Type{StructProbCircuit}, str::AbstractString, ::JpcFormat, id2vtree)
    ast = Lerche.parse(jpc_parser, str)
    Lerche.transform(StructJpcParse(id2vtree), ast)
end

function Base.parse(::Type{StructProbCircuit}, strings, format::JpcVtreeFormat)
    id2vtree = parse(Dict{String,Vtree}, strings[2], format[2])
    parse(StructProbCircuit, strings[1], format[1], id2vtree)
end

Base.read(io::IO, ::Type{StructProbCircuit}, ::JpcFormat, id2vtree) =
    parse(StructProbCircuit, read(io, String), JpcFormat(), id2vtree)

function Base.read(ios::Tuple{IO,IO}, ::Type{StructProbCircuit}, ::JpcVtreeFormat)
    circuit_str = read(ios[1], String)
    vtree_str = read(ios[2], String)
    parse(StructProbCircuit, (circuit_str,vtree_str), JpcVtreeFormat())
end

##############################################
# Write JPCs
##############################################

const JPC_FORMAT = """c this file was saved by ProbabilisticCircuits.jl
c ids of jpc nodes start at 0
c jpc nodes appear bottom-up, children before parents
c
c file syntax:
c jpc count-of-jpc-nodes
c L id-of-literal-jpc-node id-of-vtree literal
c P id-of-sum-jpc-node id-of-vtree number-of-children {child-id}+
c S id-of-product-jpc-node id-of-vtree number-of-children {child-id}+ {log-probability}+
c"""

function Base.write(io::IO, circuit::ProbCircuit, ::JpcFormat, vtreeid::Function = (x -> 0))

    labeling = label_nodes(circuit)
    map!(x -> x-1, values(labeling)) # vtree nodes are 0-based indexed

    println(io, JPC_FORMAT)
    println(io, "jpc $(num_nodes(circuit))")
    foreach(circuit) do n
        if isliteralgate(n)
            println(io, "L $(labeling[n]) $(vtreeid(n)) $(literal(n))")
        elseif isconstantgate(n)
            @assert false
        else
            t = is⋀gate(n) ? "P" : "S"
            print(io, "$t $(labeling[n]) $(vtreeid(n)) $(num_children(n))")
            if is⋀gate(n)
                for child in children(n)
                    print(io, " $(labeling[child])")
                end
            else
                @assert is⋁gate(n)
                for (child, logp) in zip(children(n), n.log_probs)
                    print(io, " $(labeling[child]) $logp")
                end
            end
            println(io)
        end
    end
    nothing
end

function Base.write(ios::Tuple{IO,IO}, circuit::StructProbCircuit, format::JpcVtreeFormat)
    vtree2id = write(ios[2], vtree(circuit), format[2])
    write(ios[1], circuit, format[1], n -> vtree2id[vtree(n)])
end


text:
using Literate

function replace_includes(str)
    included = ["lv_model.jl", "lv_symbolic.jl"]
    path = @__DIR__

    for ex in included
        content = read(joinpath(path, ex), String)
        str = replace(str, "include(\"$(ex)\")" => content)
    end

    return str
end

Literate.markdown(joinpath(@__DIR__, "lv_comp.jl"), joinpath(@__DIR__, "docs/markdown/"), preprocess = replace_includes)
Literate.notebook(joinpath(@__DIR__, "lv_comp.jl"), joinpath(@__DIR__, "docs/notebook/"), preprocess = replace_includes)

max_stars_repo_name: McCannLab/TheoryClub

max_issues_repo_name: McCannLab/TheoryClub

max_forks_repo_name: McCannLab/TheoryClub


text:
using FileIO
using LinearAlgebra
using Parameters
include("network_sims.jl")
include("random_matrices.jl")

output_dir = "results2-analysis"

N = 3000
g = 1.5
etas = [-0.3, 0.0, 0.3]
s = 0.0
J_seed = 5000
h0_seed = 8000

dt = 0.01
t_bound = 300.0
wait = 200.0
savedunits = 10
timesteps = round(Int64, (t_bound - wait) / dt)
ts = dt .* collect(1:timesteps)

activity = Array{Float64,3}(undef, timesteps, savedunits, length(etas))
eigval_parts = Array{Float64,3}(undef, N, 2, length(etas))

for (i, eta) in enumerate(etas)
    J = generate_J(N, eta, seed=J_seed)
    h0 = generate_nv(N, Float64, seed=h0_seed)
    A = sim_activity(J, g, h0, s=s, dt=dt, t_bound=t_bound)
    activity[:,:,i] = view(A, 1:savedunits,
            round(Int, wait / dt) + 1:round(Int, t_bound / dt))'

    eigval_parts[:,1,i] = real(eigvals(J))
    eigval_parts[:,2,i] = imag(eigvals(J))
end

fn = "example_activity_and_eigvals_results2.jld2"
out = Dict{String,Any}()
@pack! out = ts, activity, eigval_parts
save(joinpath(output_dir, fn), out)


text:
using Experimental
using Documenter

makedocs(;
    modules=[Experimental],
    authors="Eric Ford",
    repo="https://github.com/RvSpectML/Experimental.jl/blob/{commit}{path}#L{line}",
    sitename="Experimental.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://RvSpectML.github.io/Experimental.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/RvSpectML/Experimental.jl",
)


text:
function cfg_net(
    ;
    agent_count::Int64=10000,
    m0::Int64=100,
    new_follows::Int64=10
)
    return (
        agent_count=agent_count,
        m0=m0,
        new_follows=new_follows
    )
end

function cfg_sim(
    ;
    ticks::Int64=100,
    addfriends::String=""
    )

    return(
        ticks=ticks,
        addfriends=addfriends
    )
end

function cfg_ot(
    ;
    backfire::Float64=0.4,
    befriend::Float64=0.2,
    unfriend::Float64=0.6
)
    return (
        backfire=backfire,
        befriend=befriend,
        unfriend=unfriend
    )
end

function cfg_ag(
    ;
    own_opinion_weight::Float64=0.90,
    unfriend_rate::Float64=0.05,
    min_friends_count::Int64=5
    )

    return (
        own_opinion_weight=own_opinion_weight,
        unfriend_rate=unfriend_rate,
        min_friends_count=min_friends_count
    )
end

struct Config
    network::NamedTuple{
        (:agent_count, :m0, :new_follows),
        NTuple{3,Int64}
    }
    simulation::NamedTuple{
    (:ticks, :addfriends),
    <:Tuple{Int64, String}
    }
    opinion_threshs::NamedTuple{
        (:backfire, :befriend, :unfriend),
        NTuple{3,Float64}
    }
    agent_props::NamedTuple{
    (:own_opinion_weight, :unfriend_rate, :min_friends_count),
    <:Tuple{Float64, Float64, Int64}
    }

    # constructor
    function Config(
        ;
        network = cfg_net(),
        simulation = cfg_sim(),
        opinion_threshs = cfg_ot(),
        agent_props = cfg_ag()
    )
        new(network, simulation, opinion_threshs, agent_props)
    end
end


text:
function WeightingOp(in,adj;w=1)

	return in.*w

end

function WeightingOp(m::String,d::String,adj;w="NULL")

	if (adj==true)
		d1,h1,e1 = SeisRead(d)
		d2,h2,e2 = SeisRead(w)
		SeisWrite(m,d1[:,:].*d2[:,:],h1,e1)
	else
		d1,h1,e1 = SeisRead(m)
		d2,h2,e2 = SeisRead(w)
		SeisWrite(d,d1[:,:].*d2[:,:],h1,e1)
	end

end

function WeightingOp(m::Array{String,1},d::Array{String,1},adj;w="NULL")

	if !(isdefined(w,2))
		for j = 1 : length(m)
			WeightingOp(m[j],d[j],adj,w=w)
		end
	else
		for j = 1 : length(m)
			WeightingOp(m[j],d[j],adj,w=w[j])
		end
	end
end


text:
const OneToZero = Base.OneTo(0)

for A in (:Missing, :Nothing, :Float64, :Float32, :Float16)
  @eval begin
    function locf(::Type{$A}, data::AbstractArray{T,1}) where {T<:IntFloat}
        indices = findall($A, data)
        return if isempty(indicies)
                   data
               else
                   locf($A, data, indicies)
               end
    end

    function locf(::Type{$A}, data::AbstractArray{T,1}, indicies) where {T}
        # cannot carry forward into data[bgn]
        if indicies[1] === 1
            indicies = length(indicies) === 1 ? OneToZero : indicies[2:end]
        end

        while !isempty(indices)
           prior_indicies = indicies - 1
           data[indices] = data[prior_indicies]
           indices = findall($A, data)
        end
        return data
    end

    function locf(::Type{$A}, data::AbstractArray{T,2}) where {T<:IntFloat}
        axs1, axs2 = axes(data)
        datavec = Vector{T}(undef, axs1.stop)
        for ax in axs2
            datavec[:] = data[axs1, ax]
            isnothing(findfirst($A, datavec)) && continue
            data[axs1, ax] = locf($A, datavec)
        end
        return data
    end

    function locf(::Type{$A}, data::AbstractArray{T,3}) where {T<:IntFloat}
        axs1, axs2, axs3 = axes(data)
        datamatrix = Matrix{T}(undef, axs1.stop, axs2.stop)
        for ax in axs3
            datamatrix[:,:] = data[axs1, axs2, ax]
            data[axs1, axs2, ax] = locf($A, datamatrix)
        end
        return data
    end

    function locf(::Type{$A}, data::AbstractArray{T,4}) where {T<:IntFloat}
        axs1, axs2, axs3, axs4 = axes(data)
        dataarray = Array{T,3}(undef, axs1.stop, axs2.stop, axs3.stop)
        for ax in axs4
            dataarray[:,:,:] = dataarray[axs1, axs2, axs3, ax]
            data[axs1, axs2, axs3, ax] = locf($A, dataarray)
        end
        return data
    end

    function locf(::Type{$A}, data::AbstractArray{T,5}) where {T<:IntFloat}
        axs1, axs2, axs3, axs4, axs5 = axes(data)
        dataarray = Array{T,4}(undef, axs1.stop, axs2.stop, axs3.stop, axs4.stop)
        for ax in axs5
            dataarray[:,:,:,:] = dataarray[axs1, axs2, axs3, axs4, ax]
            data[axs1, axs2, axs3, axs4, ax] = locf($A, dataarray)
        end
        return data
    end

    function locf(::Type{$A}, data::AbstractArray{T,6}) where {T<:IntFloat}
        axs1, axs2, axs3, axs4, axs5, axs6 = axes(data)
        dataarray = Array{T,5}(undef, axs1.stop, axs2.stop, axs3.stop, axs4.stop, axs5.stop)
        for ax in axs6
            dataarray[:,:,:,:,:] = dataarray[axs1, axs2, axs3, axs4, axs5, ax]
            data[axs1, axs2, axs3, axs4, axs5, ax] = locf($A, dataarray)
        end
        return data
    end

  end
end


text:
"""
    khatri_rao(A, B)

References
* Liu, Shuangzhe, and Gõtz Trenkler (2008) Hadamard, Khatri-Rao, Kronecker and Other Matrix Products. International J. Information and Systems Sciences 4, 160–177.
* Khatri, C. G., and Rao, C. Radhakrishna (1968) Solutions to Some Functional Equations and Their Applications to Characterization of Probability Distributions. Sankhya: Indian J. Statistics, Series A 30, 167–180.
"""
function khatri_rao(A::AbstractBlockMatrix, B::AbstractBlockMatrix)
    #
    Ablksize = blocksize(A)
    Bblksize = blocksize(B)

    @assert Ablksize == Bblksize "A and B must have the same blocksize"

    kblk = []
    for iblk in blockaxes(A,1)
        kblk_j = []
        for _jblk in blockaxes(A,2)
            Ablk = A[iblk, _jblk]
            Bblk = B[iblk, _jblk]
            push!(kblk_j, kron(Ablk, Bblk))
        end
        push!(kblk, tuple(kblk_j...))
    end
    mortar(kblk...)
end

function khatri_rao(A::AbstractMatrix, B::AbstractMatrix)
    kron(A, B)
end

""""
    BlockKron(A...)

creates a lazy representation of kron(A...) with the natural
block-structure imposed. This is a component in `blockkron(A...)`.
"""
struct BlockKron{T,N,ARGS<:Tuple} <: AbstractBlockArray{T,N}
    args::ARGS
end

BlockKron{T,N}(A...) where {T,N} = BlockKron{T,N,typeof(A)}(A)
BlockKron{T}(A::AbstractVector, B::AbstractVector, C::AbstractVector...) where {T} = BlockKron{T,1}(A, B, C...)
BlockKron{T}(A, B, C...) where {T} = BlockKron{T,2}(A, B, C...)
BlockKron(A, B, C...) = BlockKron{mapreduce(eltype,promote_type,(A,B,C...))}(A, B, C...)

size(B::BlockKron) = size(Kron(B))

size(K::BlockKron, j::Int) = prod(size.(K.args, j))
size(a::BlockKron{<:Any,1}) = (size(a,1),)
size(a::BlockKron{<:Any,2}) = (size(a,1), size(a,2))

function axes(K::BlockKron{<:Any,1})
    A,B = K.args
    (blockedrange(fill(prod(size.(tail(K.args),1)), size(K.args[1],1))),)
end

function axes(K::BlockKron{<:Any,2})
    A,B = K.args
    blockedrange.((fill(prod(size.(tail(K.args),1)), size(K.args[1],1)),
                   fill(prod(size.(tail(K.args),2)), size(K.args[1],2))))
end

kron_getindex((A,)::Tuple{AbstractVector}, k::Integer) = A[k]
function kron_getindex((A,B)::NTuple{2,AbstractVector}, k::Integer)
    K,κ = divrem(k-1, length(B))
    A[K+1]*B[κ+1]
end
kron_getindex((A,)::Tuple{AbstractMatrix}, k::Integer, j::Integer) = A[k,j]
function kron_getindex((A,B)::NTuple{2,AbstractVecOrMat}, k::Integer, j::Integer)
    K,κ = divrem(k-1, size(B,1))
    J,ξ = divrem(j-1, size(B,2))
    A[K+1,J+1]*B[κ+1,ξ+1]
end

kron_getindex(args::Tuple, k::Integer, j::Integer) = kron_getindex(tuple(BlockKron(args[1:2]...), args[3:end]...), k, j)
kron_getindex(args::Tuple, k::Integer) = kron_getindex(tuple(BlockKron(args[1:2]...), args[3:end]...), k)

getindex(K::BlockKron{<:Any,1}, k::Integer) = kron_getindex(K.args, k)
getindex(K::BlockKron{<:Any,2}, k::Integer, j::Integer) = kron_getindex(K.args, k, j)

kron_viewblock((a,b)::Tuple{Any,Any}, k::Integer) = a[k]*b
kron_viewblock(args, k::Integer) = args[1][k]*BlockKron(tail(args)...)

kron_viewblock((a,b)::Tuple{Any,Any}, k::Integer, j::Integer) = a[k,j]*b
kron_viewblock(args, k::Integer, j::Integer) = args[1][k,j]*BlockKron(tail(args)...)

viewblock(K::BlockKron{<:Any,1}, k::Block{1}) = kron_viewblock(K.args, Int(k))
viewblock(K::BlockKron{<:Any,2}, kj::Block{2}) = kron_viewblock(K.args, Int.(kj.n)...)

# const SubKron{T,M1,M2,R1,R2} = SubArray{T,2,<:BlockKron{T,M1,M2},<:Tuple{<:BlockSlice{R1},<:BlockSlice{R2}}}

# BroadcastStyle(::Type{<:SubKron{<:Any,<:Any,B,Block1,Block1}}) where B =
#     BroadcastStyle(B)

# allow dispatch on memory layout
_blockkron(_, A) = BlockArray(BlockKron(A...))

""""
    blockkron(A...)

creates a blocked version of kron(A...) with the natural
block-structure imposed.
"""
blockkron(A...) = _blockkron(map(MemoryLayout,A), A)

"""
    blockvec(A::AbstractMatrix)

creates a blocked version of `vec(A)`, with the block structure used to represent the columns.
"""
blockvec(A::AbstractMatrix) = PseudoBlockVector(vec(A), Fill(size(A,1), size(A,2)))


text:
using LinearAlgebra

include("diagonal_sbp.jl")

let
  N = 20
  λ = [mod(i, 2) + 1 for i = 0:N]
  for p in (2, 4, 6)
    (_, _, _, _, _, A, _) = variable_diagonal_sbp_D2(p, N, λ; xc = (-1,1))
    (D1, _, H, _) = diagonal_sbp_D1(p, N; xc = (-1,1))
    R = A - D1' * H * Diagonal(λ) * D1
    eigenvalues = eigen(Matrix(R)).values
    println("SBP order = $p")
    @show extrema(real.(eigenvalues))
    @show extrema(imag.(eigenvalues))
    println()
  end
end


text:
module TestTypes

using Test
using RLEVectors

@testset begin

    # Initialization
    expected_run_values = [1, 2, 3]
    expected_run_ends = [2, 4, 7]
    x = RLEVector([1, 1, 2, 2, 3, 3, 3])
    y = RLEVector(expected_run_values, expected_run_ends)
    @test isa(x, RLEVector)
    @test x.runvalues == expected_run_values
    @test x.runends == expected_run_ends
    x = RLEVector(["a", "b", "c", "d", "d", "d", "e", "f"], [3, 6, 9, 12, 15, 18, 21, 24])
    y = RLEVector(["a", "b", "c", "d", "e", "f"], [3, 6, 9, 18, 21, 24])
    @test x.runvalues == y.runvalues
    @test x.runends == y.runends
    x = RLEVector(BitVector(undef, 5), collect(1:5))
    y = RLEVector([false], [5])
    @test x.runvalues == y.runvalues
    @test x.runends == y.runends
    x = RLEVector([4, 5, 6], [2, 4, 6])
    y = RLEVector([4, 4, 5, 5, 6, 6])
    @test x.runvalues == y.runvalues
    @test x.runends == y.runends
    x = RLEVector(1:100 .< 50)
    @test x.runvalues == [true, false]
    @test x.runends == [49, 100]

    @test_throws ArgumentError RLEVector([1, 2, 3], [2, 4, 3]) # Sorted runends
    @test IntegerRle([1, 2, 3], [2, 4, 6]) == RLEVector([1, 2, 3], Int32[2, 4, 6])
    @test FloatRle([1, 2, 3], [2, 4, 6]) == RLEVector(Float64[1, 2, 3], Int32[2, 4, 6])
    @test BoolRle([true, false, true], [2, 4, 6]) == RLEVector([true, false, true], Int32[2, 4, 6])
    @test StringRle(["a", "bob", "joe"], [2, 4, 6]) ==
          RLEVector(["a", "bob", "joe"], Int32[2, 4, 6])

    @test RLEVector(5, 3) == RLEVector([5, 5, 5])

    sim = similar(RLEVector([1, 2, 3], [2, 4, 6]), Int32, (9))
    @test length(sim) == 9
    @test isa(values(sim), Vector{Int32}) == true

    # Conversion
    x = RLEVector([4, 4, 5, 5, 6, 7, 8])
    @test convert(Vector, x) == [4, 4, 5, 5, 6, 7, 8]
    @test convert(Set, x) == Set([4, 5, 6, 7, 8])
    @test convert(RLEVector, [1, 1, 2, 2]) == RLEVector([1, 2], [2, 4])
    x = RLEVector([4, 4, 5, 5, 6, 7, 8])
    y = convert(RLEVector{Float64,Int32}, x)
    @test typeof(y) == RLEVector{Float64,Int32}

    # Expanding
    @test collect(RLEVector([1, 1, 2, 2, 3, 3])) == [1, 1, 2, 2, 3, 3]
    @test convert(Vector, RLEVector([1, 1, 2, 2, 3, 3])) == [1, 1, 2, 2, 3, 3]

    ## Describing
    x = RLEVector(expected_run_values, expected_run_ends)
    @test eltype(x) == eltype(expected_run_values)
    io = IOBuffer()
    @test typeof(show(io, x)) == Nothing # At least test that show does not give error
    @test typeof(show(io, RLEVector(collect(1:100)))) == Nothing # At least test that show does not give error

    # Getters and setters
    x = RLEVector([1, 2, 3], [2, 9, 22])
    @test starts(x) == [1, 3, 10]
    @test starts(x, 1) == 1
    @test starts(x, 3) == 10
    @test widths(x) == [2, 7, 13]
    @test widths(x, 1) == 2
    @test widths(x, 2) == 7
    @test widths(x, 3) == 13
    @test ends(x) == [2, 9, 22]
    @test ends(x, 3) == 22
    @test values(x) == [1, 2, 3]
    @test widths!([2, 5, 9, 22]) == [2, 3, 4, 13]
    @test widths!([4]) == [4]
    @test widths!([]) == []

    # Hashing
    x = RLEVector([1, 1, 2, 2, 3, 3, 3])
    d = Dict("bob" => x)
    typeof(d["bob"]) == typeof(x)

end # testset

end # module


text:
module Inputs

using Compat
using PyFortran90Namelists: fstring
using QuantumESPRESSOBase.Inputs: QuantumESPRESSOInput, dropdefault, groupname

import AbInitioSoftwareBase.Inputs: FormatConfig, Namelist, asstring

FormatConfig(::Union{QuantumESPRESSOInput,Namelist}) = FormatConfig(;
    delimiter = " ",
    newline = "\n",
    indent = ' '^4,
    float = "%f",
    int = "%i",
    bool = ".%.",
)

"""
    asstring(input::QuantumESPRESSOInput)

Return a `String` representing a `QuantumESPRESSOInput`, valid for Quantum ESPRESSO's input.
"""
function asstring(input::QuantumESPRESSOInput)
    newline = FormatConfig(input).newline
    iter = Iterators.map(1:nfields(input)) do i
        x = getfield(input, i)
        x === nothing ? "" : asstring(x)
    end
    return join(iter, newline) * newline  # Add a new line at the end of line to prevent errors
end
"""
    asstring(nml::Namelist)

Return a `String` representing a `Namelist`, valid for Quantum ESPRESSO's input.
"""
function asstring(nml::Namelist)
    dict = dropdefault(nml)
    config = FormatConfig(nml)
    indent, delimiter, newline = config.indent, config.delimiter, config.newline
    iter = Iterators.map(dict) do (key, value)
        if value isa AbstractVector
            data = Iterators.map(enumerate(value)) do (i, x)
                if x !== nothing
                    indent * join((string(key, '(', i, ')'), "=", fstring(x)), delimiter)
                end
            end
            join(Iterators.filter(!isnothing, data), newline)
        elseif value isa NamedTuple
            data = Iterators.map(value) do (x, y)
                indent * join((string(key, '%', x), "=", fstring(y)), delimiter)
            end
            join(data, newline)
        else
            indent * join((string(key), "=", fstring(value)), delimiter)
        end
    end
    content = join(iter, newline)
    return join(filter(!isempty, ("&" * groupname(nml), content, '/')), newline)
end
asstring(str::AbstractString) = str

include("PWscf.jl")

end


text:
# revert JuliaLang/julia#32867; avoid string interpolation
#
# NOTE: without contextual dispatch, we can only redefine methods where a GPU-specific
#       type occurs in the signature (or we'll get a "fatal precompilation failure" error)

_bcs1(a::Integer, b::Integer) = a == 1 ? b : (b == 1 ? a : (a == b ? a : throw(DimensionMismatch("arrays could not be broadcast to a common size"))))
_bcs1(a::Integer, b) = a == 1 ? b : (first(b) == 1 && last(b) == a ? b : throw(DimensionMismatch("arrays could not be broadcast to a common size")))
_bcs1(a, b::Integer) = _bcs1(b, a)
_bcs1(a, b) = Broadcast._bcsm(b, a) ? Broadcast.axistype(b, a) : (Broadcast._bcsm(a, b) ? Broadcast.axistype(a, b) : throw(DimensionMismatch("arrays could not be broadcast to a common size")))

_bcs(::Tuple{}, ::Tuple{}) = ()
_bcs(::Tuple{}, newshape::Tuple) = (newshape[1], _bcs((), Base.tail(newshape))...)
_bcs(shape::Tuple, ::Tuple{}) = (shape[1], _bcs(Base.tail(shape), ())...)
function _bcs(shape::Tuple, newshape::Tuple)
    return (_bcs1(shape[1], newshape[1]), _bcs(Base.tail(shape), Base.tail(newshape))...)
end

broadcast_shape(shape::Tuple) = shape
broadcast_shape(shape::Tuple, shape1::Tuple, shapes::Tuple...) = broadcast_shape(_bcs(shape, shape1), shapes...)

@inline combine_axes(A, B...) = broadcast_shape(axes(A), combine_axes(B...))
combine_axes(A) = axes(A)

Broadcast._axes(::Broadcasted{<:AbstractGPUArrayStyle}, axes::Tuple) = axes
@inline Broadcast._axes(bc::Broadcasted{<:AbstractGPUArrayStyle}, ::Nothing) = combine_axes(bc.args...)


text:
vtrait(X) = X |> trait |> Val

const REQUIRE = "(requires MLJBase to be loaded)"

errlight(s) = throw(InterfaceError("Only `MLJModelInterface` is loaded. " *
                                    "Import `MLJBase` in order to use `$s`."))

# ------------------------------------------------------------------------
# categorical, note: not exported to avoid clashes; this is fine because
# MLJBase loads CategoricalArrays and MLJ interfaces should use qualified
# statements.

categorical(a...; kw...) = categorical(get_interface_mode(), a...; kw...)

categorical(::LightInterface, a...; kw...) = errlight("categorical")

# ------------------------------------------------------------------------
# matrix

"""
    matrix(X; transpose=false)

If `X <: AbstractMatrix`, return `X` or `permutedims(X)` if `transpose=true`.
If `X` is a Tables.jl compatible table source, convert `X` into a `Matrix`
$REQUIRE.
"""
matrix(X; kw...) = matrix(get_interface_mode(), vtrait(X), X; kw...)

matrix(::Mode, ::Val{:other}, X::AbstractMatrix; transpose=false) =
    transpose ? permutedims(X) : X

matrix(::Mode, ::Val{:other}, X; kw...) =
    throw(ArgumentError("Function `matrix` only supports AbstractMatrix or " *
                        "containers implementing the Tables interface."))

matrix(::LightInterface, ::Val{:table}, X; kw...) = errlight("matrix")

# ------------------------------------------------------------------------
# int

"""
   int(x; type=nothing)

The positional integer of the `CategoricalString` or `CategoricalValue` `x`, in
the ordering defined by the pool of `x`. The type of `int(x)` is the reference
type of `x` $REQUIRE.

Not to be confused with `x.ref`, which is unchanged by reordering of the pool
of `x`, but has the same type.

    int(X::CategoricalArray)
    int(W::Array{<:CategoricalString})
    int(W::Array{<:CategoricalValue})

Broadcasted versions of `int`.

    julia> v = categorical([:c, :b, :c, :a])
    julia> levels(v)
    3-element Array{Symbol,1}:
     :a
     :b
     :c
    julia> int(v)
    4-element Array{UInt32,1}:
     0x00000003
     0x00000002
     0x00000003
     0x00000001

See also: [`decoder`](@ref).
"""
function int(x; type::Union{Nothing,Type{T}}=nothing) where T <: Real
    type === nothing && return int(get_interface_mode(), x)
    return convert.(T, int(get_interface_mode(), x))
end

int(::LightInterface, x) = errlight("int")

# ------------------------------------------------------------------------
# classes

"""
    classes(x)

All the categorical elements with the same pool as `x` (including `x`),
returned as a list, with an ordering consistent with the pool $REQUIRE.
Here `x` has `CategoricalValue` or `CategoricalString` type, and `classes(x)`
is a vector of the same eltype. Note that `x in classes(x)` is always true.

Not to be confused with `levels(x.pool)`. See the example below.

    julia>  v = categorical([:c, :b, :c, :a])
    4-element CategoricalArrays.CategoricalArray{Symbol,1,UInt32}:
     :c
     :b
     :c
     :a

    julia> levels(v)
    3-element Array{Symbol,1}:
     :a
     :b
     :c

    julia> x = v[4]
    CategoricalArrays.CategoricalValue{Symbol,UInt32} :a

    julia> classes(x)
    3-element CategoricalArrays.CategoricalArray{Symbol,1,UInt32}:
     :a
     :b
     :c

    julia> levels(x.pool)
    3-element Array{Symbol,1}:
     :a
     :b
     :c

"""
classes(x) = classes(get_interface_mode(), x)

classes(::LightInterface, x) = errlight("classes")

# ------------------------------------------------------------------------
# decoder

"""
    d = decoder(x)

A callable object for decoding the integer representation of a
`CategoricalString` or `CategoricalValue` sharing the same pool as `x`
$REQUIRE. (Here `x` is of one of these two types.) Specifically, one has
`d(int(y)) == y` for all `y in classes(x)`. One can also call `d` on integer
arrays, in which case `d` is broadcast over all elements.

    julia> v = categorical([:c, :b, :c, :a])
    julia> int(v)
    4-element Array{UInt32,1}:
     0x00000003
     0x00000002
     0x00000003
     0x00000001
    julia> d = decoder(v[3])
    julia> d(int(v)) == v
    true

*Warning:* It is *not* true that `int(d(u)) == u` always holds.

See also: [`int`](@ref), [`classes`](@ref).
"""
decoder(x) = decoder(get_interface_mode(), x)

decoder(::LightInterface, x) = errlight("decoder")

# ------------------------------------------------------------------------
# table

"""
    table(columntable; prototype=nothing)

Convert a named tuple of vectors or tuples `columntable`, into a table of the
"preferred sink type" of `prototype` $REQUIRE. This is often the type of
`prototype` itself, when `prototype` is a sink; see the Tables.jl
documentation. If `prototype` is not specified, then a named tuple of vectors
is returned.

    table(A::AbstractMatrix; names=nothing, prototype=nothing)

Wrap an abstract matrix `A` as a Tables.jl compatible table with the specified
column `names` (a tuple of symbols). If `names` are not specified,
`names=(:x1, :x2, ..., :xn)` is used, where `n=size(A, 2)` $REQUIRE.

If a `prototype` is specified, then the matrix is materialized as a table of
the preferred sink type of `prototype`, rather than wrapped. Note that if
`prototype` is *not* specified, then `matrix(table(A))` is essentially a no-op.
"""
table(X; kw...) = table(get_interface_mode(), X; kw...)

table(::LightInterface, X; kw...) = errlight("table")

# ------------------------------------------------------------------------
# nrows, select, selectrows, selectcols

"""
    nrows(X)

Return the number of rows for a table, abstract vector or matrix `X` $REQUIRE.
"""
nrows(X) = nrows(get_interface_mode(), vtrait(X), X)

nrows(::Mode, ::Val{:other}, X::AbstractVecOrMat) = size(X, 1)
nrows(::Mode, ::Val{:other}, X::Nothing) = 0

nrows(::Mode, ::Val{:other}, X) =
    throw(ArgumentError("Function `nrows` only supports AbstractVector or " *
                        "AbstractMatrix or containers implementing the " *
                        "Tables interface."))

nrows(::LightInterface, ::Val{:table}, X) = errlight("table")

"""
    selectrows(X, r)

Select single or multiple rows from a table, abstract vector or matrix `X`
$REQUIRE. If `X` is tabular, the object returned is a table of the
preferred sink type of `typeof(X)`, even if only a single row is selected.
"""
selectrows(X, r) = selectrows(get_interface_mode(), vtrait(X), X, r)

selectrows(::Mode, ::Val{:other}, ::Nothing, r) = nothing

selectrows(::Mode, ::Val{:other}, X::AbstractVector, r)          = X[r]
selectrows(::Mode, ::Val{:other}, X::AbstractVector, r::Integer) = X[r:r]
selectrows(::Mode, ::Val{:other}, X::AbstractVector, ::Colon)    = X

selectrows(::Mode, ::Val{:other}, X::AbstractMatrix, r)          = X[r, :]
selectrows(::Mode, ::Val{:other}, X::AbstractMatrix, r::Integer) = X[r:r, :]
selectrows(::Mode, ::Val{:other}, X::AbstractMatrix, ::Colon)    = X

selectrows(::Mode, ::Val{:other}, X, r) =
    throw(ArgumentError("Function `selectrows` only supports AbstractVector " *
                        "or AbstractMatrix or containers implementing the " * "Tables interface."))

selectrows(::LightInterface, ::Val{:table}, X, r; kw...) =
    errlight("selectrows")

"""
    selectcols(X, c)

Select single or multiple columns from a matrix or table `X` $REQUIRE. If `c`
is an abstract vector of integers or symbols, then the object returned
is a table of the preferred sink type of `typeof(X)`. If `c` is a
*single* integer or column, then an `AbstractVector` is returned.
"""
selectcols(X, c) = selectcols(get_interface_mode(), vtrait(X), X, c)

selectcols(::Mode, ::Val{:other}, ::Nothing, c) = nothing

selectcols(::Mode, ::Val{:other}, X::AbstractMatrix, r)       = X[:, r]
selectcols(::Mode, ::Val{:other}, X::AbstractMatrix, ::Colon) = X

selectcols(::Mode, ::Val{:other}, X, r) =
    throw(ArgumentError("Function `selectcols` only supports AbstractMatrix " *
                        "or containers implementing the Tables interface."))

selectcols(::LightInterface, ::Val{:table}, X, c; kw...) =
    errlight("selectcols")

"""
    select(X, r, c)

Select element(s) of a table or matrix at row(s) r and column(s) c. An object
of the sink type of X (or a matrix) is returned unless c is a single integer or
symbol. In that case a vector is returned, unless r is a single integer, in
which case a single element is returned.

See also: [`selectrows`](@ref), [`selectcols`](@ref).
"""
select(X, r, c) = select(get_interface_mode(), vtrait(X), X, r, c)

# only used here to denote "group of indices"
const MIdx = Union{AbstractArray,Colon}

select(::Mode, ::Val, X, r::MIdx, c)       = selectcols(selectrows(X, r), c)
select(::Mode, ::Val, X, r, c::MIdx)       = selectcols(selectrows(X, r), c)
select(::Mode, ::Val, X, r::MIdx, c::MIdx) = selectcols(selectrows(X, r), c)
select(::Mode, ::Val, X, r, c) = _squeeze(selectcols(selectrows(X, r), c))

_squeeze(::Nothing) = nothing
_squeeze(v) = first(v)

# ------------------------------------------------------------------------
# UnivariateFinite

const UNIVARIATE_FINITE_DOCSTRING =
    """
        UnivariateFinite(classes, p)

A discrete univariate distribution whose finite support is the elements of the
vector `classes`, and whose corresponding probabilities are elements of the
vector `p`, which must sum to one $REQUIRE. Here `classes` must have type
`AbstractVector{<:CategoricalElement}` where

        CategoricalElement = Union{CategoricalValue,CategoricalString}

and all classes are assumed to share the same categorical pool.

        UnivariateFinite(prob_given_class)

A discrete univariate distribution whose finite support is the set of keys of
the provided dictionary, `prob_given_class` $REQUIRE. The dictionary keys must
be of type `CategoricalElement` (see above) and the dictionary values specify
the corresponding probabilities.
"""
UnivariateFinite(d::AbstractDict) = UnivariateFinite(get_interface_mode(), d)
UnivariateFinite(c::AbstractVector, p) =
    UnivariateFinite(get_interface_mode(), c, p)

UnivariateFinite(::LightInterface, a...) = errlight("UnivariateFinite")


text:
module churrosJulia

include("churro.jl")

export churro

end


text:
module CA

include("./Stochastic/Stochastic.jl")

end

module CleverAlgorithms
    using CA
end


text:
module TestCases

using Docile
@document

"f/0"
f() = ()

"f/1"
f(x) = ()

"f"
f

"g/0"
g() = ()

"g/1"
g(x) = ()

"g"
g

"K"
const K = ()

"@m"
macro m()
end

"S"
abstract S

"T"
type T
end

module A

using Docile
@document

"A.f/0"
f() = ()

"A.f/1"
f(x) = ()

"A.f"
f

"A.g/0"
g() = ()

"A.g/1"
g(x) = ()

"A.g"
g

"A.K"
const K = ()

"A.@m"
macro m()
end

"A.S"
abstract S

"A.T"
type T
end

module B

using Docile
@document

"A.B.f/0"
f() = ()

"A.B.f/1"
f(x) = ()

"A.B.f"
f

"A.B.g/0"
g() = ()

"A.B.g/1"
g(x) = ()

"A.B.g"
g

"A.B.K"
const K = ()

"A.B.@m"
macro m()
end

"A.B.S"
abstract S

"A.B.T"
type T
end

end

end

end


text:
export DSOSPoly, SDSOSPoly, SOSPoly

function JuMP.value(p::GramMatrix{JuMP.VariableRef})
    GramMatrix(map(JuMP.value, p.Q), p.x)
end

for poly in (:DSOSPoly, :SDSOSPoly, :SOSPoly)
    @eval begin
        struct $poly{PB<:PolyJuMP.AbstractPolynomialBasis} <: PolyJuMP.AbstractPoly
            polynomial_basis::PB
        end
        $poly(x::AbstractVector{<:MP.APL}) = $poly(MonomialBasis(x))
    end
end

matrix_cone_type(::SOSPoly) = MOI.PositiveSemidefiniteConeTriangle
matrix_cone_type(::DSOSPoly) = DiagonallyDominantConeTriangle
matrix_cone_type(::SDSOSPoly) = ScaledDiagonallyDominantConeTriangle

const PosPoly{PB} = Union{DSOSPoly{PB}, SDSOSPoly{PB}, SOSPoly{PB}}

JuMP.variable_type(m::JuMP.Model, p::PosPoly) = PolyJuMP.polytype(m, p, p.polynomial_basis)
PolyJuMP.polytype(m::JuMP.Model, ::PosPoly, basis::PolyJuMP.MonomialBasis{MT, MV}) where {MT<:MP.AbstractMonomial, MV<:AbstractVector{MT}} = GramMatrix{JuMP.VariableRef, MT, MV}

# Sum-of-Squares polynomial

_polytype(m::JuMP.Model, ::PosPoly, x::MVT) where {MT<:MP.AbstractMonomial, MVT<:AbstractVector{MT}} = GramMatrix{JuMP.VariableRef, MT, MVT}

function moi_add_variable(model::MOI.ModelLike, set::MOI.AbstractVectorSet,
                          binary::Bool, integer::Bool)
    VB = variable_bridge_type(typeof(set), Float64)
    Q, variable_bridge = add_variable_bridge(VB, model, set)
    if binary
        for q in Q
            MOI.add_constraint(model, q, MOI.ZeroOne())
        end
    end
    if integer
        for q in Q
            MOI.add_constraint(model, q, MOI.Integer())
        end
    end
    return Q
end

function JuMP.add_variable(model::JuMP.AbstractModel,
                           v::PolyJuMP.Variable{<:PosPoly{<:PolyJuMP.MonomialBasis}},
                           name::String="")
    monos = v.p.polynomial_basis.monomials
    set = matrix_cone(matrix_cone_type(v.p), length(monos))
    Q = moi_add_variable(backend(model), set, v.binary, v.integer)
    F = JuMP.jump_function_type(model, eltype(Q))
    return build_gram_matrix(
        F[JuMP.jump_function(model, vi) for vi in Q], monos)
end


text:
# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

"""
    ZerosBridge{T} <: Bridges.Variable.AbstractBridge

`ZerosBridge` implements the following reformulation:

* ``x \\in \\{0\\}`` into the substitution rule ``x = 0``,

where `T` is the coefficient type of `0`.

## Source node

`ZerosBridge` supports:

 * [`MOI.VectorOfVariables`](@ref) in [`MOI.Zeros`](@ref)

## Target nodes

`ZerosBridge` does not create target nodes. It replaces all instances of `x`
with `0` via substitution. This means that no variables are created in the
underlying model.

## Caveats

The bridged variables are similar to parameters with zero values. Parameters
with non-zero values can be created with constrained variables in
[`MOI.EqualTo`](@ref) by combining a [`VectorizeBridge`](@ref) and this bridge.

However, functions modified by `ZerosBridge` cannot be unbridged. That is, for a
given function, we cannot determine if the bridged variables were used.

A related implication is that this bridge does not support
[`MOI.ConstraintDual`](@ref). However, if a [`MOI.Utilities.CachingOptimizer`](@ref)
is used, the dual can be determined by the bridged optimizer using
[`MOI.Utilities.get_fallback`](@ref) because the caching optimizer records the
unbridged function.
"""
struct ZerosBridge{T} <: AbstractBridge
    n::Int # Number of variables
end

const Zeros{T,OT<:MOI.ModelLike} = SingleBridgeOptimizer{ZerosBridge{T},OT}

function bridge_constrained_variable(
    ::Type{ZerosBridge{T}},
    model::MOI.ModelLike,
    set::MOI.Zeros,
) where {T}
    return ZerosBridge{T}(MOI.dimension(set))
end

function supports_constrained_variable(::Type{<:ZerosBridge}, ::Type{MOI.Zeros})
    return true
end

function MOI.Bridges.added_constrained_variable_types(::Type{<:ZerosBridge})
    return Tuple{Type}[]
end

function MOI.Bridges.added_constraint_types(::Type{<:ZerosBridge})
    return Tuple{Type,Type}[]
end

# Attributes, Bridge acting as a model
MOI.get(::ZerosBridge, ::MOI.NumberOfVariables)::Int64 = 0

MOI.get(::ZerosBridge, ::MOI.ListOfVariableIndices) = MOI.VariableIndex[]

# References
MOI.delete(::MOI.ModelLike, ::ZerosBridge) = nothing

# Attributes, Bridge acting as a constraint

function MOI.get(::MOI.ModelLike, ::MOI.ConstraintSet, bridge::ZerosBridge)
    return MOI.Zeros(bridge.n)
end

function MOI.get(
    ::MOI.ModelLike,
    ::MOI.ConstraintPrimal,
    bridge::ZerosBridge{T},
) where {T}
    return zeros(T, bridge.n)
end

function MOI.get(
    ::MOI.ModelLike,
    ::MOI.ConstraintDual,
    ::ZerosBridge{T},
) where {T}
    return error(
        "Unable to query the dual of a variable bound that was reformulated " *
        "using `ZerosBridge`. This usually arises in conic solvers when a " *
        "variable is fixed to a value. As a work-around, instead of creating " *
        "a fixed variable using variable bounds like `p == 1`, add an affine " *
        "equality constraint like `1 * p == 1` (or `[1 * p - 1,] in Zeros(1)`).",
    )
end

function MOI.get(
    ::MOI.ModelLike,
    ::MOI.VariablePrimal,
    ::ZerosBridge{T},
    ::MOI.Bridges.IndexInVector,
) where {T}
    return zero(T)
end

function MOI.Bridges.bridged_function(
    ::ZerosBridge{T},
    ::MOI.Bridges.IndexInVector,
) where {T}
    return zero(MOI.ScalarAffineFunction{T})
end

function unbridged_map(
    ::ZerosBridge,
    ::MOI.VariableIndex,
    ::MOI.Bridges.IndexInVector,
)
    # The transformation is not recoverable, so we explicitly return `nothing`
    # here. See the docstring of `unbridged_map` for details.
    return nothing
end


text:
import Core: CodeInfo, ReturnNode, MethodInstance
import Core.Compiler: argextype, singleton_type
import Base.Meta: isexpr

argextype(@nospecialize args...) = argextype(args..., Any[])
code_typed1(args...; kwargs...) = first(only(code_typed(args...; kwargs...)))::CodeInfo
get_code(args...; kwargs...) = code_typed1(args...; kwargs...).code

# check if `x` is a statement with a given `head`
isnew(@nospecialize x) = isexpr(x, :new)
isreturn(@nospecialize x) = isa(x, ReturnNode)

# check if `x` is a dynamic call of a given function
iscall(y) = @nospecialize(x) -> iscall(y, x)
function iscall((src, f)::Tuple{CodeInfo,Base.Callable}, @nospecialize(x))
    return iscall(x) do @nospecialize x
        singleton_type(argextype(x, src)) === f
    end
end
iscall(pred::Base.Callable, @nospecialize(x)) = isexpr(x, :call) && pred(x.args[1])

# check if `x` is a statically-resolved call of a function whose name is `sym`
isinvoke(y) = @nospecialize(x) -> isinvoke(y, x)
isinvoke(sym::Symbol, @nospecialize(x)) = isinvoke(mi->mi.def.name===sym, x)
isinvoke(pred::Function, @nospecialize(x)) = isexpr(x, :invoke) && pred(x.args[1]::MethodInstance)

function fully_eliminated(@nospecialize args...; retval=(@__FILE__), kwargs...)
    code = code_typed1(args...; kwargs...).code
    if retval !== (@__FILE__)
        return length(code) == 1 && isreturn(code[1]) && code[1].val == retval
    else
        return length(code) == 1 && isreturn(code[1])
    end
end

max_stars_repo_name: JuliaLabs/julia

max_issues_repo_name: JuliaLabs/julia

max_forks_repo_name: JuliaLabs/julia


text:
using Test
using IntervalArithmetic

@testset "atan2" begin include("atan2.jl") end
@testset "c-xsc" begin include("c-xsc.jl") end
@testset "fi_lib" begin include("fi_lib.jl") end
@testset "IEEE 1788" begin
    @testset "ieee1788-constructors" begin include("ieee1788-constructors.jl") end
    @testset "ieee1788-exceptions" begin include("ieee1788-exceptions.jl") end
    @testset "libieeep1788_bool" begin include("libieeep1788_bool.jl") end
    @testset "libieeep1788_cancel" begin include("libieeep1788_cancel.jl") end
    @testset "libieeep1788_class" begin include("libieeep1788_class.jl") end
    @testset "libieeep1788_elem" begin include("libieeep1788_elem.jl") end
    @testset "libieeep1788_num" begin include("libieeep1788_num.jl") end
    @testset "libieeep1788_overlap" begin include("libieeep1788_overlap.jl") end
    @testset "libieeep1788_rec_bool" begin include("libieeep1788_rec_bool.jl") end
    @testset "libieeep1788_reduction" begin include("libieeep1788_reduction.jl") end
    @testset "libieeep1788_set" begin include("libieeep1788_set.jl") end
    @testset "libieeep1788_rev" begin include("libieeep1788_rev.jl") end
    @testset "libieeep1788_mul_rev" begin include("libieeep1788_mul_rev.jl") end
end # IEEE 1788
@testset "mpfi" begin include("mpfi.jl") end
@testset "abs_rev" begin include("abs_rev.jl") end
@testset "pow_rev" begin include("pow_rev.jl") end


text:
include("test.jl")
include("test_psl.jl")


text:
# set random seed to promote repeatability in CI unit tests
using Random
Random.seed!(101)

for filename in (
        "ginsu.jl",
        "jop_prop2DAcoIsoDenQ_DEO2_FDTD.jl",
        "jop_prop2DAcoTTIDenQ_DEO2_FDTD.jl",
        "jop_prop2DAcoVTIDenQ_DEO2_FDTD.jl",
        "jop_prop3DAcoIsoDenQ_DEO2_FDTD.jl",
        "jop_prop3DAcoTTIDenQ_DEO2_FDTD.jl",
        "jop_prop3DAcoVTIDenQ_DEO2_FDTD.jl",
        "jop_sinc_regular.jl",
        "wavefield_separation.jl"
    )
    include(filename)
end


text:
# ------------------------------------------------------------------
# Copyright (c) 2017, Júlio Hoffimann Mendes <juliohm@stanford.edu>
# Licensed under the ISC License. See LICENCE in the project root.
# ------------------------------------------------------------------

"""
    PointSet(coords)

A set of points with coordinate matrix `coords`. The number of rows
of the matrix is the dimensionality of the domain whereas the number
of columns is the number of points in the set.
"""
struct PointSet{T<:Real,N} <: AbstractDomain{T,N}
  coords::Matrix{T}

  function PointSet{T,N}(coords) where {N,T<:Real}
    @assert !isempty(coords) "coordinates must be non-empty"
    new(coords)
  end
end

PointSet(coords::AbstractMatrix{T}) where {T<:Real} =
  PointSet{T,size(coords,1)}(coords)

npoints(ps::PointSet) = size(ps.coords, 2)

function coordinates!(buff::AbstractVector{T}, ps::PointSet{T,N},
                      location::Int) where {N,T<:Real}
  for i in 1:N
    @inbounds buff[i] = ps.coords[i,location]
  end
end

# ------------
# IO methods
# ------------
function Base.show(io::IO, ps::PointSet{T,N}) where {N,T<:Real}
  dims = size(ps.coords, 1)
  npts = size(ps.coords, 2)
  print(io, "$dims×$npts PointSet{$T,$N}")
end

function Base.show(io::IO, ::MIME"text/plain", ps::PointSet{T,N}) where {N,T<:Real}
  println(io, ps)
  Base.print_array(io, ps.coords)
end


text:
# ---
# title: 1079. Letter Tile Possibilities
# id: problem1079
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Backtracking
# link: <https://leetcode.com/problems/letter-tile-possibilities/description/>
# hidden: true
# ---
#
# You have `n`  `tiles`, where each tile has one letter `tiles[i]` printed on
# it.
#
# Return _the number of possible non-empty sequences of letters_ you can make
# using the letters printed on those `tiles`.
#
#
#
# **Example 1:**
#
#
#
#     Input: tiles = "AAB"
#     Output: 8
#     Explanation: The possible sequences are "A", "B", "AA", "AB", "BA", "AAB", "ABA", "BAA".
#
#
# **Example 2:**
#
#
#
#     Input: tiles = "AAABBC"
#     Output: 188
#
#
# **Example 3:**
#
#
#
#     Input: tiles = "V"
#     Output: 1
#
#
#
#
# **Constraints:**
#
#   * `1 <= tiles.length <= 7`
#   * `tiles` consists of uppercase English letters.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end