text:
module HFMod

if length(findin("C:\\Users\\Clinton\\Dropbox\\Projects\\SummerProject17",LOAD_PATH)) == 0
  push!(LOAD_PATH,"C:\\Users\\Clinton\\Dropbox\\Projects\\SummerProject17")
end

#=TODO:
1) IV=#
#2) Run specifications

using DataFrames, Distributions, StatsBase, GZip, JLD, Gadfly, CTMod
#importall CT

#Set this to an optimal value- I find logical cores *(3/4) works well
BLAS.set_num_threads(Int(round((Sys.CPU_CORES*3÷4))))

#if !isdefined(:constDef)
  const constDef = true
  #other constants here

  const DATA_PATH = pwd() * "\\data"
  const DATA_FILE_NAME = "HFPerformance_2015-2017"
  const DATA_FILE_NAME_DEAD = DATA_FILE_NAME * "_dead"
  const DATE_FORMAT_STR = "yyyymmdd"
  const NOTICE_PERIOD_ADJ = 7 #days to add to the notice (since data are monthly)
  const MIN_PERIODS_FOR_INCLUSION = 6 #number of periods required for inclusion
  const T_TEST_THRESHOLD = 2.0 #this is a parameter for identificaiton if the flows are lagged
  const PERIODS_SAME_AUM_TO_DELETE = 4 #assume 4 consecutive periods of identical aum implies stale data
  const LAGS_NOTICE_TO_CAPTURE = 3
  const LAGS_PERFORMANCE_TO_CAPTURE = 12
  const MIN_ASSETS = 10.0 #10,000,000
  const START_LAGS_FROM_NOTICE = true #start the performance lags from the notice period (not the redemption date)
  const RUN_LM_INTERACTIONS = false #run the interaction regressions (takes a while)
  const RUN_2SLS_INTERACTIONS =  false
  const RESULTS_PATH = pwd() * "\\results"
  const FOOTER_NAME = "footer.tex"
  const HEADER_NAME = "header.tex"
  const USE_AGGRESSIVE_GC = false #this enables very agressive GC which helps memory management at cost of performance
  const R_TEST_PATH = "C:\\Users\\Clinton\\Dropbox\\Projects\\Summer17RTester\\Summer17RTester"

#end

###############HELPER FUNCTIONS

#helper function to extract all positive numbers from a string
#IN: A string
#OUT: A vector of floats, possibly empty
function getPosNumsFromString(s::String)
  len::Int = length(s)
  preParseS::Vector{Char} = collect(s)

  #replace non-numerical values with delimitor Ξ
  for i ∈ 1:len
     preParseS[i] = (isnumber(s[i]) || (
      s[i]=='.' && (i>1 && isnumber(s[i-1]) || (
        i<len && isnumber(s[i+1])))))?s[i]:'Ξ'
  end

  #split the tokens, keeping only numberical strings
  sSplit::Vector{String} = split(preParseS,'Ξ',keep=false)

  #parse and return the values as an array of floats
  return broadcast((x::String)->Float64(parse(x)),sSplit)::Vector{Float64}
end

#=helper function which sorts and splits DF on symbol, extracts all positive
numeric values, averages them, and replaces the column
IN: DataFrame, column symbol, a default value, and whetehr to leave NAs
OUT: None explicit. DataFrame modified with numberic column named col,
old column is named col_str=#
function DFStringtoMeanNum!(DF::DataFrame, col::Symbol;
  default::Float64=0.0, leaveNAs::Bool=false, silent::Bool=true)

  oldCol = Symbol(col, "_str")
  rename!(DF, col, oldCol)
  DF[:,col] = default
  sort!(DF, cols=(oldCol))

  by(DF, oldCol) do DFSub::SubDataFrame
    if ~isna(DFSub[1,oldCol])
      parsed::Vector{Float64} = #parse the string
        getPosNumsFromString(DFSub[1,oldCol])
      if length(parsed)>0 #if string wasn't all garbage
        DFSub[:,col] = mean(parsed)
      elseif !silent#Give a heads up given a garbage incentive fee
        println("WARNING: Could not parse ", col, ". Value: ", DFSub[1, oldCol],
          ", Assuming numeric value of ", default)
      end
    elseif leaveNAs
      DFSub[:,col] = NA
    end
  end
end

#takes a vector and writes a string
function vec2String(v::Vector{T}, delim::W = "")::String where
  {T<:Any, W<:Union{AbstractString,Char}}

  return join(String.(v), delim)
end

#calculates the t statistic of the Pearson Correlation (Lowry 2017)
function tStatOfCor(x::Vector{Float64},y::Vector{Float64})
  r::Float64 = cor(x,y)
  return r/sqrt((1-r^2.0)/(length(x)-2))::Float64
end

#this is a type which stores all of the notice lags
struct NoticeLagSymbols
  noticeLags::Vector{Symbol}
  noticeLagsNo1st::Vector{Symbol}
  noticeLLags::Vector{Symbol}
  noticeLLagsNo1st::Vector{Symbol}
end

#Constructor for the above type
#IN: Number of lags
#OUT: A container of symbols with the appropraite lags
function NoticeLagSymbols(lags::Int)::NoticeLagSymbols
  noticeLags::Vector{Symbol} = Vector{Symbol}(lags)
  noticeLagsNo1st::Vector{Symbol} = similar(noticeLags)
  noticeLLags::Vector{Symbol} = similar(noticeLags)
  noticeLLagsNo1st::Vector{Symbol} = similar(noticeLags)

  for i ∈ 1:lags
    noticeLags[i] = Symbol("noticeLag$i")
    noticeLagsNo1st[i] = Symbol("noticeLag$(i)No1st")
    noticeLLags[i] = Symbol("noticeLLag$i")
    noticeLLagsNo1st[i] = Symbol("noticeLLag$(i)No1st")
  end

  return NoticeLagSymbols(noticeLags,noticeLagsNo1st, noticeLLags,
    noticeLLagsNo1st)
end

#this is a type which stores all of the performance lags
struct PerformanceLagSymbols
  performanceLags::Vector{Symbol}
end

#Constructor for the above type
#IN: Number of lags
#OUT: A container of symbols with the appropraite lags
function PerformanceLagSymbols(lags::Int)::PerformanceLagSymbols

  performanceLags::Vector{Symbol} = Vector{Symbol}(lags)   #preallocate

  for i ∈ 1:lags
    performanceLags[i] = Symbol("performanceLag$i")
  end

  return PerformanceLagSymbols(performanceLags)
end

contains(e::T where T<:CTExpr, s::String)::Bool  = Base.contains("$e",s)
###########Main Methods
#U for unlagged, L for lagged, C for corrected, R for corrected but with
#a lagged default

function PreProcessHFData(;FlowsLaggingProcedure::Char = 'U', newBinary::Bool = false)

  nPer::Int = 0 #numer of periods
  dtFormat::DateFormat = DateFormat(DATE_FORMAT_STR) ## a date format object

  #this allows us to capture a variable amount of lags
  noticeLS::NoticeLagSymbols = NoticeLagSymbols(LAGS_NOTICE_TO_CAPTURE)
  performanceLS::PerformanceLagSymbols = PerformanceLagSymbols(LAGS_PERFORMANCE_TO_CAPTURE)

  #this makes sure we have a binary of the data (Improves load times 3-4x)
  if !isfile("$DATA_PATH\\$DATA_FILE_NAME.jls") || newBinary
    #extract from the zip file
    gHandle::GZipStream = GZip.open("$DATA_PATH\\$DATA_FILE_NAME.gz")
    write("$DATA_PATH\\$DATA_FILE_NAME.csv", readlines(gHandle,chomp=false))
    close(gHandle)

    gHandle = GZip.open("$DATA_PATH\\$DATA_FILE_NAME_DEAD.gz")
    write("$DATA_PATH\\$DATA_FILE_NAME_DEAD.csv", readlines(gHandle,chomp=false))
    close(gHandle)

    #write the binary
    oStream::IOStream = open("$DATA_PATH\\$DATA_FILE_NAME.jls","w+")
    serialize(oStream,
        [readtable("$DATA_PATH\\$DATA_FILE_NAME.csv"); readtable("$DATA_PATH\\$DATA_FILE_NAME_DEAD.csv")])
    close(oStream)
  end

  #load the binary
  iStream::IOStream = open("$DATA_PATH\\$DATA_FILE_NAME.jls")
  HFData::DataFrame = deserialize(iStream)
  close(iStream)

  println("Initial rows: ", size(HFData,1), ", Initial columns: ", size(HFData,2))
  println("Begining pre-processing...")

  #drop some columns
  HFData = HFData[:, [:fund_id, :inception, :main_strategy,
    :sub_strategy, :leverage, :fund_assets, :firm_assets,
    :returns_denomination, :management_fee, :incentive_fee, :high_watermark,
    :hurdle_rate, :subscriptions, :redemptions,
    :advance_notice, :lockup, :fund_assets_as_of,
    :fund_assets_denomin, :sales_fee, :other_fees,
    :date_added_to_db, :fund_status, :ucitsiii, :date, :performance,
    :nav, :assets]]

  #Drop some null valued fields:
  HFData = dropNullsFromDF(HFData, [:date, :performance, :assets,
  :returns_denomination, :management_fee,:advance_notice])

  #only want active funds
  #HFData = HFData[HFData[:,:fund_status].=="Active",:]

  delete!(HFData,:fund_status) #no need for this anymore
  HFData = HFData[HFData[:,:assets] .≥ MIN_ASSETS, :]

  #parse the incentive fees and mangement fees
  DFStringtoMeanNum!(HFData, :incentive_fee)
  DFStringtoMeanNum!(HFData, :management_fee)

  #convert from percentage to decimal
  HFData[:,:incentive_fee] ./= 100.0
  HFData[:,:management_fee] ./= 100.0
  HFData[:,:performance] ./= 100.0

  #broadcast an anonymous functino to do the type conversions
  HFData[:,:inception] = ((x::Int)->Date("$x",dtFormat)).(HFData[:,:inception])
  HFData[:,:date_added_to_db] = ((x::Int)->Date("$x",dtFormat)).(HFData[:,:date_added_to_db])
  HFData[:,:date] = ((x::Int)->Date("$x",dtFormat)).(HFData[:,:date])
  HFData[:,:fund_assets_as_of] = ((x::Int)->Date("$x",dtFormat)).(HFData[:,:fund_assets_as_of])

  #sort the data by fund id and date
  sort!(HFData,cols = (:fund_id,:date))
  #delStaleAUMDat!(HFData)

  #clean out some funds with stale aum data
  by(HFData, :fund_id) do HFSub::SubDataFrame
    nPer = size(HFSub,1)

    if nPer > PERIODS_SAME_AUM_TO_DELETE #check if this procedure makes sense
      identicalAUMs::Int = 1
      i::Int =  2
      while i ≤ nPer && identicalAUMs < PERIODS_SAME_AUM_TO_DELETE
        identicalAUMs = (HFSub[i,:assets]==HFSub[i-1,:assets])?identicalAUMs+1:1
        i+=1
      end
      if identicalAUMs == PERIODS_SAME_AUM_TO_DELETE #set assets to NA to flag for deletion
        HFSub[:,:assets] = NA
      end
    end
  end

  HFData =  HFData[~isna(HFData[:,:assets]),:] #delete the identified funds

  #delete funds with gaps in their performance history, or funds with a
  #track record less than six months
  HFData[:,:mcnt]=0
  by(HFData, :fund_id) do HFSub::SubDataFrame
    mcnt::Int = size(HFSub,1)
    #now check for track record breaks
    if maximum(HFSub[:,:date]) ≠
      Dates.lastdayofmonth(minimum(HFSub[:,:date])+Dates.Month(mcnt)-Dates.Month(1))
      mcnt = 0
    end
    HFSub[:,:mcnt] = mcnt
  end
  HFData=HFData[~(HFData[:,:mcnt].<MIN_PERIODS_FOR_INCLUSION),:]

  #get the months of advance notice
  HFData[:,:months_notice] = broadcast((x::Int)->
    ((x+NOTICE_PERIOD_ADJ)÷30)::Int,HFData[:,:advance_notice])

  ###The following sections create nessecary columns with placeholders
  #calculate monthly flows and construct the data set
  HFData[:,:flows] = 0.0 #flows
  HFData[:,:redemp_notice] = 0.0 #notifications about flows
  HFData[:,:rolling_notice] = 0.0 #notifications about flows
  HFData[:,:lFlows]  = 0.0 #will hold log flows
  HFData[:,:redemp_notice_dt] = Date(1,1,1)

  #create the appropriate columns for lagged flow notice
  for i::Int ∈ 1:LAGS_NOTICE_TO_CAPTURE
    HFData[:,noticeLS.noticeLags[i]] = 0.0
    HFData[:,noticeLS.noticeLagsNo1st[i]] = 0.0
    HFData[:,noticeLS.noticeLLags[i]] = 0.0
    HFData[:,noticeLS.noticeLLagsNo1st[i]] = 0.0
  end

  #create columns for performance lag
  for i::Int ∈ 1:LAGS_PERFORMANCE_TO_CAPTURE
    HFData[:,performanceLS.performanceLags[i]] = 0.0
    HFData[:,performanceLS.performanceLags[i]] .= NA
  end

  ctr = 0
  ###This is the main section for manipulating data by fund. We will
  #adjust the lagging of assets and assign to appropriate lagged variables
  by(HFData, :fund_id) do HFSub::SubDataFrame
    nPer = HFSub[1,:mcnt]
    useLaggingProc::Bool = false

    #now calculate the flows in the lagged and unlagged flows
    flowsUnlagged::Vector{Float64} =Vector{Float64}(nPer-1)
    flowsLagged::Vector{Float64} =Vector{Float64}(nPer-1)
    for i ∈ 2:nPer
      #back out the flows:
      #F_unlagged=A(t)-A(t-1)*(1+r(t))
      flowsUnlagged[i-1] = HFSub[i,:assets] -
        (HFSub[i-1,:assets])*(1.0+HFSub[i,:performance])

      #F_lagged=A(t)-A(t-1)*(1+r(t-1))
      flowsLagged[i-1]  = HFSub[i,:assets] -
        (HFSub[i-1,:assets])*(1.0+HFSub[i-1,:performance])
    end

    if FlowsLaggingProcedure == 'U'  #assume all unlagged
      useLaggingProc = false

    elseif FlowsLaggingProcedure == 'L' #assume all lagged
      useLaggingProc = true

      #below assumes a comparison process
    elseif FlowsLaggingProcedure == 'C' || FlowsLaggingProcedure == 'R'

      #get the differencein performance between lagged and unlagged
      perfDelta::Vector{Float64} = HFSub[2:end,:performance]
        .-HFSub[1:(end-1),:performance]

      #check the correlations
      flowsUnlaggedT::Float64 = tStatOfCor(flowsUnlagged,perfDelta)
      flowsLaggedT::Float64 = tStatOfCor(flowsLagged,perfDelta)

      #='C' Logic: Assume unlagged. If unlagged T test for correlation against
      return deltas is <-2 (or some other threshold), switch to lagged unless
      the lagged T stat is greater than 2
      'R' Logic: Same as C but reverse.=#
      if flowsUnlaggedT > -1.0 * T_TEST_THRESHOLD && flowsLaggedT > 1.0 * T_TEST_THRESHOLD
        useLaggingProc = false
      elseif flowsUnlaggedT < -1.0 * T_TEST_THRESHOLD && flowsLaggedT < 1.0 * T_TEST_THRESHOLD
        useLaggingProc = true
      elseif flowsUnlaggedT * flowsLaggedT ≤ 0.0
        if FlowsLaggingProcedure == 'C'
          useLaggingProc = false
        elseif FlowsLaggingProcedure == 'R'
          useLaggingProc = true
        end
      end
    end

    #now execute the appropriate lagging procedure
    if useLaggingProc
      HFSub[end,:flows] = NA
      HFSub[1:(end-1),:flows] = flowsLagged
      HFSub[1:(end-1),:lFlows] = flowsLagged ./
        (HFSub[2:end,:assets] .-  flowsLagged)
    else
      HFSub[1,:flows] = NA
      HFSub[1,:lFlows] = NA
      HFSub[2:end,:flows] = flowsUnlagged
      HFSub[2:end,:lFlows] = flowsUnlagged ./
        (HFSub[2:end,:assets] .- flowsUnlagged)
    end

    #now we look up the notice period to determine the appropriate month
    #Only will be for negative flows
    #we will make the table two way, at the expense of redundancy
    for i::Int ∈ 1:nPer
      if !isna(HFSub[i,:flows]) && HFSub[i,:flows] < 0.0 &&
        i::Int > HFSub[i,:months_notice]
        HFSub[i-HFSub[i,:months_notice],:redemp_notice] = HFSub[i,:flows]
        HFSub[(i-HFSub[i,:months_notice]):(i),:rolling_notice] += HFSub[i,:flows]
        HFSub[i-HFSub[i,:months_notice],:redemp_notice_dt] = HFSub[i,:date]

        #record the redemption notificaion lags and performance lags
        if HFSub[i,:months_notice] > 0
          for l::Int ∈ 1:min(LAGS_NOTICE_TO_CAPTURE, HFSub[i,:months_notice])
            HFSub[i-l, noticeLS.noticeLags[l]] = HFSub[i,:flows]
            HFSub[i-l, noticeLS.noticeLLags[l]] = HFSub[i,:lFlows]
          end
        end

        #record the redemption notificaion lags, excluding the first months' notice
        if HFSub[i,:months_notice] > 1
          for l::Int ∈ 1:min(LAGS_NOTICE_TO_CAPTURE, HFSub[i,:months_notice]-1)
            HFSub[i-l, noticeLS.noticeLagsNo1st[l]] = HFSub[i,:flows]
            HFSub[i-l, noticeLS.noticeLLagsNo1st[l]] = HFSub[i,:lFlows]
          end
        end

      end
    end

    #now assign to the performance lags
    if START_LAGS_FROM_NOTICE #if true, we start from the end of the notice period
      for i::Int ∈ 2:nPer
        if i - 1 - HFSub[i,:months_notice] ≥ 1 #(make sure its positive number)
          for l::Int ∈ 1:min(i-1- HFSub[i,:months_notice],LAGS_PERFORMANCE_TO_CAPTURE)
            HFSub[i, performanceLS.performanceLags[l]] = HFSub[i - l, :performance]
          end
        end
      end
    else
      for i::Int ∈ 2:nPer
        for l::Int ∈ 1:min(i-1,LAGS_PERFORMANCE_TO_CAPTURE)
          HFSub[i, performanceLS.performanceLags[l]] = HFSub[i - l, :performance]
        end
      end
    end
  end

  #get a seperate column with net negative flows
  HFData[:,:negLFlows] = HFData[:,:lFlows]
  HFData[(~isna(HFData[:,:lFlows])) & (HFData[:,:lFlows].≥0.0),:negLFlows] .= NA

  #convert string columns to factors
  HFTypes::Vector{Type} = eltypes(HFData) #all types
  HFNames::Vector{String} = names(HFData) #all names
  for i ∈ 1:size(HFData,2) #for all columns
    if HFTypes[i] == String && !(typeof(HFData[:,i]) <: PooledDataVector)
      pool!(HFData,Symbol(HFNames[i])) #convert to factor
    end
  end

  #create a factor column for dates
  HFData[:,:fDate] = HFData[:,:date]
  pool!(HFData,:fDate)

  println("Processing complete.")
  println("Ending rows: ", size(HFData,1), ", Ending columns: ", size(HFData,2))

  #save the binary
  oStream = open("$DATA_PATH\\$(DATA_FILE_NAME)_$(FlowsLaggingProcedure)_clean.jls","w+")
  serialize(oStream, HFData)
  close(oStream)

end

#returns a table with the total flows
#High performance requires an initial sort
function aggregateOnDate(HFData::DataFrame, lags::Int=5,
  noFirstNoticeMonth::Bool=false)

  nDates::Int = size(unique(HFData[:,:date]),1)
  noticeLS::NoticeLagSymbols = NoticeLagSymbols(lags)

  #preallocate data frame
  HFByDate::DataFrame =
    DataFrame([Date, Int, Float64, Float64, Float64, Float64,Float64],
    [:date, :numFunds, :totAUM, :totFlows, :totRedemp,
    :totRollingRedemp, :totRedempNotice], nDates)
  HFByDate[:,:date] .= unique(HFData[:,:date])

  #allocate the lag columns
  for s ∈ [noticeLS.noticeLags; noticeLS.noticeLagsNo1st]

    HFByDate[:,s] = 0.0
  end

  #make a lookup table
  dateTbl = Dict(HFByDate[i,:date] => i for i::Int ∈ 1:nDates)

  #might be a better way to do this with aggregates and joins,
  #although this is reasonably fast with a low number of comparisons
  #Basically we are getting totals for each symbol by a bunch of symbols by date
  by(HFData, :date) do HFSub::SubDataFrame
    dateIdx::Int = dateTbl[HFSub[1,:date]]
    applyToDate::Vector{Bool} = HFByDate[:,:date].==HFSub[1,:date]
    HFByDate[dateIdx,:numFunds]=size(unique(HFSub[:,:fund_id]),1)
    HFByDate[dateIdx,:totAUM]=sum(HFSub[~isna(HFSub[:,:assets]),:assets])
    HFByDate[dateIdx,:totFlows]=sum(HFSub[~isna(HFSub[:,:flows]),:flows])
    HFByDate[dateIdx,:totRedemp]=
      sum(HFSub[broadcast((x)->ifelse(~isna(x)&&x<0,true, false),HFSub[:,:flows]),:flows])
    HFByDate[dateIdx,:totRedempNotice]=
      sum(HFSub[~isna(HFSub[:,:redemp_notice]),:redemp_notice])
    HFByDate[dateIdx,:totRollingRedemp]=
      sum(HFSub[~isna(HFSub[:,:rolling_notice]),:rolling_notice])
    for s ∈ [noticeLS.noticeLags; noticeLS.noticeLagsNo1st]
      HFByDate[dateIdx,s]= sum(HFSub[~isna(HFSub[:,s]),s])
    end

  end

  return HFByDate

end

#this is the main method for analysis
function AnalyzeAggregates(;FlowsLaggingProcedure::Char = 'U',
    verbose::Bool = true)

    #open the pre-processed data file
    iStream::IOStream =
        open("$DATA_PATH\\$(DATA_FILE_NAME)_$(FlowsLaggingProcedure)_clean.jls")
    HFData::DataFrame = deserialize(iStream)

    sort!(HFData,cols = (:date, :fund_id))
    HFByDate::DataFrame = aggregateOnDate(HFData, LAGS_NOTICE_TO_CAPTURE)
    performanceLS::PerformanceLagSymbols =
        PerformanceLagSymbols(LAGS_PERFORMANCE_TO_CAPTURE)

    XLMSpecs::Vector{CTExpr} = [parse("negLFlows")]
    XLMNames::Vector{Vector{Symbol}} = [[:intercept; :negLFlows]] #note only non-factor names
    #####put additional specifications here
    push!(XLMSpecs, parse("negLFlows + fDate"))
    push!(XLMNames, [:intercept; :negLFlows])

    push!(XLMSpecs, parse("negLFlows + main_strategy"))
    push!(XLMNames, [:intercept; :negLFlows])

    push!(XLMSpecs, parse("negLFlows + " *
        "$(vec2String(performanceLS.performanceLags, "+"))"))
    push!(XLMNames, [:intercept; :negLFlows; performanceLS.performanceLags])

    push!(XLMSpecs, parse("negLFlows + " *
    "$(vec2String(performanceLS.performanceLags, "+")) + fDate"))
    push!(XLMNames, [:intercept; :negLFlows; performanceLS.performanceLags])

    push!(XLMSpecs, parse("negLFlows + " *
    "$(vec2String(performanceLS.performanceLags, "+")) + main_strategy"))
    push!(XLMNames, [:intercept; :negLFlows; performanceLS.performanceLags])

    push!(XLMSpecs, parse("negLFlows + " *
    "$(vec2String(performanceLS.performanceLags, "+")) + fDate + main_strategy"))
    push!(XLMNames, [:intercept; :negLFlows; performanceLS.performanceLags])

    if RUN_LM_INTERACTIONS
        push!(XLMSpecs, parse("negLFlows + fDate * main_strategy"))
        push!(XLMNames, [:intercept; :negLFlows])
        push!(XLMSpecs, parse("negLFlows + "  *
          "$(vec2String(performanceLS.performanceLags, "+")) + " *
          "fDate * main_strategy"))
        push!(XLMNames, [:intercept; :negLFlows; performanceLS.performanceLags])
    end
    #####

    modelsLM::Vector{CTLM} = Vector{CTLM}()
    #run the LM Specs
    for i ∈ 1:length(XLMSpecs)

        if verbose
            println("Begin Spec RHS: $(XLMSpecs[i])")
        end

        push!(modelsLM, CTLM(HFData, XLMSpecs[i],
            :performance, XNames=XLMNames[i], YName = :performance))

        if verbose
            println("Linear Model X Specification: $(XLMSpecs[i])")
            println("Names: $(XLMNames[i])")
            println("β: $(modelsLM[i].β[1:min(10,end)])")
            println("σ: $(sqrt.(diag(getHomoskedΣ!(modelsLM[i])))[1:min(10,end)])")
            println("Data points: $(size(modelsLM[i].X,1))\n")
        end

    end
    ##########LM IO
    descRowNamesLM::Vector{String} = ["Performance Lags", "Date Fixed Effects",
        "Strategy Fixed Effects", "F.E. Interactions", "N ('000s)"]

    #need to print the descriptive rows
    descContentLM::Vector{Vector{String}} =
        #[fill("",length(modelsLM)) for i∈ 1:length(descRowNamesLM)]
        [Vector{String}(length(modelsLM))for i∈ 1:length(descRowNamesLM)]

    #this builds the descriptive rows. There is Probably a better way to do this,
    #but its fairly specific to the project.
    for i ∈ 1:length(XLMSpecs)
        descContentLM[1][i] = "$(contains(XLMSpecs[i], "performanceLag")?"X":"")"
        descContentLM[2][i] = "$(contains(XLMSpecs[i], "fDate")?"X":"")"
        descContentLM[3][i] = "$(contains(XLMSpecs[i], "main_strategy")?"X":"")"
        descContentLM[4][i] = "$(contains(XLMSpecs[i], "*")?"X":"")"
        descContentLM[5][i] = "$(round(modelsLM[i].N/1000.0,1))"
    end

    TableTextLM::String = texTable(modelsLM, getHomoskedΣ!, [:intercept; :negLFlows],
        caption = "Performance Flows OLS Specifications",
        colNames = [["($i)" for i∈1:length(modelsLM)]],
        contentRowNames = ["intercept", "outflows/AUM x 100"],
        descRowNames = descRowNamesLM,
        descContent = descContentLM,
        decimalDigits = 3,
        columnSepPt = -20,
        scaling = [1000.0,1000.0],
        clearMem = USE_AGGRESSIVE_GC)

    if USE_AGGRESSIVE_GC
        gc()
    end

    #####################2SLS code###################
    X2SLSSpecs::Vector{CTExpr} = [parse("negLFlows+0")]
    W2SLSSpecs::Vector{CTExpr} = [Symbol("")]
    Z2SLSSpecs::Vector{CTExpr} = [parse("noticeLag1+0")]

    W2SLSNames::Vector{Vector{Symbol}} = [[:intercept]]
    X2SLSNames::Vector{Vector{Symbol}}  = [[:negLFlows]]
    Z2SLSNames::Vector{Vector{Symbol}}  = [[:noticeLag1]]

    #####put additional 2SLS specifications here
    push!(W2SLSSpecs, parse("fDate"))
    push!(W2SLSNames, [:intercept])

    push!(W2SLSSpecs, parse("main_strategy"))
    push!(W2SLSNames, [:intercept])

    push!(W2SLSSpecs, parse("$(vec2String(performanceLS.performanceLags, "+"))"))
    push!(W2SLSNames, [:intercept; performanceLS.performanceLags])

    push!(W2SLSSpecs, parse("$(vec2String(performanceLS.performanceLags, "+"))" *
        " + fDate"))
    push!(W2SLSNames, [:intercept; performanceLS.performanceLags])

    push!(W2SLSSpecs, parse("$(vec2String(performanceLS.performanceLags, "+"))" *
        " + main_strategy"))
    push!(W2SLSNames, [:intercept; performanceLS.performanceLags])

    push!(W2SLSSpecs, parse("$(vec2String(performanceLS.performanceLags, "+"))" *
        " + fDate + main_strategy"))
    push!(W2SLSNames, [:intercept; performanceLS.performanceLags])

    if RUN_2SLS_INTERACTIONS #needs too much memory. Would need a mem-mapped array or simialr solution to do this.
        push!(W2SLSSpecs, parse("fDate * main_strategy"))
        push!(W2SLSNames, [:intercept])

        push!(W2SLSSpecs, parse("$(vec2String(performanceLS.performanceLags, "+"))" *
          " + fDate * main_strategy"))
        push!(W2SLSNames, [:intercept; performanceLS.performanceLags])
    end

    ################# run the 2SLS models

    models2SLS::Vector{CT2SLS} = Vector{CT2SLS}()

    for iX ∈ 1:length(X2SLSSpecs), iW ∈ 1:length(W2SLSSpecs), iZ ∈ 1:length(Z2SLSSpecs)
        if verbose #print the specification to run
            println("Spec X: $(X2SLSSpecs[iX])
            Spec W: $(W2SLSSpecs[iW])
            Spec Z: $(Z2SLSSpecs[iZ])")
        end

        #build the regression model
        push!(models2SLS, CT2SLS(HFData, X2SLSSpecs[iX], W2SLSSpecs[iW],
            :performance, Z2SLSSpecs[iZ], XNames=X2SLSNames[iX],
            WNames = W2SLSNames[iW], YName = :performance, ZNames = Z2SLSNames[iZ]))
        println(models2SLS[end].ZNames)
        println(models2SLS[end].Π1)

        if verbose
            println("IV Model Specification:
                X: $(X2SLSSpecs[iX])
                W: $(W2SLSSpecs[iW])
                Z: $(Z2SLSSpecs[iZ])")
            println("WNames: $(W2SLSNames[iW])")
            println("β: $(models2SLS[i].δ2[1:min(10,end)])")
            println("σ: $(sqrt.(diag(getHomoskedΣ!(models2SLS[i])))[1:min(10,end)])")
            println("Data points: $(size(models2SLS[i].X,1))\n")
        end

    end
##################IO Code for 2SLS

    descRowNames2SLS::Vector{String} = ["Performance Lags", "Date Fixed Effects",
        "Strategy Fixed Effects", "N ('000s)"]

    #need to print the descriptive rows
    descContent2SLS::Vector{Vector{String}} =
        #[fill("",length(modelsLM)) for i∈ 1:length(descRowNamesLM)]
        [Vector{String}(length(models2SLS))for i∈ 1:length(descRowNames2SLS)]

    #this builds the descriptive rows. There is Probably a better way to do this,
    #but its fairly specific to the project.
    for i ∈ 1:length(W2SLSSpecs)
        descContent2SLS[1][i] = "$(contains(W2SLSSpecs[i], "performanceLag")?"X":"")"
        descContent2SLS[2][i] = "$(contains(W2SLSSpecs[i], "fDate")?"X":"")"
        descContent2SLS[3][i] = "$(contains(W2SLSSpecs[i], "main_strategy")?"X":"")"
        descContent2SLS[4][i] = "$(round(models2SLS[i].N/1000.0,1))"
    end
    #note this only works for a single focal var
    stage12SLS::Vector{CTLM} = [get1stStage(m)[1] for m ∈ models2SLS]
    TableText1stStage::String = texTable(stage12SLS, getHomoskedΣ!, [:intercept; :noticeLag1],
        caption = "Performance Flows IV 1st Stage (Focal: Outflows)",
        colNames = [["($i)" for i∈1:length(models2SLS)]],
        contentRowNames = ["intercept", "Notifications (frac of fund x1000)"],
        descRowNames = descRowNames2SLS,
        descContent = descContent2SLS,
        decimalDigits = 3,
        columnSepPt = -20,
        scaling = [1000.0,1000.0],
        clearMem = USE_AGGRESSIVE_GC)

    TableText2SLS::String = texTable(models2SLS, getHomoskedΣ!, [:intercept; :negLFlows],
        caption = "Performance Flows 2SLS Specifications",
        colNames = [["Z=1 Month Notification of Flows"],["($i)" for i∈1:length(models2SLS)]],
        contentRowNames = ["intercept", "outflows/AUM x 1000"],
        descRowNames = descRowNames2SLS,
        descContent = descContent2SLS,
        decimalDigits = 3,
        columnSepPt = -20,
        scaling = [1.0,1000.0],
        widthColNames = [[length(W2SLSSpecs)],ones(Int,length(W2SLSSpecs))],
        clearMem = USE_AGGRESSIVE_GC)

    writeTables2File([TableTextLM, TableText1stStage, TableText2SLS], HEADER_NAME, FOOTER_NAME, path=RESULTS_PATH,
        outName = "RegTables_$FlowsLaggingProcedure.tex")
    #=plOut::Plot = plot(melt(HFByDate[:,[:date,
    #:totFlows,  :totRedemp, :noticeLag1,
    :totRedempNotice,:noticeLag1,:noticeLag2]],:date),
        x=:date,y=:value,color=:variable,Geom.line,
        Guide.title("Total Fund Flows over Time"),
        Guide.xlabel("Date"),Guide.ylabel("Flows"))
      draw(PNG("$(DATA_FILE_NAME)_$(FlowsLaggingProcedure)_flows.png", 7.5inch, 6inch), plOut)
      =#
    close(iStream)
end

function verify2R(;FlowsLaggingProcedure::Char = 'U')

    #read the seria
    iStream::IOStream =
        open("$DATA_PATH\\$(DATA_FILE_NAME)_$(FlowsLaggingProcedure)_clean.jls")
    HFData::DataFrame = deserialize(iStream)
    close(iStream)

    #write the csv
    writetable("$R_TEST_PATH\\RHFData.csv",HFData)
end

@time begin
#PreProcessHFData(FlowsLaggingProcedure='C', newBinary=false)
AnalyzeAggregates(FlowsLaggingProcedure='C', verbose=false)
verify2R(FlowsLaggingProcedure='C')
end

##Full Pass:

#=@time for c ∈ ['C', 'R', 'L', 'U']
  PreProcessHFData(FlowsLaggingProcedure=c, newBinary=true)
  AnalyzeAggregates(FlowsLaggingProcedure=c)
end=#

end

max_stars_repo_path: Old/SummerProject17/SAVE- HF Data Functions-preNLRedux.jl


text:
struct Terminator end

mutable struct InPort <: AbstractPort
    line::Union{AbstractSignalLine,Terminator}
    blk::AbstractBlock
    InPort() = new(Terminator())
end

mutable struct OutPort <: AbstractPort
#     line::Union{AbstractSignalLine,Terminator}
    blk::AbstractBlock
    OutPort() = new() #, Terminator())
end

"""
IO
"""

function Base.show(io::IO, port::InPort)
    println(io, "InPort()")
end

function Base.show(io::IO, port::OutPort)
    println(io, "OutPort()")
end

"""
Add ports
"""

function _addport!(blk::AbstractBlock, port::InPort)
    port.blk = blk
    push!(blk.inport, port)
    nothing
end

function _addport!(blk::AbstractBlock, port::OutPort)
    port.blk = blk
    push!(blk.outport, port)
    nothing
end

macro createblock(expr, innum, outnum)
    ex = quote
        b = $(expr)
        for t = 1:$innum
            _addport!(b, InPort())
        end
        for t = 1:$outnum
            _addport!(b, OutPort())
        end
        b
    end
    esc(ex)
end

"""
toexpr

This is used to make the expression.
The expression is used for making both system equation and output equation.
"""

function _toexpr(port::Terminator)
    :terminator
end

function _toexpr(port::InPort)
    _toexpr(port.line)
end

function _toexpr(port::Vector{InPort})
    _toexpr(port[1].line)
end

function _toexpr(port::OutPort)
    _toexpr(port.blk)
end

function _toexpr(port::Vector{OutPort})
    _toexpr(port[1].blk)
end

ext: jl


text:
using Test

m = page_model()
include("../src/components/N2Oforcing.jl")

add_comp!(m, n2oforcing)

set_param!(m, :n2oforcing, :c_N2Oconcentration, readpagedata(m,"test/validationdata/c_n2oconcentration.csv"))
set_param!(m, :n2oforcing, :c_CH4concentration, readpagedata(m,"test/validationdata/c_ch4concentration.csv"))

##running Model
run(m)

forcing=m[:n2oforcing,:f_N2Oforcing]
forcing_compare=readpagedata(m,"test/validationdata/f_n2oforcing.csv")

@test forcing ≈ forcing_compare rtol=1e-3


text:
using RealSense

include("rs-utils.jl")

err = Ref{Ptr{rs2_error}}(C_NULL)

ctx = rs2_create_context(RS2_API_VERSION, err)
checkerror(err)

device = get_a_realsense_device(ctx)
print_device_information(device)

@info "Please choose a sensor and then a stream that will be used as the origin of extrinsic transformation:"
from_sensor = get_a_sensor_from_a_device(device)
from = choose_a_streaming_profile(from_sensor)

@info "Please choose a sensor and then a stream that will be used as the target of extrinsic transformation:"
to_sensor = get_a_sensor_from_a_device(device)
to = choose_a_streaming_profile(to_sensor)

extrinsicsRef = Ref{rs2_extrinsics}(rs2_extrinsics(ntuple(x->0,Val{9}),(0,0,0)))
rs2_get_extrinsics(from, to, extrinsicsRef, err)
checkerror(err)

extrinsics = extrinsicsRef[]
@info "Rotation Matrix: "
@info extrinsics.rotation[1:3]
@info extrinsics.rotation[4:6]
@info extrinsics.rotation[7:9]
@info "Translation Vector: " extrinsics.translation


text:
@testset "Check identifiability of `ODESystem` object" begin
    @parameters a01 a21 a12
    @variables t x0(t) x1(t) y1(t) [output=true]
    D = Differential(t)

    eqs = [
        D(x0) ~ -(a01 + a21) * x0 + a12 * x1,
        D(x1) ~ a21 * x0 - a12 * x1,
        y1 ~ x0
    ]
    de = ODESystem(eqs, t, name=:Test)

    correct = Dict(a01 => :nonidentifiable, a21 => :nonidentifiable, a12 => :nonidentifiable)

    @test isequal(correct, assess_identifiability(de; measured_quantities=[y1~x0]))
    # --------------------------------------------------------------------------
    @parameters a01 a21 a12
    @variables t x0(t) x1(t) y1(t) [output=true]
    D = Differential(t)

    eqs = [
        D(x0) ~ -(a01 + a21) * x0 + a12 * x1,
        D(x1) ~ a21 * x0 - a12 * x1,
        y1 ~ x0
    ]
    de = ODESystem(eqs, t, name=:Test)

    correct = Dict(a01 => :nonidentifiable, a21 => :nonidentifiable, a12 => :nonidentifiable)

    @test isequal(correct, assess_identifiability(de))

    # --------------------------------------------------------------------------

    @parameters a01 a21 a12
    @variables t x0(t) x1(t) y1(t) [output=true]
    D = Differential(t)

    eqs = [
        D(x0) ~ -(a01 + a21) * x0 + a12 * x1,
        D(x1) ~ a21 * x0 - a12 * x1,
        y1 ~ x0
    ]
    de = ODESystem(eqs, t, name=:Test)
    funcs_to_check = [a01, a21, a12, a01 * a12, a01 + a12 + a21]
    correct = Dict(a12 => :nonidentifiable, a01 + a12 + a21 => :globally, a01 * a12 => :globally, a21 => :nonidentifiable, a01 => :nonidentifiable)
    @test isequal(correct, assess_identifiability(de; funcs_to_check = funcs_to_check))

    # --------------------------------------------------------------------------

    @parameters a01 a21 a12
    @variables t x0(t) x1(t) y1(t)
    D = Differential(t)

    eqs = [
        D(x0) ~ -(a01 + a21) * x0 + a12 * x1,
        D(x1) ~ a21 * x0 - a12 * x1
    ]
    measured_quantities = [y1 ~ x0]
    de = ODESystem(eqs, t, name=:Test)
    funcs_to_check = [a01, a21, a12, a01 * a12, a01 + a12 + a21]
    correct = Dict(a12 => :nonidentifiable, a01 + a12 + a21 => :globally, a01 * a12 => :globally, a21 => :nonidentifiable, a01 => :nonidentifiable)
    @test isequal(correct, assess_identifiability(de; measured_quantities=measured_quantities, funcs_to_check=funcs_to_check))

    # --------------------------------------------------------------------------
    @parameters mu bi bw a xi gm k
    @variables t S(t) I(t) W(t) R(t) y(t)

    eqs = [
        D(S) ~ mu - bi * S * I - bw * S * W - mu * S + a * R,
        D(I) ~ bw * S * W + bi * S * I - (gm + mu) * I,
        D(W) ~ xi * (I - W),
        D(R) ~ gm * I - (mu + a) * R
    ]
    de = ODESystem(eqs, t, name=:TestSIWR)
    measured_quantities = [y ~ k * I]
    # check all parameters (default)
    @test isequal(true, all(assess_local_identifiability(de; measured_quantities=measured_quantities)))

    # check specific parameters
    funcs_to_check = [mu, bi, bw, a, xi, gm, gm + mu, k, S, I, W, R]
    correct = [true for _ in funcs_to_check]
    @test isequal(correct, assess_local_identifiability(de; measured_quantities=measured_quantities, funcs_to_check=funcs_to_check))

    # checking ME identifiability
    funcs_to_check = [mu, bi, bw, a, xi, gm, gm + mu, k]
    correct = [true for _ in funcs_to_check]
    @test isequal((correct, 1), assess_local_identifiability(de;  measured_quantities=measured_quantities, funcs_to_check=funcs_to_check, p=0.99, type=:ME))

    # --------------------------------------------------------------------------
    @parameters mu bi bw a xi gm k
    @variables t S(t) I(t) W(t) R(t) y(t) [output=true]

    eqs = [
        D(S) ~ mu - bi * S * I - bw * S * W - mu * S + a * R,
        D(I) ~ bw * S * W + bi * S * I - (gm + mu) * I,
        D(W) ~ xi * (I - W),
        D(R) ~ gm * I - (mu + a) * R,
        y ~ k * I
    ]
    de = ODESystem(eqs, t, name=:TestSIWR)
    # check all parameters (default)
    @test isequal(true, all(assess_local_identifiability(de)))

    @test isequal(true, all(assess_local_identifiability(de; measured_quantities=[y~k*I])))

    # check specific parameters
    funcs_to_check = [mu, bi, bw, a, xi, gm, gm + mu, k, S, I, W, R]
    correct = [true for _ in funcs_to_check]
    @test isequal(correct, assess_local_identifiability(de; funcs_to_check=funcs_to_check))

    # checking ME identifiability
    funcs_to_check = [mu, bi, bw, a, xi, gm, gm + mu, k]
    correct = [true for _ in funcs_to_check]
    @test isequal((correct, 1), assess_local_identifiability(de; funcs_to_check=funcs_to_check, p=0.99, type=:ME))

    # --------------------------------------------------------------------------
    @parameters mu bi bw a xi gm k
    @variables t S(t) I(t) W(t) R(t) y(t)

    eqs = [
        D(S) ~ 2.0 * mu - bi * S * I - bw * S * W - mu * S + a * R,
        D(I) ~ bw * S * W + bi * S * I - (gm + mu) * I,
        D(W) ~ xi * (I - 0.6 * W),
        D(R) ~ gm * I - (mu + a) * R,
    ]
    de = ODESystem(eqs, t, name=:TestSIWR)
    measured_quantities = [y ~ 1.57 * I * k]
    funcs_to_check = [mu, bi, bw, a, xi, gm, mu, gm + mu, k, S, I, W, R]
    correct = [true for _ in funcs_to_check]
    @test isequal(correct, assess_local_identifiability(de; measured_quantities=measured_quantities, funcs_to_check=funcs_to_check))

    # checking ME identifiability
    funcs_to_check = [bi, bw, a, xi, gm, mu, gm + mu, k]
    correct = [true for _ in funcs_to_check]
    @test isequal((correct, 1), assess_local_identifiability(de; measured_quantities=measured_quantities, funcs_to_check=funcs_to_check, p=0.99, type=:ME))

    # ----------

    @parameters a01 a21 a12
    @variables t x0(t) x1(t) y1(t)
    D = Differential(t)
    using SpecialFunctions

    eqs = [
        D(x0) ~ -(a01 + a21) * SpecialFunctions.erfc(x0) + a12 * x1,
        D(x1) ~ a21 * x0 - a12 * x1
    ]

    de = ODESystem(eqs, t, name=:Test)
    measured_quantities = [y1 ~ x0]
    funcs_to_check = [a01, a21, a12, a01 * a12, a01 + a12 + a21]
    correct = Dict(a01 => :nonidentifiable, a21 => :nonidentifiable, a12 => :nonidentifiable, a01 * a12 => :globally, a01 + a12 + a21 => :globally)
    @test_throws ArgumentError assess_identifiability(de; measured_quantities=measured_quantities,  funcs_to_check=funcs_to_check)
    # ----------
    @parameters a b c
    @variables t x1(t) x2(t) y(t)
    D = Differential(t)

    eqs = [
        D(x1) ~ -a * x1 + x2 * b / (x1 + b / (c^2 - x2)),
        D(x2) ~ x2 * c^2 + x1
    ]
    de = ODESystem(eqs, t, name=:Test)
    measured_quantities = [y ~ x2]
    correct = Dict(a => :globally, b => :globally, c => :locally)
    to_check = [a, b, c]
    @test isequal(correct, assess_identifiability(de; measured_quantities=measured_quantities, funcs_to_check=to_check))
end


text:
# Interaction with Mouse
# Original implementation is here:
# https://simondanisch.github.io/ReferenceImages/gallery/interaction_with_mouse/index.html

using Makie

function main()
    scene = Scene(raw = true, camera = cam2d!, resolution = (500, 500))
    r = LinRange(0, 3, 4)
    the_time = Node(time())
    last_open = false

    @async while true
        global last_open
        the_time[] = time()
        last_open && !isopen(scene) && break
        last_open && isopen(scene)
        sleep(1 / 30)
    end

    pos = lift(scene.events.mouseposition, the_time) do mpos, t
        mouse = to_world(scene, Point2f0(mpos))
        Point2f0.([mouse, mouse])
    end

    scatter!(scene, pos)

    RecordEvents(hbox(scene), "output")
end

main()

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
using Testfast

# Run once before all tests
@before_all begin
  before_all_test = 0
  setup_test = 0
  teardown_test = 0
  after_all_test = 0

  before_all_test += 1
end

# Run before each test
@setup begin
  setup_test += 1
end

# Run after each test
@teardown begin
  teardown_test += 1
end

# Run once after all tests
@after_all begin
  after_all_test += 1
end

Testfast.run!()


text:
module ACWT
export
    # Autocorrelation wavelet transform
    acdwt,
    acdwt!,
    acwpt,
    acwpt!,
    acwpd,
    acwpd!,
    # Inverse autocorrelation wavelet transform
    iacdwt,
    iacdwt!,
    iacwpt,
    iacwpt!,
    iacwpd,
    iacwpd!,
    # Transform on all signals
    acdwtall,
    acwptall,
    acwpdall,
    iacdwtall,
    iacwptall,
    iacwpdall

using ..Utils
using LinearAlgebra, Wavelets

# ========== Autocorrelation Discrete Wavelet Transform ==========
@doc raw"""
	acdwt(x, wt[, L])

Performs a discrete autocorrelation wavelet transform for a given signal `x`.
The signal can be 1D or 2D. The wavelet type `wt` determines the transform type.
Refer to Wavelet.jl for a list of available methods.

# Arguments
- `x::AbstractArray{T} where T<:Number`: Original signal, preferably of size 2ᴷ where ``K
    \in \mathbb{N}``.
- `wt::OrthoFilter`: Orthogonal wavelet filter.
- `L::Integer`: (Default: `maxtransformlevels(x)`) Number of levels of decomposition.

# Returns
`::Array{T}`: Output from ACDWT on `x`.

# Examples
```julia
using Wavelets, WaveletsExt

# Setup
x = generatesignals(:heavysine)
wt = wavelet(WT.haar)

# ACDWT
acdwt(x, wt)
acdwt(x, wt, 4) # level 4 decomposition
```

**See also:** [`acdwt_step`](@ref), [`iacdwt`](@ref), [`acdwt!`](@ref)
"""
function acdwt(x::AbstractArray{T},
               wt::OrthoFilter,
               L::Integer = maxtransformlevels(x)) where T<:Number
    # Sanity check
    @assert 1 ≤ ndims(x) ≤ 2
    @assert L ≤ maxtransformlevels(x) || throw(ArgumentError("Too many transform levels (length(x) < 2^L"))
    @assert L ≥ 1 || throw(ArgumentError("L must be ≥ 1"))
    # Setup
    sz = size(x)
    N = ndims(x)+1
    k = N==2 ? L+1 : 3*L+1
    xw = Array{T,N}(undef, (sz...,k))
    # Compute transforms
    acdwt!(xw, x, wt, L)
    return xw
end

# 1D ACDWT without allocation
@doc raw"""
    acdwt!(xw, x, wt, L)

Same as `acdwt` but without array allocation.

# Arguments
- `xw::AbstractArray{T}`: An allocated array of dimension `(n,L+1)` or `(n,m,3L+1)` to write
  the outputs of `x` onto.
- `x::AbstractArray{T} where T<:Number`: Original signal, preferably of size 2ᴷ where ``K
  \in \mathbb{N}``.
- `wt::OrthoFilter`: Orthogonal wavelet filter.
- `L::Integer`: (Default: `maxtransformlevels(x)`) Number of levels of decomposition.

# Returns
`xw::Array{T,2}`: Output from ACDWT on `x`.

# Examples
```julia
using Wavelets, WaveletsExt

# Setup
x = generatesignals(:heavysine, 7)
wt = wavelet(WT.haar)

# ACDWT
xw = Matrix{Float64}(undef, (128,5))
acdwt!(xw, x, wt, 4)
```

**See also:** [`acdwt`](@ref)
"""
function acdwt!(xw::AbstractArray{T,2},
                x::AbstractVector{T},
                wt::OrthoFilter,
                L::Integer = maxtransformlevels(x)) where T<:Number
    # Sanity check
    @assert L <= maxtransformlevels(x) || throw(ArgumentError("Too many transform levels (length(x) < 2^L"))
    @assert L >= 1 || throw(ArgumentError("L must be >= 1"))

    # Setup
    Pmf, Qmf = make_acreverseqmfpair(wt)
    xw[:,end] = x

    # ACDWT
    for d in 0:(L-1)
        @inbounds v = xw[:, L-d+1]       # Parent node
        w₁ = @view xw[:, L-d]            # Scaling coefficients
        w₂ = @view xw[:, L-d+1]          # Detail coefficients
        @inbounds acdwt_step!(w₁, w₂, v, d, Qmf, Pmf)
    end
    return xw
end
# 2D ACDWT without allocation
function acdwt!(xw::AbstractArray{T,3},
                x::AbstractMatrix{T},
                wt::OrthoFilter,
                L::Integer = maxtransformlevels(x)) where T<:Number
    # Sanity check
    @assert L ≤ maxtransformlevels(x) ||
        throw(ArgumentError("Too many transform levels"))
    @assert L ≥ 1 || throw(ArgumentError("L must be ≥ 1"))
    @assert size(xw,3) == 3*L+1

    # Setup
    n, m = size(x)
    g, h = make_acreverseqmfpair(wt)
    xw[:,:,end] = x
    temp = Array{T,3}(undef, (n,m,2))

    # ACDWT
    for d in 0:(L-1)
        @inbounds v = xw[:,:, 3*(L-d)+1]       # Parent node
        @inbounds w₁ = @view xw[:,:, 3*(L-d)-2]          # Scaling + Scaling coefficients
        @inbounds w₂ = @view xw[:,:, 3*(L-d)-1]          # Detail + Scaling coefficients
        @inbounds w₃ = @view xw[:,:, 3*(L-d)]            # Scaling + Detail coefficients
        @inbounds w₄ = @view xw[:,:, 3*(L-d)+1]          # Detail + Detail coefficients
        @inbounds acdwt_step!(w₁, w₂, w₃, w₄, v, d, h, g, temp)
    end
    return xw
end

## 2D ##
"""
    hacdwt(x, wt[, L=maxtransformlevels(x,2)])

Computes the column-wise discrete autocorrelation transform coeficients for 2D signals.

**See also:** [`vacdwt`](@ref)
"""
# function hacdwt(x::AbstractArray{T,2},
#                wt::OrthoFilter,
#                L::Integer=maxtransformlevels(x,2)) where T<:Number
#     nrow, ncol = size(x)
#     W = Array{T,3}(undef,nrow,L+1,ncol)
#     for i in 1:ncol
#         @inbounds W[:,:,i] = acdwt(x[:,i],wt,L)
#     end
#     return W
# end

"""
    vacdwt(x, wt[, L=maxtransformlevels(x)])

Computes the row-wise discrete autocorrelation transform coeficients for 2D signals.

**See also:** [`hacdwt`](@ref)
"""
# function vacdwt(x::AbstractArray{T,2},
#                wt::OrthoFilter,
#                L::Integer=maxtransformlevels(x,1)) where T<:Number
#     nrow, ncol = size(x)
#     W = Array{T,3}(undef,ncol,L+1,nrow)
#     for i in 1:nrow
#         W[:,:,i] = acdwt(x[i,:],wt,L)
#     end
#     return W
# end

# function acdwt(x::AbstractArray{T,2}, wt::OrthoFilter,
#               Lrow::Integer=maxtransformlevels(x,1),
#               Lcol::Integer=maxtransformlevels(x,2)) where T<:Number
#     nrow, ncol = size(x)
#     W3d = hacdwt(x,wt,Lcol)
#     W4d = Array{T,4}(undef,Lcol+1,ncol,Lrow+1,nrow)
#     for i in 1:Lcol+1
#         @inbounds W4d[i,:,:,:] = vacdwt(W3d[:,i,:],wt,Lrow)
#     end
#     W4d = permutedims(W4d, [4,2,3,1])
#     return W4d
# end

"""
    iacdwt(xw[, wt])

Performs the inverse autocorrelation discrete wavelet transform. Can be used for both the 1D
and 2D case.

!!! note
    The inverse autocorrelation transform does not require any wavelet filter, but an
    optional `wt` positional argument is included for the standardization of syntax with
    `dwt` and `sdwt`, but is ignored during the reconstruction of signals.

# Arguments
- `xw::AbstractArray{T} where T<:Number`: ACDWT-transformed array.
- `wt::Union{OrthoFilter, Nothing}`: (Default: `nothing`) Orthogonal wavelet filter.

# Returns
`::Array{T}`: Inverse transformed signal.

# Examples
```julia
using Wavelets, WaveletsExt

# Setup
x = generatesignals(:heavysine)
wt = wavelet(WT.haar)

# ACDWT
xw = acdwt(x, wt)

# IACDWT
x̃ = iacdwt(xw)
```

**See also:** [`acdwt`](@ref), [`iacdwt_step!`](@ref)
"""
function iacdwt(xw::AbstractArray{T}, wt::Union{OrthoFilter,Nothing} = nothing) where
                T<:Number
    @assert ndims(xw) == 2 || ndims(xw) == 3
    sz = size(xw)[1:(end-1)]
    N = ndims(xw)-1
    x = Array{T,N}(undef, sz)
    iacdwt!(x, xw)
    return x
end

"""
    iacdwt!(x, xw[, wt])

Similar to `iacdwt` but without array allocation.

# Arguments
- `x::AbstractArray{T} where T<:Number` or `x::AbstractArray{T,2} where T<:Number`:
  Allocation for reconstructed signal.
- `xw::AbstractArray{T} where T<:Number` or `xw::AbstractArray{T,4}`: ACDWT-transformed
  array.
- `wt::Union{OrthoFilter, Nothing}`: (Default: `nothing`) Orthogonal wavelet filter.

# Returns
`x::Array{T}`: Inverse transformed signal.

# Examples
```julia
using Wavelets, WaveletsExt

# Setup
x = generatesignals(:heavysine)
wt = wavelet(WT.haar)

# ACDWT
xw = acdwt(x, wt)

# IACDWT
x̃ = similar(x)
iacdwt!(x̃, xw)
```

**See also:** [`iacdwt`](@ref), [`acdwt`](@ref)
"""
function iacdwt!(x::AbstractVector{T},
                 xw::AbstractArray{T,2},
                 wt::Union{OrthoFilter,Nothing} = nothing) where T<:Number
    # Setup
    _, k = size(xw)
    L = k-1
    for i in eachindex(x)
        @inbounds x[i] = xw[i,1]
    end

    # IACDWT
    for d in reverse(0:(L-1))
        w₁ = copy(x)
        @inbounds w₂ = @view xw[:,L-d+1]
        @inbounds iacdwt_step!(x, w₁, w₂)
    end
    return x
end
# 2D IACDWT
function iacdwt!(x::AbstractMatrix{T},
                 xw::AbstractArray{T,3},
                 wt::Union{OrthoFilter, Nothing} = nothing) where
                 T<:Number
    # Sanity check
    n, m, k = size(xw)
    L = (k-1) ÷ 3

    # Setup
    temp = Array{T,3}(undef, (n,m,2))
    for i in eachindex(x)
        @inbounds x[i] = xw[i]
    end

    # ISDWT
    for d in reverse(0:(L-1))
        w₁ = copy(x)
        @inbounds w₂ = @view xw[:,:,3*(L-d)-1]
        @inbounds w₃ = @view xw[:,:,3*(L-d)]
        @inbounds w₄ = @view xw[:,:,3*(L-d)+1]
        @inbounds iacdwt_step!(x, w₁, w₂, w₃, w₄, temp)
    end
    return x
end
# function iacdwt!(x::AbstractArray{T,2},
#                  xw::AbstractArray{T,4},
#                  wt::Union{OrthoFilter,Nothing} = nothing) where T<:Number
#     @assert size(x,1) == size(xw,1)
#     @assert size(x,2) == size(xw,2)

#     nrow, ncol, _, Lcol = size(xw)
#     W4d = permutedims(xw,[4,2,3,1])
#     W3d = Array{T,3}(undef, nrow, Lcol, ncol)
#     for i in 1:Lcol
#         for j in 1:nrow
#             @inbounds W3d[j,i,:] = iacdwt(W4d[i,:,:,j])
#         end
#     end
#     @views for (i, xᵢ) in enumerate(eachcol(x))
#         @inbounds iacdwt!(xᵢ, W3d[:,:,i])
#     end
#     return x
# end

# ========== Autocorrelation Wavelet Packet Transform ==========
@doc raw"""
    acwpt(x, wt[, L])

Computes `L` levels of autocorrelation wavelet packet transforms (ACWPT) on `x`.

# Arguments
- `x::AbstractArray{T} where T<:Number`: Original signal, preferably of size 2ᴷ or (2ᴷ,2ᴹ)
  where ``K, M \in \mathbb{N}``.
- `wt::OrthoFilter`: Orthogonal wavelet filter.
- `L::Integer`: (Default: `maxtransformlevels(x)`) Number of levels of decomposition.

# Returns
`::Array{T}`: Output from ACWPT on `x`.

# Examples
```julia
using Wavelets, WaveletsExt

# Setup
x = generatesignals(:heavysine)
wt = wavelet(WT.haar)

# ACWPT
xw = acwpt(x, wt)
```

**See also:** [`iacwpt`](@ref), [`acdwt`](@ref), [`acwpd`](@ref)
"""
function acwpt(x::AbstractArray{T},
               wt::OrthoFilter,
               L::Integer = maxtransformlevels(x)) where T<:Number
    # Sanity check
    @assert 1 ≤ ndims(x) ≤ 2
    @assert L ≤ maxtransformlevels(x) ||
            throw(ArgumentError("Too many transform levels (length(x) < 2ᴸ)"))
    @assert L ≥ 1 || throw(ArgumentError("L must be ≥ 1"))
    # Setup
    sz = size(x)
    N = ndims(x)+1
    k = N==2 ? 1<<L : 1<<(2*L)
    xw = Array{T,N}(undef, (sz..., k))
    # Transform
    acwpt!(xw, x, wt, L)
    return xw
end

@doc raw"""
    acwpt!(xw, x, wt[, L])

Same as `acwpt` but without array allocation.

# Arguments
- `xw::AbstractArray{T} where T<:Number`: Allocation for transformed signal.
- `x::AbstractArray{T} where T<:Number`: Original signal, preferably of size 2ᴷ or (2ᴷ,2ᴹ)
  where ``K, M \in \mathbb{N}``.
- `wt::OrthoFilter`: Orthogonal wavelet filter.
- `L::Integer`: (Default: `maxtransformlevels(x)`) Number of levels of decomposition.

# Returns
`xw::Array{T}`: Output from ACWPT on `x`.

# Examples
```julia
using Wavelets, WaveletsExt

# Setup
x = generatesignals(:heavysine)
wt = wavelet(WT.haar)

# ACWPT
xw = Array{Float64,2}(undef, (128,128))
acwpt!(xw, x, wt)
```

**See also:** [`acwpt`](@ref), [`iacwpt`](@ref)
"""
function acwpt!(xw::AbstractArray{T,2},
                x::AbstractVector{T},
                wt::OrthoFilter,
                L::Integer = maxtransformlevels(x)) where T<:Number
    # Sanity check
    @assert L ≤ maxtransformlevels(x) ||
            throw(ArgumentError("Too many transform levels (length(x) < 2ᴸ)"))
    @assert L ≥ 1 || throw(ArgumentError("L must be ≥ 1"))
    @assert size(x,1) == size(xw,1)
    @assert size(xw,2) == 1<<L

    # Setup
    Pmf, Qmf = make_acreverseqmfpair(wt)
    for i in eachindex(x)
        @inbounds xw[i,1] = x[i]
    end

    # ACWPT for L levels
    for d in 0:(L-1)
        nn = 1<<d
        for b in 0:(nn-1)
            np = (1<<L)÷nn
            nc = np÷2
            j₁ = (2*b)*nc + 1   # Parent and (scaling) child index
            j₂ = (2*b+1)*nc + 1 # Detail child index
            v = xw[:,j₁]
            w₁ = @view xw[:,j₁]
            w₂ = @view xw[:,j₂]
            # Overwrite output of ACDWT directly onto xw
            @inbounds acdwt_step!(w₁, w₂, v, d, Qmf, Pmf)
        end
    end
    return xw
end
# 2D ACWPT
function acwpt!(xw::AbstractArray{T,3},
                x::AbstractMatrix{T},
                wt::OrthoFilter,
                L::Integer = maxtransformlevels(x)) where T<:Number
    @assert L ≤ maxtransformlevels(x) ||
        throw(ArgumentError("Too many transform levels."))
    @assert L ≥ 1 || throw(ArgumentError("L must be ≥ 1."))
    @assert size(xw,3) == 4^L
    @assert size(xw,1) == size(x,1)
    @assert size(xw,2) == size(x,2)

    # Setup
    n, m = size(x)
    g, h = make_acreverseqmfpair(wt)
    temp = Array{T,3}(undef, (n,m,2))
    for i in eachindex(x)
        @inbounds xw[i] = x[i]
    end

    # ACWPT for L levels
    for d in 0:(L-1)
        nn = 4^d            # Number of nodes at current level
        for b in 0:(nn-1)
            np = (4^L)÷nn       # Number of slices occupied by parent
            nc = np÷4           # Number of slices occupied by children
            j₁ = (4*b)*nc + 1   # Parent and scaling+scaling child index
            j₂ = (4*b+1)*nc + 1 # Scaling+detail child index
            j₃ = (4*b+2)*nc + 1 # Detail+scaling child index
            j₄ = (4*b+3)*nc + 1 # Detail+detail child index
            v = xw[:,:,j₁]
            w₁ = @view xw[:,:,j₁]
            w₂ = @view xw[:,:,j₂]
            w₃ = @view xw[:,:,j₃]
            w₄ = @view xw[:,:,j₄]
            # Overwrite output of SDWT directly onto xw
            @inbounds acdwt_step!(w₁, w₂, w₃, w₄, v, d, h, g, temp)
        end
    end
    return xw
end

"""
    iacwpt(xw[, wt])

Computes the inverse autocorrelation wavelet packet transform (IACWPT) on `xw`.

!!! note
    The inverse autocorrelation transform does not require any wavelet filter, but an
    optional `wt` positional argument is included for the standardization of syntax with
    `wpt` and `swpt`, but is ignored during the reconstruction of signals.

# Arguments
- `xw::AbstractArray{T} where T<:Number`: ACWPT-transformed array.
- `wt::Union{OrthoFilter, Nothing}`: (Default: `nothing`) Orthogonal wavelet filter.

# Returns
`::Array{T}`: Inverse transformed signal.

# Examples
```julia
using Wavelets, WaveletsExt

# Setup
x = generatesignals(:heavysine)
wt = wavelet(WT.haar)

# ACWPT
xw = acwpt(x, wt)

# IACWPT
x̃ = iacwpt(xw)
```

**See also:** [`iacdwt`](@ref), [`acwpt`](@ref)
"""
function iacwpt(xw::AbstractArray{T}, wt::Union{OrthoFilter, Nothing} = nothing) where
                T<:Number
    @assert 2 ≤ ndims(xw) ≤ 3
    # Setup
    sz = size(xw)[1:(end-1)]
    N = ndims(xw)-1
    x = Array{T,N}(undef, sz)
    # Transform
    iacwpt!(x, xw)
    return x
end

# 1D IACWPT without allocation
"""
    iacwpt!(x, xw[, wt])

Same as `iacwpt` but without array allocation.

# Arguments
- `x::AbstractArray{T} where T<:Number`: Allocated array for output.
- `xw::AbstractArray{T} where T<:Number`: ACWPD-transformed array.
- `wt::Union{OrthoFilter, Nothing}`: (Default: `nothing`) Orthogonal wavelet filter.

# Returns
`x::Array{T}`: Inverse transformed signal.

# Examples
```julia
using Wavelets, WaveletsExt

# Setup
x = generatesignals(:heavysine)
wt = wavelet(WT.haar)

# ACWPT
xw = acwpt(x, wt)

# IACWPT
x̂ = similar(x)
iacwpt!(x̂, xw)
```

**See also:** [`iacwpt`](@ref)
"""
function iacwpt!(x::AbstractVector{T},
                 xw::AbstractArray{T,2},
                 wt::Union{OrthoFilter,Nothing} = nothing) where T<:Number
    # Sanity check
    n, m = size(xw)
    @assert isdyadic(m) || throw(ArgumentError("Number of columns of xw is not dyadic."))
    @assert ndyadicscales(m) ≤ maxtransformlevels(n) ||
            throw(ArgumentError("Number of nodes in `xw` is more than possible number of nodes at any depth for signal of length `n`"))

    # Setup
    temp = copy(xw)                         # Temp. array to store intermediate outputs
    L = ndyadicscales(m)                    # Number of decompositions from xw

    # IACWPT
    for d in reverse(0:(L-1))
        nn = 1<<d
        for b in 0:(nn-1)
            np = (1<<L)÷nn
            nc = np÷2
            j₁ = (2*b)*nc + 1                           # Parent and (scaling) child index
            j₂ = (2*b+1)*nc + 1                         # Detail child index
            @inbounds v = d==0 ? x : @view temp[:,j₁]   # Parent node
            @inbounds w₁ = temp[:,j₁]                   # Scaling node
            @inbounds w₂ = @view temp[:,j₂]             # Detail node
            # Overwrite output of iSWPT directly onto temp
            @inbounds iacdwt_step!(v, w₁, w₂)
        end
    end
    return x
end
# 2D IACWPT
function iacwpt!(x::AbstractMatrix{T},
                 xw::AbstractArray{T,3},
                 wt::Union{OrthoFilter,Nothing} = nothing) where
                 T<:Number
    # Sanity check
    n, m, k = size(xw)
    @assert isinteger(log(4,k)) || throw(ArgumentError("Size of dimension 3 is not a power of 4."))
    @assert size(x,1) == n
    @assert size(x,2) == m
    @assert log(4,k) ≤ maxtransformlevels(x)

    # Setup
    xwₜ = copy(xw)
    temp = Array{T,3}(undef, (n,m,2))       # Temporary array
    L = log(4,k) |> Int                     # Number of decomposition levels of xw

    # ISWPT
    for d in reverse(0:(L-1))
        nn = 4^d
        for b in 0:(nn-1)
            np = (4^L)÷nn                   # Number of slices occupied by parent
            nc = np÷4                       # Number of slices occupied by children
            j₁ = (4*b)*nc + 1               # Parent and scaling+scaling child index
            j₂ = (4*b+1)*nc + 1             # Scaling+detail child index
            j₃ = (4*b+2)*nc + 1             # Detail+scaling child index
            j₄ = (4*b+3)*nc + 1             # Detail+detail child index
            v = d==0 ? x : @view xwₜ[:,:,j₁]
            w₁ = xwₜ[:,:,j₁]
            w₂ = @view xwₜ[:,:,j₂]
            w₃ = @view xwₜ[:,:,j₃]
            w₄ = @view xwₜ[:,:,j₄]
            # Overwrite output of iSWPT directly onto temp
            @inbounds iacdwt_step!(v, w₁, w₂, w₃, w₄, temp)
        end
    end
    return x
end

# ========== Autocorrelation Wavelet Packet Decomposition ==========
@doc raw"""
    acwpd(x, wt[, L])

Performs a discrete autocorrelation wavelet packet transform for a given signal `x`. The
wavelet type `wt` determines the transform type. Refer to Wavelet.jl for a list of available
methods.

# Arguments
- `x::AbstractArray{T} where T<:Number`: Original signal, preferably of size 2ᴷ or (2ᴷ,2ᴹ)
  where ``K, M \in \mathbb{N}``.
- `wt::OrthoFilter`: Orthogonal wavelet filter.
- `L::Integer`: (Default: `maxtransformlevels(x)`) Number of levels of decomposition.

# Returns
`::Array{T}`: Output from ACWPD on `x`.

# Examples
```julia
using Wavelets, WaveletsExt

# Setup
x = generatesignals(:heavysine)
wt = wavelet(WT.haar)

# ACWPD
acwpd(x, wt)

acwpd(x, wt, 4)
```

**See also:** [`iacwpd`](@ref), [`acdwt`](@ref), [`acdwt_step`](@ref)
"""
function acwpd(x::AbstractArray{T},
               wt::OrthoFilter,
               L::Integer = maxtransformlevels(x)) where T<:Number
    # Sanity check
    @assert 1 ≤ ndims(x) ≤ 2
    @assert L ≤ maxtransformlevels(x) ||
        throw(ArgumentError("Too many transform levels (length(x) < 2^L"))
    @assert L ≥ 1 || throw(ArgumentError("L must be ≥ 1"))
    # Setup
    sz = size(x)                            # Signal size
    N = ndims(x)+1
    k = N==2 ? 1<<(L+1)-1 : sum(4 .^(0:L))  # Total number of nodes
    xw = Array{T}(undef, (sz...,k))
    # Transform
    acwpd!(xw, x, wt, L)
    return xw
end

# 1D ACWPD
@doc raw"""
    acwpd!(xw, x, wt[, L])

Same as `acwpd` but without array allocation.

# Arguments
- `xw::AbstractArray{T} where T<:Number`: Allocated array for output.
- `x::AbstractArray{T} where T<:Number`: Original signal, preferably of size 2ᴷ where ``K
  \in \mathbb{N}``.
- `wt::OrthoFilter`: Orthogonal wavelet filter.
- `L::Integer`: (Default: `maxtransformlevels(x)`) Number of levels of decomposition.

# Returns
`xw::Array{T}`: Output from ACWPD on `x`.

# Examples
```julia
using Wavelets, WaveletsExt

# Setup
x = generatesignals(:heavysine)
wt = wavelet(WT.haar)

# ACWPD
xw = Matrix{Float64}(undef, (128, 255))
acwpd!(xw, x, wt)
acwpd!(xw, x, wt, 7)
```

**See also:** [`acwpd!`](@ref)
"""
function acwpd!(xw::AbstractArray{T,2},
                x::AbstractVector{T},
                wt::OrthoFilter,
                L::Integer=maxtransformlevels(x)) where T<:Number
    # Sanity check
    @assert L <= maxtransformlevels(x) ||
        throw(ArgumentError("Too many transform levels (length(x) < 2^L"))
    @assert L >= 1 || throw(ArgumentError("L must be >= 1"))

    # Setup
    Pmf, Qmf = make_acreverseqmfpair(wt)
    n₀ = 1<<(L+1)-1                     # Total number of nodes
    n₁ = n₀ - (1<<L)                    # Total number of nodes excluding leaf nodes
    xw[:,1] = x

    # ACWPD
    for i in 1:n₁
        d = floor(Int, log2(i))
        j₁ = getchildindex(i,:left)
        j₂ = getchildindex(i,:right)
        @inbounds v = @view xw[:,i]
        @inbounds w₁ = @view xw[:,j₁]
        @inbounds w₂ = @view xw[:,j₂]
        @inbounds acdwt_step!(w₁, w₂, v, d, Qmf, Pmf)
    end
    return xw
end
# 2D ACWPD
function acwpd!(xw::AbstractArray{T,3},
                x::AbstractMatrix{T},
                wt::OrthoFilter,
                L::Integer = maxtransformlevels(x)) where T<:Number
    # Sanity check
    n, m, k = size(xw)
    @assert k == sum(4 .^(0:L))
    @assert n == size(x,1)
    @assert m == size(x,2)

    # Setup
    g, h = make_acreverseqmfpair(wt)
    temp = Array{T,3}(undef, (n,m,2))
    n₁ = k - (4^L)                    # Total number of nodes excluding leaf nodes
    for i in eachindex(x)
        @inbounds xw[i] = x[i]
    end

    # ACWPD
    for i in 1:n₁
        d = getdepth(i,:quad)
        j₁ = getchildindex(i,:topleft)
        j₂ = getchildindex(i,:topright)
        j₃ = getchildindex(i,:bottomleft)
        j₄ = getchildindex(i,:bottomright)
        v = @view xw[:,:,i]
        w₁ = @view xw[:,:,j₁]
        w₂ = @view xw[:,:,j₂]
        w₃ = @view xw[:,:,j₃]
        w₄ = @view xw[:,:,j₄]
        acdwt_step!(w₁, w₂, w₃, w₄, v, d, h, g, temp)
    end
end

"""
    iacwpd(xw, L)
    iacwpd(xw[, wt, L])
    iacwpd(xw, tree)
    iacwpd(xw, wt, tree)

Performs the inverse autocorrelation discrete wavelet packet transform, with respect to a
decomposition tree.

!!! note
    The inverse autocorrelation transform does not require any wavelet filter, but an
    optional `wt` positional argument is included for the standardization of syntax with
    `wpt` and `swpt`, but is ignored during the reconstruction of signals.

!!! note
    This function might not be very useful if one is looking to reconstruct a raw
    decomposed signal. The purpose of this function would be better utilized in applications
    such as denoising, where a signal is decomposed (`swpd`) and thresholded
    (`denoise`/`denoiseall`) before being reconstructed.

# Arguments
- `xw::AbstractArray{T} where T<:Number`: ACWPD-transformed array.
- `wt::Union{OrthoFilter, Nothing}`: (Default: `nothing`) Orthogonal wavelet filter.
- `L::Integer`: (Default: `maxtransformlevels(x)`) Number of levels of decomposition used
  for reconstruction.
- `tree::BitVector`: Binary tree for inverse transform to be computed accordingly.

# Returns
`::Array{T}`: Inverse transformed signal.

# Examples
```julia
using Wavelets, WaveletsExt

# Setup
x = generatesignals(:heavysine)
wt = wavelet(WT.haar)

# ACWPD
xw = acwpd(x, wt)

# IACWPD
x̂ = iacwpd(xw, 4)
x̂ = iacwpd(xw, wt, 4)
x̂ = iacwpd(xw, maketree(x))
x̂ = iacwpd(xw, wt, maketree(x))
```

**See also:** [`acwpd`](@ref)
"""
function iacwpd(xw::AbstractArray{T},
                wt::Union{OrthoFilter, Nothing} = nothing,
                L::Integer = maxtransformlevels(xw,1)) where T<:Number
    @assert 2 ≤ ndims(xw) ≤ 3
    return iacwpd(xw, L)
end

function iacwpd(xw::AbstractArray{T}, L::Integer) where T<:Number
    # Sanity check
    @assert 2 ≤ ndims(xw) ≤ 3
    @assert L ≤ maxtransformlevels(xw,1) ||
            throw(ArgumentError("Too many transform levels (length(x) < 2ᴸ)"))
    @assert L ≥ 1 || throw(ArgumentError("L must be ≥ 1"))
    sz = size(xw)[1:(end-1)]
    return iacwpd(xw, maketree(sz..., L, :full))
end

function iacwpd(xw::AbstractArray{T},
                wt::Union{OrthoFilter, Nothing},
                tree::BitVector) where T<:Number
    @assert 2 ≤ ndims(xw) ≤ 3
    return iacwpd(xw, tree)
end

function iacwpd(xw::AbstractArray{T}, tree::BitVector) where T<:Number
    @assert 2 ≤ ndims(xw) ≤ 3
    sz = size(xw)[1:(end-1)]
    x = Array{T}(undef, sz)
    iacwpd!(x, xw, tree)
    return x
end

"""
    iacwpd!(x, xw, L)
    iacwpd!(x, xw[, wt, L])
    iacwpd!(x, xw, tree)
    iacwpd!(x, xw, wt, tree)

Same as `iacwpd` but with no array allocation.

# Arguments
- `x::AbstractArray{T} where T<:Number`: Allocated array for output.
- `xw::AbstractArray{T} where T<:Number`: ACWPD-transformed array.
- `wt::OrthoFilter`: Orthogonal wavelet filter.
- `L::Integer`: (Default: `maxtransformlevels(x)`) Number of levels of decomposition used
  for reconstruction.
- `tree::BitVector`: Binary tree for inverse transform to be computed accordingly.

# Returns
`x::Array{T}`: Inverse transformed signal.

# Examples
```julia
using Wavelets, WaveletsExt

# Setup
x = generatesignals(:heavysine)
wt = wavelet(WT.haar)

# ACWPD
xw = acwpd(x, wt)

# IACWPD
x̂ = similar(x)
iacwpd!(x̂, xw, 4)
iacwpd!(x̂, xw, wt, 4)
iacwpd!(x̂, xw, maketree(x))
iacwpd!(x̂, xw, wt, maketree(x))
```

**See also:** [`iacwpd`](@ref)
"""
function iacwpd!(x::AbstractArray{T},
                 xw::AbstractArray{T},
                 wt::Union{OrthoFilter, Nothing} = nothing,
                 L::Integer = maxtransformlevels(x)) where T<:Number
    iacwpd!(x, xw, L)
    return x
end

function iacwpd!(x::AbstractArray{T}, xw::AbstractArray{T}, L::Integer) where
                 T<:Number
    @assert ndims(x) == ndims(xw)-1
    @assert size(x) == size(xw)[1:(end-1)]
    @assert L ≤ maxtransformlevels(x) ||
            throw(ArgumentError("Too many transform levels."))
    @assert L ≥ 1 || throw(ArgumentError("L must be ≥ 1."))
    iacwpd!(x, xw, maketree(size(x)..., L, :full))
    return x
end

function iacwpd!(x::AbstractArray{T},
                 xw::AbstractArray{T},
                 wt::Union{OrthoFilter, Nothing},
                 tree::BitVector) where T<:Number
    iacwpd!(x, xw, tree)
    return x
end
# 1D IACWPD
function iacwpd!(x::AbstractVector{T}, xw::AbstractMatrix{T}, tree::BitVector) where
                 T<:Number
    # Sanity check
    @assert size(x,1) == size(xw,1)
    @assert isvalidtree(x, tree)

    # Setup
    tmp = copy(xw)

    # iACWPD
    for (i, haschild) in Iterators.reverse(enumerate(tree))
        # Current node has child => Compute one step of inverse transform
        if haschild
            d = floor(Int, log2(i))                 # Parent depth
            j₁ = getchildindex(i,:left)             # Scaling child index
            j₂ = getchildindex(i,:right)            # Detail child index
            @inbounds v = i==1 ? x : @view tmp[:,i] # Parent node
            @inbounds w₁ = @view tmp[:,j₁]          # Scaling child node
            @inbounds w₂ = @view tmp[:,j₂]          # Detail child node
            # Inverse transform
            @inbounds iacdwt_step!(v, w₁, w₂)
        end
    end
    return x
end
# 2D IACWPD
function iacwpd!(x::AbstractMatrix{T}, xw::AbstractArray{T,3}, tree::BitVector) where
                 T<:Number
    # Sanity check
    @assert size(x) == size(xw)[1:(end-1)]
    @assert isvalidtree(x, tree)

    # Setup
    n, m, k = size(xw)
    xwₜ = copy(xw)                           # Temp. array to store intermediate outputs
    temp = Array{T,3}(undef, (n,m,2))

    # IACWPD
    for (i, haschild) in Iterators.reverse(enumerate(tree))
        # Current node has child => Compute one step of inverse transform
        if haschild
            d = getdepth(i,:quad)                       # Parent depth
            j₁ = getchildindex(i,:topleft)              # Scaling child index
            j₂ = getchildindex(i,:topright)             # Detail child index
            j₃ = getchildindex(i,:bottomleft)           # Detail child index
            j₄ = getchildindex(i,:bottomright)          # Detail child index
            @inbounds v = i==1 ? x : @view xwₜ[:,:,i]    # Parent node
            @inbounds w₁ = @view xwₜ[:,:,j₁]             # Scaling child node
            @inbounds w₂ = @view xwₜ[:,:,j₂]             # Detail child node
            @inbounds w₃ = @view xwₜ[:,:,j₃]             # Detail child node
            @inbounds w₄ = @view xwₜ[:,:,j₄]             # Detail child node
            # Inverse transform
            @inbounds iacdwt_step!(v, w₁, w₂, w₃, w₄, temp)
        end
    end
    return x
end

# function iacwpt(xw::AbstractArray{<:Number,2}, tree::BitVector)
#   @assert isvalidtree(xw[:,1], tree)
#   v₀ = iacwpt(xw,tree,2)
#   v₁ = iacwpt(xw,tree,3)
#   return (v₀ + v₁) / √2
# end

include("acwt_utils.jl")
include("acwt_one_level.jl")
include("acwt_all.jl")

end # end module


text:
export
SumProductRule,
@sumProductRule

"""
A non-specific sum-product update
"""
abstract type SumProductRule{factor_type} <: MessageUpdateRule end

"""
`internalSumProductSchedule()` generates a sum-product message passing schedule
for the inner graph of a `CompositeFactor`. This schedule produces the sum-product
message out of the specified `outbound_interface`.
"""
function internalSumProductSchedule(cnode::CompositeFactor,
                                    outbound_interface::Interface,
                                    inferred_outbound_types::Dict{Interface, <:Type})

    # Collect types of messages towards the CompositeFactor
    msg_types = Dict{Interface, Type}()
    for (idx, terminal) in enumerate(cnode.terminals)
        (cnode.interfaces[idx] === outbound_interface) && continue # don't need incoming msg on outbound interface
        msg_types[terminal.interfaces[1]] = inferred_outbound_types[cnode.interfaces[idx].partner]
    end

    terminal_interfaces = [term.interfaces[1] for term in cnode.terminals]
    target_terminal = cnode.interface2terminal[outbound_interface]
    internal_outbound_interface = target_terminal.interfaces[1].partner

    # Generate a feasible summary propagation schedule
    schedule = summaryPropagationSchedule(  Variable[],
                                            limit_set=edges(cnode.inner_graph),
                                            target_sites=[internal_outbound_interface],
                                            breaker_sites=terminal_interfaces)

    # Assign the sum-product update rule to each of the schedule entries
    for entry in schedule
        entry.message_update_rule = SumProductRule{typeof(entry.interface.node)}
    end

    # Sanity check
    if schedule[end].interface !== internal_outbound_interface
        error("The last schedule entry should correspond to the message going out of the CompositeFactor")
    end

    # Type inference for internal messages
    inferUpdateRules!(schedule, inferred_outbound_types=msg_types)

    return schedule
end

function inferUpdateRule!(entry::ScheduleEntry,
                          rule_type::Type{T},
                          inferred_outbound_types::Dict{Interface, <:Type}
                         ) where T<:SumProductRule
    # Collect inbound types
    inbound_types = collectInboundTypes(entry, rule_type, inferred_outbound_types)

    # Find applicable rule(s)
    applicable_rules = Type[]
    for rule in leaftypes(entry.message_update_rule)
        if isApplicable(rule, inbound_types)
            push!(applicable_rules, rule)
        end
    end

    # Select and set applicable rule
    if isempty(applicable_rules)
        if isa(entry.interface.node, CompositeFactor)
            # No 'shortcut rule' available for CompositeFactor.
            # Try to fall back to msg passing on the internal graph.
            entry.internal_schedule = internalSumProductSchedule(entry.interface.node, entry.interface, inferred_outbound_types)
            entry.message_update_rule = entry.internal_schedule[end].message_update_rule
        else
            error("No applicable $(rule_type) update for $(typeof(entry.interface.node)) node with inbound types: $(join(inbound_types, ", "))")
        end
    elseif length(applicable_rules) > 1
        error("Multiple applicable $(rule_type) updates for $(typeof(entry.interface.node)) node with inbound types: $(join(inbound_types, ", ")): $(join(applicable_rules, ", "))")
    else
        entry.message_update_rule = first(applicable_rules)
    end

    return entry
end

function collectInboundTypes(entry::ScheduleEntry,
                             ::Type{T},
                             inferred_outbound_types::Dict{Interface, <:Type}
                            ) where T<:SumProductRule
    inbound_message_types = Type[]
    for node_interface in entry.interface.node.interfaces
        if node_interface === entry.interface
            push!(inbound_message_types, Nothing)
        elseif isDeltaConstrained(node_interface.partner)
            push!(inbound_message_types, Message{PointMass})
        else
            push!(inbound_message_types, inferred_outbound_types[node_interface.partner])
        end
    end

    return inbound_message_types
end

"""
`@sumProductRule` registers a sum-product update rule by defining the rule type
and the corresponding methods for the `outboundType` and `isApplicable`
functions. If no name (type) for the new rule is passed, a unique name (type)
will be generated. Returns the rule type.
"""
macro sumProductRule(fields...)
    # Init required fields in macro scope
    node_type = :unknown
    outbound_type = :unknown
    inbound_types = :unknown
    name = :auto # Triggers automatic naming unless overwritten

    # Loop over fields because order is unknown
    for arg in fields

        (arg.args[1] == :(=>)) || error("Invalid call to @sumProductRule")

        if arg.args[2].value == :node_type
            node_type = arg.args[3]
        elseif arg.args[2].value == :outbound_type
            outbound_type = arg.args[3]
            (outbound_type.head == :curly && outbound_type.args[1] == :Message) || error("Outbound type for SumProductRule should be a Message")
        elseif arg.args[2].value == :inbound_types
            inbound_types = arg.args[3]
            (inbound_types.head == :tuple) || error("Inbound types should be passed as Tuple")
        elseif arg.args[2].value == :name
            name = arg.args[3]
        else
            error("Unrecognized field $(arg.args[2].value) in call to @sumProductRule")
        end
    end

    # Assign unique name if not set already
    if name == :auto
        # Added hash ensures that the rule name is unique
        msg_types_hash = string(hash(vcat([outbound_type], inbound_types)))[1:6]
        name = Symbol("SP$(node_type)$(msg_types_hash)")
    end

    # Build validators for isApplicable
    input_type_validators = Expr[]

    push!(input_type_validators, :(length(input_types) == $(length(inbound_types.args))))
    for (i, i_type) in enumerate(inbound_types.args)
        if i_type != :Nothing
            # Only validate inbounds required for message update
            push!(input_type_validators, :(ForneyLab.matches(input_types[$i], $i_type)))
        end
    end

    expr = quote
        struct $name <: SumProductRule{$node_type} end
        ForneyLab.outboundType(::Type{$name}) = $outbound_type
        ForneyLab.isApplicable(::Type{$name}, input_types::Vector{<:Type}) = begin
            $(reduce((current, item) -> :($current && $item), input_type_validators, init = :true))
        end
    end

    return esc(expr)
end

"""
Collect and construct SP update code for each inbound.
"""
collectInbounds(entry::ScheduleEntry, ::Type{T}) where T<:SumProductRule = collectSumProductNodeInbounds(entry.interface.node, entry)

"""
Construct the inbound code that computes the message for `entry`. Allows for
overloading and for a user the define custom node-specific inbounds collection.
Returns a vector with inbounds that correspond with required interfaces.
"""
function collectSumProductNodeInbounds(::FactorNode, entry::ScheduleEntry)
    interface_to_schedule_entry = current_inference_algorithm.interface_to_schedule_entry

    inbounds = Any[]
    for node_interface in entry.interface.node.interfaces
        inbound_interface = ultimatePartner(node_interface)
        if node_interface === entry.interface
            # Ignore inbound message on outbound interface
            push!(inbounds, nothing)
        elseif isClamped(inbound_interface)
            # Hard-code outbound message of constant node in schedule
            push!(inbounds, assembleClamp!(inbound_interface.node, Message))
        else
            # Collect message from previous result
            push!(inbounds, interface_to_schedule_entry[inbound_interface])
        end
    end

    return inbounds
end


text:
using ExtXYZ
using Test
using PyCall

"""
Check two dictionaries `seq1` and `seq2` are approximately equal

Returns true if all keys match, all non-float values match, and all pairs
of float values `(v1, v2)` satisfy `isapprox(v1, v2)`
"""
function Base.isapprox(seq1::AbstractDict, seq2::AbstractDict)
    for (k1,v1) in seq1
        k1 ∈ keys(seq2) || (println("key $k1 missing from seq2"); return false)
        if v1 isa AbstractDict
            return isapprox(v1, seq2[k1])
        elseif v1 isa Array{AbstractFloat} || v1 isa AbstractFloat
            v1 ≈ seq2[k1]  || (println("key $k1: $v1 !≈ $(seq2[k1])"); return false)
        else
            v1 == seq2[k1] || (println("key $k1: $v1 != $(seq2[k1])"); return false)
        end
    end
    return true
end

@testset "ExtXYZ.jl" begin
    test_xyz(frame) = """2
cutoff=5.50000000 pbc=[T, T, T] nneightol=1.20000000 config_type=isolated_atom matrix=[[$frame.00000000, 1.00000000, 2.00000000], [3.00000000, 4.00000000, 5.00000000], [6.00000000, 7.00000000, 8.00000000]] dft_energy=-158.$(frame)4496821 Lattice="1.00000000 2.00000000 3.00000000 4.00000000 5.00000000 6.00000000 7.00000000 $(frame).00000000 9.00000000" gap_energy=-157.72725320  Properties=species:S:1:pos:R:3:n_neighb:I:1:dft_force:R:3:gap_force:R:3:map_shift:I:3
Si        10.00000000      11.00000000      $frame.00000000          0         0.10000000       0.20000000       0.30000000         0.1$(frame)000000       0.22000000       0.330000000         -1       -2       0
Si        13.00000000      14.00000000      $(frame+1).00000000          0         0.10000000       0.20000000       0.30000000         0.2$(frame)000000       0.22000000       0.330000000         -1       -2       0
"""

    infile = "test.xyz"
    open(infile, "w") do io
        for frame=1:10
            print(io, test_xyz(frame))
        end
    end
    outfile = "dump.xyz"

    try
        @testset "read" begin
            seq1 = read_frames(infile)

            @test all([s["cell"][3,2] for s in seq1] .== 1.0:10.0)
            @test all([s["arrays"]["pos"][3,1] for s in seq1] .== 1.0:10.0)
            @test all([s["arrays"]["pos"][3,2] for s in seq1] .== 2.0:11.0)
            @test all([s["info"]["matrix"][1,1] for s in seq1] .== 1.0:10.0)
            @test all([s["info"]["matrix"][2,1] for s in seq1] .== 3.0)

            @test all([s["info"]["dft_energy"] for s in seq1] .== [parse(Float64, "-158.$(frame)4496821") for frame=1:10])
            @test all([s["arrays"]["gap_force"][1,1] for s in seq1] .== [parse(Float64, "0.1$(frame)") for frame=1:10])
            @test all([s["arrays"]["gap_force"][1,2] for s in seq1] .== [parse(Float64, "0.2$(frame)") for frame=1:10])

            seq2 = read_frames(infile, 4:10)
            @test all(seq1[4:10] .== seq2)

            frame4 = read_frame(infile, 4)
            @test frame4 == seq1[4]

            frame1 = read_frame(infile)
            @test frame1 == seq1[1]

            f = open(infile, "r")
            seq3 = read_frames(f)
            close(f)

            @test all(seq1 .== seq3)
        end

        @testset "convert" begin
            nat, info, arrays = ExtXYZ.cfopen(infile, "r") do fp
                ExtXYZ.read_frame_dicts(fp)
            end
            cinfo = convert(Ptr{ExtXYZ.DictEntry}, info)
            carrays = convert(Ptr{ExtXYZ.DictEntry}, arrays)

            ninfo = convert(Dict{String,Any}, cinfo)
            narrays = convert(Dict{String,Any}, carrays)

            @test ninfo == info
            @test narrays == arrays
        end

        @testset "write" begin
            seq1 = read_frames(infile)
            write_frames(outfile, seq1)
            seq2 = read_frames(outfile)
            @test all(seq1 .≈ seq2) # use custom isapprox(), since expect some loss of precision on round-trip
        end

        @testset "iread" begin
            seq1 = read_frames(infile)
            seq2 = collect(iread_frames(infile))
            @test all(seq1 .≈ seq2)
        end

        @testset "iwrite" begin
            seq1 = read_frames(infile)
            ch = Channel()
            job = @async write_frames(outfile, ch)
            for frame in seq1
                put!(ch, frame)
            end
            close(ch)
            wait(job)
            seq2 = read_frames(outfile)
            @test all(seq1 .≈ seq2)
        end

        try
            ase_io = pyimport("ase.io")

            @testset "ASE" begin
                seq = read_frames(infile)
                ase_seq = ase_io.read(infile * "@:")

                for (frame, ase_atoms) in zip(seq, ase_seq)
                    @test frame["N_atoms"] ≈ length(ase_atoms)
                    @test frame["arrays"]["pos"] ≈ ase_atoms.positions'
                    @test frame["arrays"]["gap_force"] ≈ ase_atoms.arrays["gap_force"]'
                    @test frame["arrays"]["n_neighb"] ≈ ase_atoms.arrays["n_neighb"]
                    @test frame["cell"] ≈ ase_atoms.cell.array
                end
            end
        catch e
            if e isa PyError
                println("ASE not installed, skipping ASE comparison tests")
            else
                throw(e)
            end
        end

    finally
        rm(infile, force=true)
        rm(outfile, force=true)
    end
end


text:
methods(+)

@which 3 + 3

@which 3.0 + 3.0

@which 3 + 3.0

import Base: +

"hello " + "world!"

@which "hello " + "world!"

+(x::String, y::String) = string(x, y)

"hello " + "world!"

@which "hello " + "world!"

foo(x, y) = println("duck-typed foo!")
foo(x::Int, y::Float64) = println("foo con entero y flotante!")
foo(x::Float64, y::Float64) = println("foo con dos flotantes!")
foo(x::Int, y::Int) = println("foo con dos enteros!")

foo(1, 1)

foo(1., 1.)

foo(1, 1.0)

foo(true, false)

max_stars_repo_name: miguelraz/JuliaTutorials

max_issues_repo_name: miguelraz/JuliaTutorials

max_forks_repo_name: miguelraz/JuliaTutorials


text:
__precompile__()

module AutoUrban

using AutomotiveSimulator
using AutomotiveVisualization
using EzXML
using Printf
import Records: QueueRecord, Frame
using LinearAlgebra
using Distributions
import Cairo

export
	LatLonAccelDirection,
	AccelSteeringDirection,
    NextState

include("actions.jl")

export
    initialize_XML,
    get_lane_length,
    handle_junctions,
    convert_roadway!,
    Connection,
    Junction,
    connect_two_lane_general!,
    connect_two_seg!,
    connect_two_seg_general!,
    add_connection!,
    add_junction!,
    gen_connected_lanes,
    gen_intersection,
    gen_loop_roadway,
    add_line!,
    in_lanes,
    get_max_curvature

include("roadway/roadway_generation.jl")
include("roadway/roadway_perception.jl")
include("roadway/roadway_intersection.jl")
include("roadway/convert2xodr.jl")

export
    AccSteerDriver,
    IDMDriver,
    UrbanDriver,
    MultiPtsDriver,
    MultiPtsTurningDriver,
    MOBILDriver,
    excute_action!

include("drivermodels/acc_steer_drivers.jl")
include("drivermodels/IDMDriver.jl")
include("drivermodels/MOBILDriver.jl")
include("drivermodels/urban_drivers.jl")
include("drivermodels/multiPts_driver.jl")
include("drivermodels/multiPts_functions_GD.jl")
include("drivermodels/waypoints_setting_functions.jl")
include("drivermodels/multiPtsTurning_driver.jl")
include("drivermodels/render_drivers.jl")

end

max_stars_repo_name: sisl/UrbanDrivingSimulation

max_issues_repo_name: sisl/UrbanDrivingSimulation

max_forks_repo_name: sisl/UrbanDrivingSimulation


text:
using Plots

function lineplot(st::SearchTrace, xdata::Symbol, ydata::Symbol; kwargs...)
    return plot(st.df[:, xdata], st.df[:, ydata]; kwargs...)
end

function objective_plot(st::SearchTrace, xsym::Symbol=:it, ysym::Symbol=:v_s; kwargs...)
    return lineplot(st, xsym, ysym; kwargs...)
end

max_stars_repo_name: Michiel-VL/PolyaViz

max_issues_repo_name: Michiel-VL/PolyaViz

max_forks_repo_name: Michiel-VL/PolyaViz


text:
using Documenter, HermitianDiagonalization

makedocs(
    modules = [HermitianDiagonalization],
    format = :html,
    checkdocs = :exports,
    sitename = "HermitianDiagonalization.jl",
    pages = Any["index.md"]
)

deploydocs(
    repo = "github.com/pablosanjose/HermitianDiagonalization.jl.git",
)

max_stars_repo_name: pablosanjose/HermitianDiagonalization

max_issues_repo_name: pablosanjose/HermitianDiagonalization

max_forks_repo_name: pablosanjose/HermitianDiagonalization


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

Time domain filtering of 1-D data tables.

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

Parameters
----------

- **A** : **single_vec** : -- Str --   Flags = m[conf]|vector

    Specify a single, primary vector instead of reading tables.
    [`-A`](http://gmt.soest.hawaii.edu/doc/latest/gmtvector.html#a)
- **C** : **cartesian** : -- Str or [] --        Flags = [i|o]

    Select Cartesian coordinates on input and output.
    [`-C`](http://gmt.soest.hawaii.edu/doc/latest/gmtvector.html#c)
- **E** : **geod2geoc** : -- Bool or [] --

    Convert input geographic coordinates from geodetic to geocentric and output geographic
    coordinates from geocentric to geodetic.
    [`-E`](http://gmt.soest.hawaii.edu/doc/latest/gmtvector.html#e)
- **N** : **normalize** : -- Bool or [] --

    Normalize the resultant vectors prior to reporting the output.
    [`-N`](http://gmt.soest.hawaii.edu/doc/latest/gmtvector.html#n)
- **S** : **secondary_vec** : -- Str or List --    Flags = [vector]

    Specify a single, secondary vector in the same format as the first vector.
    [`-S`](http://gmt.soest.hawaii.edu/doc/latest/gmtvector.html#s)
- **T** : **transform** : -- List or Str --     Flags = a|d|D|paz|s|r[arg|R|x]

    Specify the vector transformation of interest.
    [`-T`](http://gmt.soest.hawaii.edu/doc/latest/gmtvector.html#t)
- $(GMT.opt_V)
- $(GMT.opt_write)
- $(GMT.opt_append)
- $(GMT.opt_b)
- $(GMT.opt_d)
- $(GMT.opt_e)
- $(GMT.opt_f)
- $(GMT.opt_g)
- $(GMT.opt_h)
- $(GMT.opt_i)
- $(GMT.opt_o)
- $(GMT.opt_swap_xy)
"""
function gmtvector(cmd0::String="", arg1=[]; kwargs...)

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

	d = KW(kwargs)

	cmd = parse_V_params("", d)
	cmd, = parse_b(cmd, d)
	cmd, = parse_d(cmd, d)
	cmd, = parse_e(cmd, d)
	cmd, = parse_f(cmd, d)
	cmd, = parse_g(cmd, d)
	cmd, = parse_h(cmd, d)
	cmd, = parse_i(cmd, d)
	cmd, = parse_o(cmd, d)
	cmd, = parse_swap_xy(cmd, d)

	cmd = add_opt(cmd, 'A', d, [:A :single_vec])
	cmd = add_opt(cmd, 'C', d, [:C :cartesian])
	cmd = add_opt(cmd, 'E', d, [:E :geod2geoc])
	cmd = add_opt(cmd, 'N', d, [:N :normalize])
	cmd = add_opt(cmd, 'S', d, [:S :secondary_vec])
	cmd = add_opt(cmd, 'T', d, [:T :transform])

	cmd, got_fname, arg1 = find_data(d, cmd0, cmd, 1, arg1)
	return common_grd(d, cmd, got_fname, 1, "gmtvector", arg1)		# Finish build cmd and run it
end

# ---------------------------------------------------------------------------------------------------
gmtvector(arg1=[], cmd0::String=""; kw...) = gmtvector(cmd0, arg1; kw...)


text:
"""Implementation of (simplified) PRM

ref:
- Kavraki, L. E., Svestka, P., Latombe, J. C., & Overmars, M. H. (1996).
  Probabilistic roadmaps for path planning in high-dimensional configuration spaces.
  IEEE transactions on Robotics and Automation
- Karaman, S., & Frazzoli, E. (2011).
  Sampling-based algorithms for optimal motion planning.
  The international journal of robotics research (IJRR)
"""

"""
    PRM_direct(
        config_init::Vector{State},
        config_goal::Vector{State},
        connect::Function,
        collide::Function,
        check_goal::Function;
        epsilon::Union{Float64,Nothing} = 0.2,
        TIME_LIMIT::Union{Nothing,Real} = nothing,
        VERBOSE::Int64 = 0,
    )::Tuple{
        Union{Nothing,Vector{Vector{Node{State}}}},  # solution
        Vector{Vector{Node{State}}},  # roadmaps
    } where {State<:AbsState}

solve the problem directly by PRM
"""
function PRM_direct(
    config_init::Vector{State},
    config_goal::Vector{State},
    connect::Function,
    collide::Function,
    check_goal::Function;
    epsilon::Union{Float64,Nothing} = 0.2,
    TIME_LIMIT::Union{Nothing,Real} = nothing,
    VERBOSE::Int64 = 0,
)::Tuple{
    Union{Nothing,Vector{Vector{Node{State}}}},  # solution
    Vector{Vector{Node{State}}},  # roadmaps
} where {State<:AbsState}

    t_s = now()
    elapsed() = elapsed_sec(t_s)
    timeover() = TIME_LIMIT != nothing && elapsed() > TIME_LIMIT

    N = length(config_init)

    # define utilities
    _sampler = MRMP.gen_uniform_sampling(config_init[1])
    sampler() = map(i -> _sampler(), 1:N)
    connect_C(C_from::Vector{State}, C_to::Vector{State}) = begin
        (isnothing(epsilon) || dist(C_from, C_to) <= epsilon) &&
        all(i -> connect(C_from[i], C_to[i], i), 1:N) &&
        begin
            for i = 1:N, j = i+1:N
                if collide(C_from[i], C_to[i], C_from[j], C_to[j], i, j)
                    return false
                end
            end
            return true
        end
    end

    get_roadmaps() = begin
        roadmaps = map(i -> Vector{Node{State}}(), 1:N)
        for (k, C) in enumerate(V)
            for (i, q) in enumerate(C)
                push!(roadmaps[i], Node{State}(q, k, E[k]))
            end
        end
        roadmaps
    end

    # vertices
    V = [config_init, config_goal]

    # edges
    E = [[], []]

    # distance, used with Dijkstra method
    D = [0.0, nothing]

    # special case
    if connect_C(config_init, config_goal)
        return (
            [
                map(q -> Node{State}(q, 1, []), config_init),
                map(q -> Node{State}(q, 2, []), config_goal),
            ],
            get_roadmaps(),
        )
    end

    iter = 0
    while !timeover()
        iter += 1

        VERBOSE > 1 &&
            iter % 100 == 0 &&
            @printf("\r\t%6.4f sec, iteration: %02d", elapsed(), iter)

        # sampling
        C = sampler()
        while !timeover() && any(i -> !connect(C[i], i), 1:N)
            C = sampler()
        end
        timeover() && break

        # insert new vertex
        push!(V, C)
        push!(D, nothing)
        L = length(V)
        # update neighbors
        neigh_from_C = collect(filter(k -> connect_C(C, V[k]), 1:L-1))
        neigh_to_C = collect(filter(k -> connect_C(V[k], C), 1:L-1))
        push!(E, neigh_from_C)
        foreach(k -> push!(E[k], L), neigh_to_C)

        isempty(neigh_from_C) && isempty(neigh_to_C) && continue

        # Dijkstra method
        # setup open list
        OPEN = PriorityQueue{Int64,Float64}()
        enqueue!(OPEN, 1, D[1])

        # parent
        P = fill(0, L)

        while !isempty(OPEN)
            # pop
            k = dequeue!(OPEN)

            # check goal
            if check_goal(V[k])
                # backtrack
                solution = Vector{Vector{Node{State}}}()
                l = k
                while l != 0
                    pushfirst!(solution, map(q -> Node{State}(q, l, []), V[l]))
                    l = P[l]
                end
                VERBOSE > 1 && @printf("\r\t%6.4f sec, iteration: %02d", elapsed(), iter)
                VERBOSE > 0 && @info @sprintf("\t%6.4f sec: found solution\n", elapsed())
                return (solution, get_roadmaps())
            end

            # expand
            for l in E[k]
                d = D[k] + dist(V[k], V[l])
                if isnothing(D[l]) || d < D[l]
                    # update connection
                    D[l] = d
                    P[l] = k
                    haskey(OPEN, l) && delete!(OPEN, l)
                    enqueue!(OPEN, l, d)
                end
            end
        end
    end

    VERBOSE > 0 && @info @sprintf("\n\t%6.4f sec: failed to find solution\n", elapsed())
    return (nothing, get_roadmaps())
end

"""
    PRM!(
        roadmap::Vector{Node{State}},
        connect::Function,
        num_vertices::Int64 = 100,
        rad::Union{Nothing,Float64} = nothing;
        TIME_LIMIT::Union{Nothing,Real} = nothing,
    )::Vector{Node{State}} where {State<:AbsState}

update known probabilistic roadmap
"""
function PRM!(
    roadmap::Vector{Node{State}},
    connect::Function,
    num_vertices::Int64 = 100,
    rad::Union{Nothing,Float64} = nothing;
    TIME_LIMIT::Union{Nothing,Real} = nothing,
)::Vector{Node{State}} where {State<:AbsState}

    t_s = now()
    elapsed() = elapsed_sec(t_s)
    timeover() = TIME_LIMIT != nothing && elapsed() > TIME_LIMIT

    sampler = MRMP.gen_uniform_sampling(roadmap[1].q)

    # sampling vertices
    while length(roadmap) < num_vertices && !timeover()
        q = sampler()
        # exclude sample not in free space
        if connect(q)
            push!(roadmap, Node{State}(q, length(roadmap) + 1, []))
        end
    end

    # identify neighbors
    for j = 1:num_vertices, k = 1:num_vertices
        timeover() && break
        j == k && continue
        q_from = roadmap[j].q
        q_to = roadmap[k].q
        if (isnothing(rad) || dist(q_from, q_to) <= rad) && connect(q_from, q_to)
            push!(roadmap[j].neighbors, k)
        end
    end

    return roadmap
end

"""
    PRM(
        q_init::State,
        q_goal::State,
        args...;
        kwargs...,
    )::Vector{Node{State}} where {State<:AbsState}

generate one PRM
"""
function PRM(
    q_init::State,
    q_goal::State,
    args...;
    kwargs...,
)::Vector{Node{State}} where {State<:AbsState}
    return PRM!(
        [Node{State}(q_init, 1, []), Node{State}(q_goal, 2, [])],
        args...;
        kwargs...,
    )
end

"""
    PRMs(
        config_init::Vector{State},
        config_goal::Vector{State},
        connect::Function,
        num_vertices::Int64 = 100;
        rad::Union{Nothing,Float64} = nothing,
        rads::Union{Vector{Nothing},Vector{Float64}} = fill(rad, length(config_init)),
        TIME_LIMIT::Union{Nothing,Real} = nothing,
    )::Vector{Vector{Node{State}}} where {State<:AbsState}

generate multiple PRMs
"""
function PRMs(
    config_init::Vector{State},
    config_goal::Vector{State},
    connect::Function,
    num_vertices::Int64 = 100;
    rad::Union{Nothing,Float64} = nothing,
    rads::Union{Vector{Nothing},Vector{Float64}} = fill(rad, length(config_init)),
    TIME_LIMIT::Union{Nothing,Real} = nothing,
)::Vector{Vector{Node{State}}} where {State<:AbsState}

    t_s = now()
    elapsed() = elapsed_sec(t_s)

    N = length(config_init)
    return map(
        i -> PRM(
            config_init[i],
            config_goal[i],
            (args...) -> connect(args..., i),
            num_vertices,
            rads[i];
            TIME_LIMIT = (isnothing(TIME_LIMIT) ? nothing : TIME_LIMIT - elapsed()),
        ),
        1:N,
    )
end

"""
    PRMs!(
        roadmaps::Vector{Vector{Node{State}}},
        connect::Function,
        num_vertices::Int64 = 100;
        rad::Union{Nothing,Float64} = nothing,
        rads::Union{Vector{Nothing},Vector{Float64}} = fill(rad, length(roadmaps)),
        TIME_LIMIT::Union{Nothing,Real} = nothing,
    )::Vector{Vector{Node{State}}} where {State<:AbsState}

update PRMs
"""
function PRMs!(
    roadmaps::Vector{Vector{Node{State}}},
    connect::Function,
    num_vertices::Int64 = 100;
    rad::Union{Nothing,Float64} = nothing,
    rads::Union{Vector{Nothing},Vector{Float64}} = fill(rad, length(roadmaps)),
    TIME_LIMIT::Union{Nothing,Real} = nothing,
)::Vector{Vector{Node{State}}} where {State<:AbsState}

    t_s = now()
    elapsed() = elapsed_sec(t_s)

    N = length(roadmaps)
    return map(
        i -> PRM!(
            roadmaps[i],
            (args...) -> connect(args..., i),
            num_vertices,
            rads[i];
            TIME_LIMIT = (isnothing(TIME_LIMIT) ? nothing : TIME_LIMIT - elapsed()),
        ),
        1:N,
    )
end